Equilibrium and Disorder-induced behavior in Quantum Light-Matter Systems

We analyze equilibrium properties of coupled-doped cavities described by the Jaynes-Cummings- Hubbard Hamiltonian. In particular, we characterize the entanglement of the system in relation to the insulating-superfluid phase transition. We point out the existence of a crossover inside the superfluid phase of the system when the excitations change from polaritonic to purely photonic. Using an ensemble statistical approach for small systems and stochastic-mean-field theory for large systems we analyze static disorder of the characteristic parameters of the system and explore the ground state induced statistics. We report on a variety of glassy phases deriving from the hybrid statistics of the system. On-site strong disorder induces insulating behavior through two different mechanisms. For disorder in the light-matter detuning, low energy cavities dominate the statistics allowing the excitations to localize and bunch in such cavities. In the case of disorder in the light- matter coupling, sites with strong coupling between light and matter become very significant, which enhances the Mott-like insulating behavior. Inter-site (hopping) disorder induces fluidity and the dominant sites are strongly coupled to each other.


Introduction
Quantum phase transitions are a remarkable zero-temperature phenomenon driven by quantum fluctuations [1]. Such transitions have been studied in many-body quantum systems where each quantum phase can be unambiguously defined. However, recent results show evidence that interesting aspects and important traces of the physics of novel quantum phase transitions may already be observed in the limit of a very few interacting sites [2][3][4][5][6]. This is especially clear in hybrid light-matter systems, such as coupled electromagnetic cavities doped with two-level impurities, where a Mott-insulating-to-superfluid crossover has been predicted for as few as six or seven sites [2]. These systems feature a composite fermion-boson excitation in each site, hence the term hybrid, and quantities such as the variance in energy for each site have been used as markers for the transition between different phases [7]. However, entanglement, a unique quantum correlation with no classical analogue which has been related to fundamental features of quantum phase transitions [8], may be regarded as a more adequate order parameter [4,9]. In this work, we study the system entanglement to show how such quantum correlations relate to the behaviors the system may present.
One possible Hamiltonian describing doped and coupled cavities is the Jaynes-Cummings-Hubbard (JCH) model [2,10], which, in some limiting approximations, mimics the more typical and simpler Bose-Hubbard one. The similarities between both models have prompted the use of the latter as a basis for the analysis of quantum phase transitions in the former for both a large [11] and a very small number of sites [12]. However, the analogy to this simpler model ignores the internal structure of each site, which prevents one from exploring the increased complexity of the JCH system. The implementation of such systems has been proposed in different quantum optical setups such as planar lattices of 3 one-mode cavities each containing one quantum dot [13], photonic crystal microcavities [14], circuit quantum electrodynamics with a finite system approach [15] and in trapped ions [16]. One of the greatest advantages of all these setups is the combination of highly controllable experimental conditions and the large effective size of each site that allows for the design of mesoscopic simulators of condensed matter systems.
In many cases, the JCH system is naturally disturbed by noise that usually takes the system out of equilibrium. However, even in equilibrium, disordered imperfections in the system preparation may induce transitions that drastically change quantum phases and their correlations. Disorder may manifest itself in very different and even opposite effects. The lattice imperfections, which differ from site to site, may suppress quantum coherence, inducing the spatial localization of quantum states and destroying the system fluidity, which leads to compressible, though non-fluid, glassy phases [17]. However, disorder may induce fluidity under certain circumstances [6,18,19]. The hybrid nature of the system also leads to interesting effects under the action of disorder, as we show in the following sections.
In this paper, we show that the entanglement between different constituents of the JCH system can be used not only to characterize the already known phase transitions (also present in the Bose-Hubbard model), but also, and more importantly, to identify new behavior involving the nature, either hybrid or bosonic, induced by the more complex JCH interaction. We address small and large quantum systems, extremes that present similarities and differences that are of great interest: while a few sites are experimentally feasible in a controllable way, phase transitions are better defined in large samples. We also analyze the entanglement and disorderinduced effects of the JCH Hamiltonian. For the analysis of the small system, we go deeply into the statistics of the ensembles induced by disorder, since in principle a physical observer could carry out spectroscopic measurements of the system structure and obtain disordered pure states (or at least quasi-pure) pertaining to the induced ensembles. For the analysis of the large system, we resort to stochastic mean-field theory (SMFT), which was recently developed in [19,20] and allows us to study on-site statistics. We show how the statistics of the system changes under the various ways in which disorder may set in and also show the disorder-induced phase transitions.
The analysis of the clean system is developed in section 2 with one subsection for the small limit and another for the large limit. The disordered small system is addressed in section 3, and the disordered large system is addressed in section 5 after a brief introduction to SMFT in section 4. Section 6 concludes the paper.

The hybrid system: the Jaynes-Cummings-Hubbard Hamiltonian
The system studied here features a chain of sites, each of which contains composite excitations, also known as polaritons, created by the interaction of a boson and a fermion. A typical experimental proposal for these systems is devised in resonant cavities, the bosons being the photons that occupy the cavity mode and the fermions being two-level electronic transitions of the on-site dopants, as depicted in figure 1 and described in [2,21]. The Hamiltonian of these coupled doped cavities (with two-level impurities) is the so-called JCH model and it is given by with n being the number of sites and S i being the intra-site Jaynes-Cummings interaction between the dopant and the resonator: The ith site annihilation operators are a i and σ i for the bosonic and fermionic species, respectively. In equation (1), T (i, j) describes the photon hopping, or tunneling, between nearest neighboring sites The coupling strength between the two-level system and the cavity in the ith site is given by g i , and the photon tunneling strength between nearest cavities is A (i,i+1) . The photon frequency at the ith site is ω i and ν i is the transition frequency of the dopant of the respective site; thus we define the ith site detuning i = ν i − ω i . The polaritons are eigenstates of the intra-site (Jaynes-Cummings) Hamiltonian and are given by |n+ = sin(θ n )|g |n + cos(θ n )|e |n − 1 and |n− = cos(θ n )|g |n − sin(θ n )|e |n − 1 , with tan(2θ n ) = −g √ n/ . The states |n are photon number states and |g and |e are the ground and excited states of the dopant inside the cavity. Finally, the number of particles operator in the ith site is given by

Behavior of small sample systems
Recent works show that the system described in the last section undergoes a Mott-superfluid phase transition when going from small hopping to large hopping or from negative detuning to positive detuning [2,10]. In the first case, the transition is induced because the hopping strength 5 circumvents the photon blockage regime (nonlinearity due to the dopant-cavity interaction) and in the second case the excitations are directly driven from mainly electronic (electrons are not able to hop) to mainly photonic, hence the fluidity. This phase transition can be witnessed by single-site properties. For example, when the system is isolated and the average occupation number per site N i is one (the same number of excitations and sites), the variance of N i for any given site is a good order parameter [2,3]: in the Mott phase each site has a single particle and there is no number fluctuation, whereas in the superfluid phase the on-site number variance is maximum. This analysis begins to fail when one takes into account interactions with the environment and spatial fluctuations that may not preserve the total number of particles in the system. For instance, dissipation introduces variances of the occupation number in each site and var(N i ) may overestimate fluidity. Furthermore, although the measure of var(N i ) hints at the type of excitation that dominates each phase, it cannot reveal this fundamental property in detail because it does not fully distinguish between photons, electrons and polaritons. We proceed to show that the entanglement between different constituents not only reproduces previous results but actually allows for the identification of a new crossover that the previous analysis did not reveal.
In order to quantify the entanglement between the various components of the system, we choose the negativity measure [22], which is very convenient to calculate. It should be remembered that null negativity does not necessarily imply null entanglement; in fact, null negativity means null or bound entanglement. However, any non-zero value of negativity guarantees some form of distillable entanglement in the system that will prove to be enough for distinguishing the different quantum phases of the system as a whole. Given the quantum state of any two constituents A and B of the system, their negativity can be found by partially transposing their reduced density matrix R = ρ T A AB and then summing up the moduli of the negative eigenvalues of R.
We begin with the case of small systems, where the diagonalization of the Hamiltonian is computable by looking at properties of the lowest energy state |G with the constraint of having an equal number of excitations and sites H|G = E n |G . In other words, |G is the lowestenergy eigenstate of the Hamiltonian having n = n i N i , with N i being the number operator at the ith site. Consider now a cluster of two sites, which is the smallest possible such system. Even for this very basic unit cell, the entanglement between the sites clearly presents the signatures of Mott and superfluid phases that were found for much larger systems in previous works. In the Mott phase (with one polariton per site) there is no entanglement between sites with the Mott insulating state being |G(MI) = |1− |1− (for two sites). In the superfluid phase and when the excitations become mainly photonic, the sites become entangled, with the superfluid state for two sites (described in [3]) given by |G(SF) = |g |g [ 1 It should be kept in mind that for the finite system analysis there is no phase transition, only a smoother crossover, even though the phase transition terminology is commonly adopted. The behavior of the system (phase-like diagram), quantified by the entanglement between different constituents, is depicted in figure 2, where we show the entanglement between sites, the in-site entanglement and the entanglement between atoms.
The site-site entanglement shows the phase crossover as partially presented in [4]. When the site-site entanglement is negligible the system resembles a Mott insulator and when the site-site entanglement is non-negligible the system presents superfluid-like behavior. Thus the site-site entanglement indicates the regimes in which the system is insulating and superfluid with small and large values of entanglement, respectively (figure 2(a)). In order to quantify the polaritonic behavior we can look at the in-site entanglement that measures how correlated are the photonic field and the electronic transition in a given cavity (or a site) (figure 2(b)). In the Mottlike regime (small site-site entanglement) the in-site entanglement is significant, indicating that the system presents polaritonic behavior. Deep in the superfluid-like regime (large site-site entanglement) the in-site entanglement is small, indicating a predominant photonic behavior. However, during the crossover (as a function of either A or ) entanglement presents a nonmonotonic behavior, with a region where it is maximum. Such a non-monotonic increase, which is even more pronounced in the atom-atom entanglement (figure 2(c) and [23]), suggests that as the system size increases and reaches the thermodynamic limit a phase transition should be verifiable, i.e. since at the point of the phase transition there are fluctuations over all length scales, more degrees of freedom interact with each other such that entanglement can exist between more degrees of freedom. Furthermore, the regime in which in-site and site-site entanglement coexist corresponds to a polaritonic superfluid, rather then just a photonic superfluid.

Large sample systems and the polariton-photon crossover
We can also obtain a wider view of the system phase diagram varying the number of polaritons in the system. In order to do that we can couple the system to a chemical reservoir of polaritonic particles with chemical potential µ (at zero temperature), such that the system is in equilibrium with this reservoir. The chemical potential can be explicitly included in the system Hamiltonian For large (infinite) systems we adopt the mean-field approach, in which we treat a small cluster of sites interacting with a mean field, i.e., a classical approximation of the rest of the chain (which we refer to as an environment). This approach gives a factorable approximation of the non-factorable tunneling (or hopping) term by approximating the operators for their mean values plus a small fluctuation (in this case a quantum fluctuation) a = a + δa. The mean-field Hamiltonian for the cluster becomes with a = α being the mean-field order parameter that has to be self-consistently determined by minimizing the ground-state energy. The phase diagram as a function of the chemical potential and hopping frequency is shown in figure 3 for large systems.
Varying the chemical reservoir we can see the Mott lobes (each lobe corresponding to plateaus of different integer numbers of polaritons) in the infinite system (figure 3). The meanfield parameter is null in the Mott phase and is positive in the superfluid phase. Only in this case we compute a site purity (one minus the purity more precisely) as the estimate of the entanglement between such a site and the rest of the chain (in this case, the cluster). Although 8 this entanglement is not strictly zero in all of the Mott phase it still gives a fair account of the phase diagram and the lobe structure. We have considered a four-site cluster, which is a rather small cluster, even though it already requires considerable computational effort. Larger clusters would increase the precision of the site-cluster entanglement. The site-cluster entanglement is maximum in the lobe borders (the middle panel of figure 3), which, again, indicates strong polaritonic fluidity in the vicinity of the phase transition.
Now, looking at the in-site entanglement (the bottom of figure 3), which can be regarded as the very essence of the polaritons, we can see the whole picture with the overlay of the Mott-superfluid and hybrid-boson crossover. The highest in-site entanglement is in the Mott lobes and the lobe structure can also be defined by this quantity. Outside the lobes fluidity sets in; however, the in-site entanglement is still very high, indicating that the system has not yet undergone the hybrid-boson crossover despite having changed from insulating to superfluid. Farther away from the lobes and deeper into the fluid phase, we finally observe the in-site entanglement vanishing, indicating that the system finally turns bosonic.

Disordered small quantum systems
Every physical system presents imperfections (disorder); that is, the system parameters may vary from site to site. There are many possible origins of disorder, for instance, imprecisions in the system manufacturing process, thermal fluctuations and fluctuations induced by other uncontrollable electromagnetic sources in the system environment. One way to study the effect of disorder is to describe the parameters of each site as a stochastic variable ξ i and the Hamiltonian becomes dependent on the stochastic parameters H{ξ i }. Naturally, the system ground state becomes dependent on the values assumed by the system parameters |G → |G({ξ i }) and there emerges a new state, an average state that contains the statistics of the effects induced by the static disorder, with d p({ξ i }) being the distribution measure of the disorder, such that it gives all moments of the site parameters ξ k i = d p({ξ i })ξ k i . We choose to analyze only uncorrelated disorder such that the global measure is a product of local measures d p( The magnitude of disorder is then given by the distribution width or the mean square deviation δ. We can then characterize the average properties of the system given that it presents disorder. For instance, we can look at the entanglement description of the phase diagram; only now we average the entanglement over the pure state ensemble generated by the different values assumed by the system parameters. We remark that the pure state ensemble given in equation (5) is a physically realizable ensemble [24], since in principle the experimentalist can carry out spectroscopic measurements and obtain the values assumed by the parameters in that particular sample system and then prepare the system ground state. We can define the reduced states with the trace being performed over the environment of A and B. For instance, if we are looking at the atom-atom entanglement, then we trace out the field, so the field would be the environment in this case. Therefore, we can define the average

entanglement between any constituents A and B as
which is physically realizable since the ensemble of ground states is also physically realizable [25,26]. In what follows, in this section we consider only two sites of the JCH Hamiltonian.

Disorder in the matter-light detuning
Now we can describe the effects induced by disorder in each of the parameters individually, detuning , hopping A and matter-light coupling g. Let us begin by analyzing disorder only in the cavity-atom detuning; thus {ξ i } = { i } (see figures 4 and 5). As can be seen in figure 4, the average entanglement between sites seems to decrease over the entire phase diagram in comparison with the clean case of figure 2(a). The decrease of site-site entanglement indicates that the excitations tend to localize through an Anderson-like mechanism. For instance, starting from the system in the superfluid phase, when we increase the detuning disorder the distribution of the single-site number occupation P( N 1 ) is broadened and then becomes a two-peaked distribution (see the top panel of figure 5). In the regime in which the distribution P( N 1 ) presents two peaks the system is fully localized, such that one of the peaks corresponds to all excitations in cavity one and the other corresponds to zero excitations in the cavity. This extreme regime of localization can be regarded as a bosonic bunching: the large disorder in the cavity line width allows for realizations in which the cavity has a very low frequency such that it is energetically favorable to fit more than one excitation in one cavity instead of distributing the excitations over the sites. For two sites and two excitations the state can be approximately given by |2 |g |0 |g or |0 |g |2 |g with the atoms in their ground states. The parameters in figure 5 are such that the system is in the superfluid phase in the clean limit. In this case, as expected, the distribution of the site-site entanglement P(E[SS]) shows a peak at the maximum value (the bottom panel of the figure). However, as disorder increases a second peak emerges close to the minimum value corresponding to localized states. Thus the ensemble presents both superfluid and insulating states for intermediate values of disorder. The presence of the two kinds of state can be regarded as a precursor of a glassy phase [5], which we will show to be true with the large system analysis. Furthermore, as disorder increases even further the system becomes fully localized and the site-site entanglement distribution becomes single peaked at very small values of entanglement.

Disorder in the photon hopping
Disorder in the photon hopping generates a different effect (see figures 6 and 7). Carrying out the same analysis as for disorder in the detuning, we see that the average site-site entanglement increases in the region where a Mott phase exists in the clean limit, while it remains practically unaltered in the superfluid phase. Fluctuations in the hopping may actually induce a glassy fluid phase [6,18,19]; that is, disorder allows realizations in which the hopping is stronger than the photon blockade and those realizations may prevail. The A-disorder may also suppress the polaritonic behavior which can be seen as a decrease in the atomic population. We can also look at the distribution for the total atomic excitation Z = i σ † i σ i as a function of disorder starting from the system at the Mott phase in the clean limit. The distribution P( Z ) is very  asymmetric and as disorder increases it concentrates at the extreme values assumed by Z . One of the extremes corresponds to polaritonic superfluidity and the other to photonic superfluidity, with the latter prevailing in the limit of very large disorder. This can be corroborated by the entanglement distribution P (E[SS]).

Disorder in the matter-light coupling
Finally, we analyze disorder in the matter light coupling g (see figure 8). A first look at the g-disordered average entanglement (not shown here since it resembles very closely the -disorder case) suggests that disorder in the Jaynes-Cummings coupling also induces localization; that is, disorder allows a great number of meaningful realizations in which the sites are almost unentangled and the excitations tend to bunch. However, the g-disorder-induced distribution of the site occupation number P( N 1 ) (top panel of figure 8) is very different from that induced by disorder. In the current case the P( N 1 ) distribution shows three peaks, the extreme ones corresponding to bunching similar to the -disorder case and a middle one that corresponds to states in which the excitations are still equally distributed among the sites. Disorder in the matter-light coupling induces states with Mott-like features in which the site-site entanglement vanishes (see the middle panel of the figure). In fact, this is the meaning of the middle peak in P( N 1 ): some sites undergo a superfluid-insulating transition through a Mott-like mechanism. Since the distribution of the atomic occupation P( Z ) (bottom panel of figure 8) is narrowly centered at an appreciable (although not extreme) value we may conclude that the system nature becomes mainly polaritonic, and thus presenting both Mott states (middle 13 peak in P( N 1 ) and polaritonic bunching (the extreme peaks in P( N 1 ). This suggests that the system behaves very similarly to an Anderson-Mott insulator [27].
The presence of superfluid, Mott and Anderson-like states for intermediate disorder in the coupling g suggests that glassy phases would be induced in larger systems of the JCH type. In fact, a different situation was analyzed in [28], in which there is disorder in the number of impurities per cavity. Since the number of atoms fluctuates, the intensity with which light couples to matter also fluctuates, and it was shown in [28] that such disorder induces glassy phases. Interestingly, even small versions of the quantum system present evidence for many diverse phases and behaviors expected only for larger quantum systems, which we discuss in the next section.

Disordered large quantum systems
To address the physics of large disordered quantum systems, we apply a recently developed technique, namely SMFT [20]. This method has been shown to provide appropriate descriptions of the effects of disorder without overestimating coherence and fluidity and has already been successfully applied to the disordered Bose-Hubbard model [20]. The main reason for the effectiveness of the method is self-consistently determining the probability distribution for the mean-field parameter P(α) (instead of α itself) through an iterative process.

Stochastic mean-field theory
Firstly, we describe how to account for any on-site disorder (with constant photon hopping A); afterwards we describe the special case of hopping disorder. In the mean-field description every site has a number z of nearest neighbors. The mean-field Hamiltonian for the kth site depends only on the scaled sum The probability distribution for η (we drop the site index for convenience) can be found from a simple and fundamental relation known as the convolution theorem which reduces to the Fourier transforms ϕ(β) = dα P(α) e iβα and The first step in the algorithm is to choose a trial distribution for α (different from a delta centered at α = 0) and the desired distribution (in our case a Gaussian) for the disordered parameter ξ (the detuning or matter-light coupling). Then we assume all α j to be independent of each other and we determine the self-consistent distribution such that a = G[ξ, η]|a|G[ξ, η] , with dq(η) = dηQ(η). The procedure is iterated until we observe convergence, that is, until P (i) (α) in the ith step is statistically close to P (i+1) (α). Finally, the average state of the site is given by the disorder-induced ensemble To account for disorder in the photon hopping we must add another step to the procedure. It is convenient to work with the variable φ = Aα whose probability distribution is given by (we use a subindex to distinguish the various distributions) Then we can determine Q(η) through the usual Fourier transforms ϕ(β) = dφ P φ (φ) e iβφ , (13) and the procedure follows as for the case of on-site disorder.

Disorder-induced transitions
The results presented in this section will allow us to conclude that the asymptotic effects of disorder (very large disorder) in the thermodynamical limit are very similar to the effects in small samples of the system. However, there are some significant quantitative differences; for instance, there are in fact phase transitions induced by disorder in the thermodynamical limit. It should also be pointed out that our approach is slightly different in this section. From now on, we work with single-site mean-field theory rather than cluster-mean-field theory, and we follow this strategy to avoid higher computational demands. This limits the applicability of the method and quantities such as the site-cluster entanglement are no longer addressable. Nonetheless, we are able to increase the local effective dimension of the oscillator to 20 states. Another difference between the approaches for small and large samples is that in the first case we fix the number of excitations in the system and it does not change as disorder increases. This is not the case in the present section, and in fact, the total number of excitations may change as a function of disorder.
Using the SMFT approach we are able to perform an analysis of the thermodynamical limit. The method allows us to recover the probability distributions (under the single-site meanfield approximation) for the various quantities we analyze for characterizing the system. The average mean-field parameter, for instance, can be readily evaluated as α = α P(α) dα. We follow the same ordering of the presentation of the results: First, we show the results for the detuning disorder, then for the hopping disorder and finally for the matter-light coupling. Given the unlimited nature of the disorder distribution we analyze (Gaussian), it follows that the insulating phases we present below are of a glassy nature. Such phases have non-vanishing number variance as opposed to the Mott-insulating phases [19]. Therefore glassy insulators can be characterized by vanishing superfluidity and non-vanishing compressibility (which can be related to the number variance). However, we do not show the compressibility of the system, since the result can be readily anticipated. The compressibility increases in the insulating regions; thus glassy phases are established.

Disorder in the matter-light detuning
As shown in figure 9, the net effect of disorder in the detuning is to induce insulating behavior, indicated by the destruction of the fluid phase surrounding the Mott lobes; in fact, the lobe structure disappears for significant amounts of disorder. The in-site entanglement shows that the system remains in superpositions of light and matter excitations for intermediate values of disorder. However, as we increase disorder the in-site components are either highly entangled or unentangled with higher probability. We can see the distributions as functions of disorder in figure 10. The transition from fluid to insulating is evident in the distribution of the meanfield parameter. All this corroborates the small system predictions and once again the system is mainly photonic for strong disorder. Interestingly, in the present limit the distribution of cavity excitation (P N in figure 10) is a series of delta functions (with different weights) centered at integer values of the mean occupation, which is in agreement with the insulating and bunched behavior.

Disorder in the photon hopping
The effects that hopping disorder produces in the system are the opposite of those produced by detuning disorder, as is the case for small systems. In the current case, disorder induces a fluid phase, as we can see in figure 11, and decreases the in-site entanglement, indicating that the system becomes more photonic in nature. As in other situations analyzed, by looking at the distributions of the physical quantities as a function of disorder (not shown), we were able to observe the transition from insulating to fluid behavior as disorder in the hopping parameter increased.

Disorder in the matter-light coupling
Finally, the effects of disorder in the light-matter coupling once again resemble those induced by the detuning disorder. Even though both the detuning and coupling disorders  induce insulating behavior, it should be pointed out that they do it through very different physical mechanisms.
In the detuning case cavities may be at lower frequencies in many sites, which allows for more photons to localize (even several photons per cavity). In the coupling case, the strength of the polaritonic nature may be increased in cavities that are strongly coupled to their corresponding matter components inducing Mott behavior. And as we can see in Adding the information provided by the entanglement and atomic population distributions, we have found that a fraction of the sites assume the Mott behavior and the rest are localized or even empty. Thus, we corroborate the small system analysis that suggests that the system behaves very similarly to the Anderson-Mott insulator [27]. It is, however, strikingly interesting that in one case (disorder in the detuning) it is the low-frequency sites that dominate the resulting statistical behavior and in the other case (disorder in the matter-light coupling) it is the strongly coupled sites that dominate, even though the distribution of the disorder parameter is Gaussian and unbiased.

Conclusions
We were able to characterize the phase diagram and the Mott-superfluid transition of small and large samples of the JCH Hamiltonian using entanglement measures between the various possible partitions of the components of the system. In particular, we showed that these nonlocal measures identify more clearly where the transition happens. Furthermore, and more importantly, we also showed that entanglement measures distinguish which type of excitation dominates each phase, which in turn allowed us to identify a crossover that is particular to this 18 hybrid system and does not have a purely bosonic analogue. This behavior splits the superfluid regime into two, the first one dominated by polaritons and the second one purely photonic. For the disordered system, we have shown that the simple statistical treatment of small systems can be quite instructive and allows us to draw conclusions that can be corroborated by the large system limit. We have shown that disorder in both the light-matter detuning and light-matter coupling induce insulating phases; however, they do it through very different physical mechanisms. The former allows for photonic localization and bunching and the latter induces Mott behavior in a fraction of the sites that prevents the bunching. Furthermore, the cavity-cavity coupling disorder induces a glassy fluid phase. The rich in-site structure of the system leads to these diverse disordered phases with very different statistics and physical meanings.
A great deal of work remains to be done on the characterization of the JCH system. As a valuable point we suggest that an appropriate and efficient method should be applied to the study of the large hopping limit (with and without disorder) in which the mean-field approach adopted here is limited by the truncation of the state space.
Finally, it is worth mentioning once again the mesoscopic aspect of the systems proposed to implement the JCH Hamiltonian as well as the increasing ability to manipulate the different parameters of these systems, sometimes even at an individual level. These properties suggest that it will be possible to carry out a thorough experimental investigation of the effects of disorder and its relation to phase transitions and entanglement in many-body physics in the near future.