The Vacua of Dipolar Cavity Quantum Electrodynamics

1)Interesting model and results, rich phase diagram (role of frustration, order-by disorder due to cavity, collective subradiant superradiant phases, etc..)

2) Solid numerical work

3)well written

1)physics in the photonic sector could be more discussed (see changes)

2)few key points in the text could be better explained (see changes)

3) figures are a bit hard to read (too much information), captions do not always help

This manuscript by Schuler et al discusses the ground-state phase diagram of a two-dimensional model of interacting dipoles (two-level systems) coupled to a cavity mode, relevant for CavityQED platforms in the ultra-strong coupling limit. The model can be seen as intermediate between Ising and Dicke models, and in this sense shares some similarities with the Rabi-Hubbard model of CavityQED.

The results, obtained with exact diagonalization (after truncating the photon Hilbert space) reveal a very rich phase diagram, in particular on the triangular lattice which is in my opinion one of the highlight of this work.

I think the results are relevant and worth to be published in SciPost Physics. At the same time I have a number of comments on the manuscript that I would kindly request the authors to consider. I append them below.

1)The authors discuss the broken symmetry phases in terms of dipolar order, but do not touch upon whether the photon field also breaks the symmetry and becomes coherent, i.e. <a>\neq 0 in the thermodynamic limit, as I would expect. Could the authors clarify this point?
2) One could then ask whether photonic and dipolar order always come hand by hand (as often in Dicke-like models) or whether one could have exotic situations where in one sector the symmetry is broken but in the other the order is still fluctuating. Is this scenario conceivable in the present model, for example in presence of frustration, i.e. on the triangular lattice? This would be pretty neat I think.
3)Following up on the points above, I notice that in terms of numerics the focus is mainly on photon number, which is ok but not really the relevant photonic "order parameter". Clearly for finite size systems the average photon field would be always zero but I wonder whether there could be interesting signatures of criticality in the low frequency retarded photon Green's function, which should be easy to compute with ED. In principle this should diverge a zero frequency at the critical point.
4)Few points in the manuscript could be expanded/clarified. For example the discussion on pag6 on "collective subradiant", it is not clear the connection with the g->infty limit. Similarly, the discussion below Eq.9 could be expanded and the connection between the effective Hamiltonian and the emergence of a super-radiant phase in the top right corner of the phase diagram.
5)Figures are a bit dense and the caption, although long and detailed does not always help. I could not find a discussion of the central panel in figures 2 and 4. In particular Figures 2 and 4 could be perhaps split in multiple figures?

1 - Detailed introduction to models and approximation.
2- Clear definition of the numerical methods used.
3- In-depth analysis of the phases and of their features.

1- Lack of scaling towards the thermodynamic limit.

In "The Vacua of Dipolar Cavity Quantum Electrodynamics", Schuler and coauthors numerically investigate the ground-state properties of a lattice of $N$ dipoles interacting with a single bosonic mode.
The authors take into account both the dipole-dipole interactions (specifically, they consider nearest-neighbors interactions) and the collective interaction of the dipoles with a photonic mode.

The article is extremely interesting from two points of view. On the one hand, as correctly pointed out by the authors, this model describes a cavity QED system in ultrastrong coupling. On the other, it can be seen as a generalization of an Ising model, where local and collective terms also compete with those of a bosonic field.

The article is well-written and deserves publications. I have, however, a major remarks which I think the authors should address.
Having answered these two points, I can strongly recommend the publication of this paper.

RMK) The authors provide in very few cases the scaling towards the thermodynamic limit. For a finite-size system, no phase transition can be observed. Since the authors numerically investigate the QED model, i.e., they always consider finite-size models, it is impossible to derive the phase diagram. I would argue that only by studying how the various quantities change in their scaling towards the thermodynamic limit $N \to \infty$, and how they change in proximity of the phase transition, the authors can safely claim the presence of criticality.

I have also these minor remarks:
1) Figures 2(a) and 4(a).
-In the left panels, it is difficult to distinguish between dotted and dashed lines.
-In Figure 4(a) the "3SL Superradiant" has the apostrophes in the wrong direction.
-The right panels are also very difficult to interpret. I understand what the authors want to do (the vertical axis is for $J$ and the horizontal one is for the quantities in legend) but maybe it can be made slightly more clear by adding the vertical axis labels. That of Figure 2(a) has also two arrows in the vertical axis.
2) In Appendic F, the authors claim "In contrast to phase transitions with an abrupt change in the behavior of the order parameter (in the thermodynamic limit), crossovers between two regimes of states with different physical properties show a smooth change in some of the observables.". This definition is somehow wrong. I would say that it is not the fact that there is "a smooth change in some of the observables" but "a smooth change in all of the observables". Indeed, the phase transition is a nonanalytical change of the ground state. As such, it can be witnessed by the discontinuity of some observables. Thus, finding a discontinuous observable (the order parameter) is sufficient to claim criticality. Vice versa, to claim that there is a crossover one must show that the wave function is continuous, and therefore that all the observables are continuous. Finding a smooth observable does not prove the presence of a crossover.

Moreover, I think that the overall discussion could benefit from addressing these two points (however, they are not necessary for the correct comprehension of the paper):
1) There has been recent controversy on the possibility of using the two-level approximation in the context of ultrastrongly light-matter coupled systems (and some of the authors are part of those article series). To cite some of the papers: * Phys. Rev. A 98, 053819 (2018) * Nat. Commun. 10, 499 (2019) *Nat. Phys. 15, 803–808 (2019) *arXiv:2005.06499. The authors begin the Section 2 by considering the two level approximation of the anharmonic dipoles. I think the author should justify this approximation given this recent (and still ongoing) debate.
2) Concerning the presence of a crossover between the paraelectric and the collective subradiant phase. In Fig. 2 (a) the boundary seems to move to the right (as discussed by the authors in App. F). In my opinion, however, there are two different questions which need to be answered: 1) Is this really a crossover? 2) What is the position of the crossover (or of the phase transition) for $N \to \infty$? For example, for $J_{i,j}=0$ several different techniques could be used to argue the presence of a crossover (or of a phase transition) between the paraelectric and collective superradiant phase. For example, bosonization of the dipoles via Holstein–Primakoff transformation, or explicitly using the permutational invariance of the spin model, to efficiently diagonalize the Hamiltonian for much larger values of $N$. Another way to understand this feature could be the semiclassical approximation of the bosonized approximation or a Gutzwiller mean-field study of the system.

1- Provide the scaling with $N$ to justify the claims that a crossovers/phase transitions is taking place.
2- Correct the figures.
3- Change the definition of the crossover in Appendix F
4- [FACULTATIVE] Comment on the use of the two-level approximation.
5- [FACULTATIVE] Discuss more in detail the paraelectric to the collective subradiant phase transition for $J=0$, in order to understand the presence of a crossover/phase transition and the position of the boundary.

1. Studies a relevant problem of the interaction between cavity-mediated and short range interactions.

2. Identifies several novel and unanticipated forms of order - collective subradiant phase, and three sublattice superradiant states.

3. Shows how some of this physics can be understood by an effective model in the strong matter-light coupling limit.

4. Provides a good analysis of finite size effects and fluctuations near the phase boundaries.

In light of SciPost acceptance criterion 6, "Provide (directly in appendices, or via links to external repositories) all reproducibility-enabling resources: explicit details of experimental protocols, datasets and processing methods, processed data and code snippets used to produce figures, etc.;", I note that while this theoretical paper does clearly describe the methods that are used, the authors do not provide a link to a code repository or the direct numerical results of the simulations.

The manuscript by Schuler et al. discusses the phase diagram of a model of dipoles interacting both through a global cavity-mediated interaction, and through local (screened) electrostatic interactions. Such a model has been derived by some of the authors as the multipolar gauge form of effective Hamiltonian in previous work. The key purpose of this paper is to explore the consequences of this Hamiltonian and the nature of its ground state through finite-size exact diagonalization.

The work finds novel phases arise due to the competition of the two interactions, particularly in the context of the frustrated triangular lattice. I believe the results on the appearance of three sublattice superradiant states from the competition of frustrated Ising interactions and long range interactions are significant and could be considered both groundbreaking and as potentially opening up a new area of research. The results are clearly presented, and the manuscript sets this work in context of most other relevant work in this field (see below). For these reasons, I believe the work should be published. There are a few minor issues the authors should consider before publication noted below.

1. On page 4, when commenting on the importance of including electrostatic screening for a consistent treatment of a cavity, it may be relevant to cite [Andras Vukics and Peter Domokos, Phys. Rev. A 86, 053807 (2012)], which made a similar point.

2. On page 6, above Eq. 4, there is a statement that the correlations functions are related to second order moments of the phase order parameters, while the equation written in (4) appear to relate them to first order moments of order parameters. This should be clarified.

3. On page 11, there is a comment about the apparent connection between the current results and supersolidity. In this context, there seems a far more immediate connection to make to experiments on cold atoms in cavities with Raman pumping, where there is competition between long-range cavity mediated interactions and short ranged interactions and hopping. e.g. [Klinder et al., Phys. Rev. Lett. 115, 230403 (2015); Landig et al, Nature 532, 476 (2016)]. While the nature of interactions there differs from the current problem, there may be relevant connections to draw.

4. On page 11, "has been constraint" should read "has been constrained"

5. On pages 15/16, in the text of section E and the caption of figure 8, there are references to the states (1,1-1) having peaks at the "edges of the hexagonal boundaries". From the figure, it appears this phrasing is referring to there being peaks at the corners of the hexagonal boundary. If correct, stating these are at the corners would be clearer.

6. Consider providing open data for the results of the numerical simulation, or access to a code repository to allow reproduction of these simulations.