Accelerated optimization problem search using Bose–Einstein condensation

We investigate a computational device that harnesses the effects of Bose–Einstein condensation to accelerate the speed of finding the solution of optimization problems. Many computationally difficult problems, including NP-complete problems, can be formulated as a ground state search problem. In a Bose–Einstein condensate, below the critical temperature, bosonic particles have a natural tendency to accumulate in the ground state. Furthermore, the speed of attaining this configuration is enhanced as a result of final state stimulation. We propose a physical device that incorporates these basic properties of bosons into the optimization problem, such that an optimized solution is found by a simple cooling of the physical temperature of the device. Using a semiclassical model to calculate the equilibration time for reaching the ground state, we found that this can be sped up by a factor of N, where N is the boson number per site. This allows for the annealing times for reaching a particular error to be systematically decreased by increasing the boson number per site.


Introduction
Quantum computation promises to offer great increases in speed over current computers due to the principle of superposition, where information can be processed in a massively parallel way [1]. One of the problems that is currently encountered is that such superposition states, which include highly entangled states, decohere very quickly due to the unavoidable coupling to the environment. On the other hand, quantum indistinguishability [2] of particles, another fundamental principle of quantum mechanics, remains relatively unexplored in the context of information processing. Bosonic indistinguishability is the mechanism responsible for fascinating physical phenomena such as Bose-Einstein condensation [3]. Can such quantum indistinguishability be used as a resource for quantum information processing? In this paper, we explore this possibility and show that it may be used to accelerate the time required to find the solution of an optimization problem.
The basic concept of our work is as follows. As is well known from quantum optics and condensed matter physics, identical bosonic particles have the property that they 'like' to occupy the same quantum state. For example, in Bose-Einstein condensation a macroscopic number of bosons all occupy the ground state of the system. This property of bosons suggests that they may be useful in the context of optimization problems in the following sense. Typically, optimization problems are formulated in terms of a cost function to be minimized. Now if the cost function happens to be the Hamiltonian of a particular system, and the system is composed of bosons, then we may expect that the bosonic effect will cause the system to preferentially occupy the ground state. Since the ground state is the minimum energy state of the Hamiltonian, which is in this case the cost function, the state of the bosons gives the solution of the optimization problem.
Another bosonic effect that is also beneficial in solving optimization problems is that bosons have an amplification factor for transitions between states. For example, the transition amplitude for a boson entering a state already occupied by N bosons is amplified by a factor of √ N + 1. In the context of lasers and Bose-Einstein condensates (BECs), this effect is known as final state stimulation, where the rate of cooling into a state occupied by N bosons is increased 3 by a factor of N + 1, according to Fermi's golden rule [5]. This means that the cooling time required for entering the ground state is reduced when a large number of bosons are present in the system. This suggests that the time for cooling into the ground state of a Hamiltonian encoding an optimization problem may be reduced when using a large number of bosons.
The above concepts are formalized to a specific proposal in the following sections. This paper is organized as follows. In section 2, we discuss the basic proposal of the system and the bosonic Ising model Hamiltonian, the fundamental Hamiltonian of the system. The protocol for running the proposed device is discussed. In section 3, we discuss the properties of the bosonic Ising model at thermal equilibrium and show the first bosonic property discussed above, the enhanced occupation of the ground state. In section 4, we discuss the reduced cooling time of the proposed device due to the final state stimulation effect discussed above. The performance of the device using a thermal annealing procedure is analyzed in section 5. In section 6, we discuss possible methods for implementing the bosonic Ising Hamiltonian using a measurement feedback approach. We summarize our results in section 7.

The proposed system
We first specify the optimization problem that we would like our proposed system to solve. The Hamiltonian, or the cost function to be minimized, is where J i j is a real symmetric matrix that specifies the connections between the sites i, j, and σ i = ±1 is a spin variable. The task is then to find the minimal energy spin configuration {σ i } of the Hamiltonian (1). Even though the Hamiltonian (1) is classical in the sense that all terms in the sum commute with each other, finding the ground state is computationally difficult due to the exponential number (2 M ) of spin configurations. The difficulty of the problem lies in the fact that typically there is little structure in the energy landscape that can be exploited to find the global minimum. The problem of finding the solution of the Hamiltonian (1) is, in fact, known to be NP-complete, since it can trivially be mapped onto the MAX-CUT problem [4]. Furthermore, it can in principle encode an arbitrary NP-complete problem by a polynomial time mapping procedure [6]; thus the potential application of a reliable method of solving (1) is extremely broad. For these reasons, it has previously been considered for applications in simulated annealing [7] and quantum annealing [8].
In order to solve (1) we consider a system such as that shown in figure 1. Each spin σ i in H P is associated with a trapping site containing N bosonic particles. The bosons can occupy one of two spin states, which we label as σ = ±1. Due to bosonic indistinguishability there are N + 1 possible states given by where a iσ is the annihilation operator for spin σ on site i. On each site i we may define a total spin operator S i = a † i+ a i+ − a † i− a i− taking eigenvalues S i |k = (−N + 2k)|k . We consider a semiclassical picture where no superposition between the states is assumed, but transitions between the states can occur due to thermal relaxation. One of the main effects that we would like to take advantage of is bosonic indistinguishability, so BECs are a natural system for an J ij Figure 1. The proposed system. Each site of the Ising Hamiltonian is encoded as a trapping site containing N bosons. The bosons can occupy one of two states σ = ±1, depicted as either red or blue. The sites interact with each other via Hamiltonian (3). The system is assumed to be coupled to a heat bath such that excess energy can be dissipated.
experimental implementation of such a system. Recent experimental advances have allowed for the realization of two-component BECs [9][10][11], which are a possible candidate for a single trap in figure 1. Another system that is promising is exciton-polaritons in semiconductor microcavities, which were recently observed to undergo Bose-Einstein condensation [12][13][14]. Exciton-polaritons possess a spin of σ = ±1, which serves as the two quantum states of the BEC. Systems that undergo BEC are natural choices for the implementation of such a device, since such systems possess large bosonic amplification factors, the primary resource that we wish to take advantage of. Now assume that the sites are controlled such that the system follows the bosonic Ising Hamiltonian where J i j is the same matrix as in H P , which specifies the computational problem. A typical run using such a device proceeds as follows (see table 1). Initially, each site is prepared with a random population of σ = ±1 spins. The system is then cooled in the presence of Hamiltonian (3) for a time τ a according to the annealing schedule T (t), where T is the temperature. The fact that we assume cooling to take place in the system means that the bosonic system is explicitly coupled to a heat bath which absorbs excess energy from the system. After the annealing procedure is finished, the spin of each site is read out by taking a majority vote on each site. The process is repeated many times in order to take statistics on the final result. The procedure is essentially identical to a simulated annealing procedure, the sole difference being the use of the bosonic Ising model. We shall see that the use of a large number of bosons on each site incorporates the two effects mentioned in the introduction, namely the enhanced ground state probability occupation and accelerated cooling. Table 1. Procedure for the accelerated optimization scheme using the system shown in figure 1.
Step Procedure 1 Prepare the system in a random initial configuration 2 Apply Hamiltonian (3) with an annealing schedule T (t) for a duration t = [0, τ a ] 3 Read off spins of each site S i and record values 4 Repeat steps 1-3 a large number of times until the desired statistical accuracy is achieved In this paper, we shall consider the class of problem Hamiltonians The energy landscape of such a Hamiltonian can be understood by noting that this may be rewritten in terms of the total spin variable Assuming J λ > 0, the first term of (5) is minimized by a configuration of all up or all down spins, while the second is minimized by all up spins. In terms of the computational problem to be solved, the solution is therefore the configuration with all up spins. Therefore the energy landscape consists of two nearly degenerate ground states separated in energy by E = 2Mλ. The highest energy configuration is the state with σ tot = 0. To move between the two lowest energy states, an energy barrier of ∼J needs to be overcome.
One potential pitfall of using (3) in place of (1) is that the nature of the original problem to be solved is fundamentally changed. For example, what guarantees that the solution of (3) has any resemblance to the solution of (1)? In fact, it can be shown that the ground state spin configuration of (3) is equivalent to the original Ising model Hamiltonian H P . We defer the proof to appendix A. The argument relies on the fact that there is a linear transformation of the energy (by a factor of N 2 ) between Hamiltonians (1) and (3). Thus there is no reordering of the energy levels between the two Hamiltonians. This property allows us to use the modified Hamiltonian (3) while still extracting the results for the original problem Hamiltonian (1).

The Bosonic Ising model at thermal equilibrium
We now turn to examining the physical properties of the bosonic Ising Hamiltonian (3). Consider the scenario where the Hamiltonian is turned on for a very long time such that the system is at thermal equilibrium with the heat bath at a temperature T . According to the discussion in the introduction, we expect that the bosons have a higher ground state occupation probability than classical distinguishable spins.
For simplicity, let us first consider the minimal two-site case. Performing the bosonic mapping procedure of the previous section, we obtain the Hamiltonian where the summation is taken corresponding to indistinguishable bosonic particles or classical distinguishable particles (see appendix B). Physically, the comparison illustrates the difference between using two-component BECs, as discussed in section 2, and an ensemble of classical spins that do not possess quantum indistinguishability. We see that the bosonic case has a larger average spin for N > 1 and all temperatures, corresponding to a spin configuration closer to the ground state. As the particle number is increased, the temperature required to reach a particular S i increases. For the bosonic case, the required temperature increases linearly with N , while for distinguishable particles it behaves as a constant for large N .
This behavior is in fact true for any number of sites, not only the two-site Hamiltonian illustrated here (see appendix B). The concentration of particles in the ground state configuration for bosons is precisely the same effect that is responsible for the formation of a BEC. Since the ground state corresponds to the solution of the computational problem, this corresponds to an enhanced probability of obtaining the correct answer at thermal equilibrium. From a practical perspective this results in an improved signal-to-noise ratio for the bosons in comparison to distinguishable particles at the same temperature.
In our calculations, we have neglected the presence of an on-site particle interaction ∝ S 2 i that has the qualitative effect of degrading the spin signal-to-noise ratio. For atomic systems the effect of such unwanted interactions may be minimized by the application of a suitable Feshbach resonance to control the interactions in the case of an implementation using atoms. For example, by tuning the interactions between the two bosonic species to be the same as the interactions within a species, such a term may may be entirely eliminated [20]. For exciton-polaritons, the natural interaction is rather weak such that using sufficiently large traps and a strong feedback signal, the term can be made negligible.

Equilibration time
We now turn to the time taken to reach thermal equilibrium, after initially preparing the system with random populations of σ = ±1 particles on each site. The time evolution is performed by an extension of the semiclassical methods presented by Glauber [15] to the bosonic Ising model (3), also accounting for transitions beyond first order in perturbation theory. Given the where k i ranges from 0 to N , and we have defined the vector k = (k 1 , k 2 , . . . , k M ). The state can thus be imagined as a point within an M-dimensional hypercube. At t = 0, there is an equal probability p k of all states |k . The probability distribution then evolves according to where δk i = (0, . . . , 0, δk i , 0, . . . , 0) is a vector in the direction of the ith axis of the M-dimensional hypercube. The w(k, δk i ) is a weight factor for the process |k → |k + δk i , containing a transition rate factor from Fermi's golden rule and a coefficient to ensure that the system evolves to the correct thermal equilibrium distribution, according to detailed balance [15]. Explicit expressions for the weight factors are given in appendix C. The time evolution equations are evolved numerically for small boson numbers and a kinetic Monte Carlo routine for larger boson numbers. To see the accelerated cooling effect in action, let us first examine first a minimal single-site system. Figure 3(a) shows the cooling of the system for N = 1 and N = 50 particles. Firstly, we see that as the number of particles is increased, the time taken to reach equilibrium is reduced. This reduced cooling time effect is precisely the final state stimulation effect, originating from the amplified transition amplitudes in Fermi's golden rule. Secondly, we see that at equilibrium the average spin for the N = 50 case is very close to the ground state value. This effect is the same as that discussed in the previous section, where a large number of bosons enhances the ground state occupation probability. In the long time limit t → ∞ the system converges to precisely the thermal equilibrium distribution discussed in the previous section.
The single-site case is in fact simple enough to study by analytic means. For low temperatures, the time dependence of the single-site case can be approximated by the rate equations where n i are the populations on levels i = 1, 2 and α is the thermalization rate. Analytically solving for large N under the condition n 1 + n 2 = N gives giving an equilibration time of τ ∼ 2/α N . We thus see that there is an accelerated cooling effect that is ∝ N by using a large number of bosons.
As is well known from past studies of simulated annealing, the presence of local minima slows down the time for equilibration dramatically. To illustrate the behavior of the system in the presence of local minima, we examine the four-site Ising model problem (4) that has an energy barrier separating the two nearly degenerate ground states ↑↑↑↑ and ↓↓↓↓. We define the error probability as the probability of failing to obtain the correct ground state configuration ↑↑↑↑ after a single measurement of the total spin where Z is the partition function. The summation is over all configurations with the same sign of spin as the ground state {σ i } of the original Ising model problem (1). In general, a higher temperature gives a larger error rate due to the higher probability of occupying excited states with a different spin configuration to the ground state. The effect of the local minimum on the equilibration time can be seen from the N = 1 case shown in figure 3(b), where the time rapidly increases as → 0 (i.e. T → 0). This is in contrast to the single-site case where there is no such dramatic increase in the equilibration time, and is a manifestation of the local minimum in the problem. A higher temperature reduces the equilibration time due to thermal fluctuations by more efficiently transferring populations between the local and global minima. Figure 3(b) shows that as the boson number is increased, there is a speedup at constant error of several orders of magnitude. There is a small odd/even effect due to the definition of the error , which does not include the S i = 0 configuration in the summation. The curves approach zero equilibration time ∝ 1/N for large N ( figure 3(c)). In 9 all our numerical simulations, we have found that bosons are able to speed up the equilibration times, with a ∝ 1/N speedup for large N , in a similar way to the single-site example.
The origin of the speedup can be understood in the following simple way. The use of many bosons increases the energy scale of the Hamiltonian from ∼ J i j to ∼ N J i j . Due to bosonic statistics, the coupling of the spins to the environment is increased by a factor of ∼ N . Thus by constructing a system out of bosons, we have increased the energy scale of the entire problem by a factor of N , which results in a speedup of N . Spin flips due to random thermal fluctuations also occur on a timescale that is faster by a factor of N , resulting in a faster escape time out of local minima. We have confirmed this numerically in our simulations by measuring the average spin flip time τ flip , defined as the average time for the total spin S i on a particular site to change sign. The results in figure 3(d) show clearly a ∝ 1/N behavior, signifying that with an increased number of bosons, there are many more random spin flips per unit time. We can thus picture the system moving through the space of spin configurations with a speedup proportional to N . Since in a typical BEC there are a large number of bosons, the speedup can be quite large in practice. For example, in atomic systems typically N ∼ 10 5 or more bosons can be put in a single trap.

Annealing
In the previous section, we calculated the equilibration time at a fixed temperature. As discussed in section 2, in a practical setting an annealing procedure would be used to gradually decrease the temperature of the system. We note that in practice it may be more convenient to gradually increase the overall energy scale of the Hamiltonian (by multiplying by a factor 1/T (t)) at constant temperature, rather than varying the temperature. Both are equivalent procedures that give identical results. We expect that the enhanced cooling of the system, incorporating the bosons as found in the previous section, should improve the performance of the system, including annealing.
We employ an exponential annealing schedule according to where T 0 (N ) is the temperature corresponding to an error probability of 0 for boson number N , and τ a is the annealing time. Different initial temperatures are necessary so that systems with differing N have equivalent starting points in terms of the error probability. After the system is evolved for an annealing time of τ a , and the error probability is calculated according to equation (12). The exponential annealing schedule (13) is typically chosen for studies in simulated annealing [7], and therefore we expect this to work reasonably well in our bosonic case. Figure 4(a) shows the average spin on a single site during the annealing process for various values of τ a . We see that for a particular annealing time the average spin monotonically approaches the ground state value of S /N = 1. Towards the end of the annealing procedure, S plateaus and the system no longer responds to further cooling. This can be interpreted to be a freezing of the system due to the equilibration time exceeding the time available in the annealing procedure. This is expected in any annealing procedure when the annealing schedule (13) is in a low-temperature range toward the end of the annealing time. As expected, for longer annealing times the average spin increasingly approaches the ground state value of S /N = 1.
The error probabilities post-annealing for various N are shown in figure 4(b). The error probabilities monotonically decrease with increasing boson number N and increasing annealing times. The improvement with N can be directly attributed to the enhanced cooling of the system, as discussed in the previous section. Early in the annealing procedure, the system is typically at a high enough temperature such that the spin configuration can continuously adapt to the cooling imposed by (13). However, at some point the annealing rate becomes too fast for the system to respond, a timescale determined by the equilibration time τ . The enhanced cooling due to the many bosons in the system allows for a lower error probability to be reached, as seen in figure 4(b). Figure 4(c) shows the annealing time required to achieve a particular error probability post-annealing. We see that several orders of magnitude improvement in the annealing times are possible even for the relatively small numbers of bosons that we use in the simulation. This confirms our initial conjecture that the decreased equilibration times lead to reduced annealing times to reach a particular error probability.

Implementation of the bosonic Ising model Hamiltonian
Up to this point, we have not discussed how such a Hamiltonian (3) could be implemented in practice. In general, the form of the interactions in (3) is long ranged and highly connected between all sites. This makes implementation in a real device using natural interactions between particles difficult. Here we propose one possible solution for implementing the interaction Hamiltonian (3) via a system of measurement and feedback control. The total spin S i on each site is continuously measured, and an appropriate control signal is continuously fed back into the system by applying a local dc field on another site. The basic concept of the scheme is as follows. Say at a particular instant a spin measurement of all the sites is made, giving the result {S j (t)}. Then, at that moment a local field B i (t) = j J i j S j (t) is applied on site i. Since the Zeeman energy due to this field is H = i B i S i , a simple substitution yields (3). Although J i j has a large connectivity and is long ranged, by using such a feedback method to induce the interactions there is no restriction on the kind of interactions J i j that can be produced in principle.
The above argument can be formalized using the framework of quantum feedback control. We start with the Wiseman-Milburn feedback master equation [16] where L 0 is a Liouville superoperator describing the internal dynamics of the system, D[C]ρ = C C † − {C † C, }/2 is the Lindblad superoperator, C is the measurement operator due to the meter coupling, is the detector efficiency, M is the measurement superoperator and c is the density matrix of the system conditional on prior measurement outcomes. We consider Markovian feedback and the system is acted on by a Hamiltonian H fb (t) = I (t)F, where I (t) is the feedback current due to the measurement outcome [17].
We now define each of the variables in the master equation for our specific implementation. Our system consists of a set of cross-coupled systems such as that shown in figure 5. First consider one particular site i. Without the feedback control, we assume that there is no internal dynamics of the system L 0 = 0. The meter measures the z-component of the spin; thus we where γ is the rate constant representing the measurement strength, n i− is the number operator counting the number of down spins on site i, and we have assumed N bosons per site. The second term of the master equation thus describes a dephasing term originating from the measurement of the z-component of the spin. The back-action of the z-measurement gives a measurement superoperator Mρ = S z i ρ + ρ S z i . As a result of the feedback, on each site we apply a field in the z-direction such that F ∝ S z i . Now consider the complete feedback system as a whole. Consider applying a feedback Hamiltonian of the form where I j (t) is the current resulting from the measurement of site j and is an overall constant. Inserting these expressions into the feedback master equation gives The first term is a coherent evolution of the system according to the Hamiltonian (3) [18]. Each site is coupled to an external heat bath which dissipates energy according to a coupling ± dependent on the total current j =i J i j I j entering site i and temperature T . These terms ensure that the system approaches the correct thermal equilibrium distribution in the long time limit. The diagonal components of the master equation then evolve according to (9). The fourth term is a dephasing term originating from the measurement on each site, as well as a contribution from the feedback circuit noise.
In the semiclassical context, (16) produces the correct form of dissipation, which is sufficient for the operation of our device. Therefore, the introduction of additional dephasing terms into the system should not adversely affect the system since we do not assume any coherence across the sites from the outset. However, in a more general scenario where coherence exists between the sites, such dephasing terms will act to reduce coherence in the device, which may affect the performance. Another potential problem is the existence of a finite delay time in transferring information between sites. This can potentially cause oscillatory or chaotic behavior depending on the parameter values chosen, as is known from classical control theory [26]. In the cases where there is coherence across the system, another question is whether such a quantum feedback approach would be able to produce the correct coherent dynamics. We leave a full investigation of the effect of such a measurement feedback approach for various situations to future work.

Summary and conclusions
We have investigated a device that uses an array of two-component BECs to solve optimization problems. The key quantum effect that is used in the system is indistinguishability, rather than quantum superposition or entanglement. The proposed system has several advantages over an equivalent system that does not include indistinguishability. At thermal equilibrium, the system has a greater probability of occupying the ground state, due to bosonic statistics preferentially occupying the low-energy states. The equilibration time is found to be reduced by a factor of 1/N , due to the effect of final state stimulation between the bosons. A key reason for the speedup is the enhanced spin flip rate proportional to N , which allows for the system to escape 13 local minima at this enhanced rate. The improved cooling time allows for faster annealing times to be performed at the same error probability. The speedup tends to be more pronounced for small error probabilities, with several orders of magnitude' improvement observed in our numerical simulations. We have also proposed a possible method for implementing the bosonic Ising model via a measurement feedback approach. This is rather crucial for a physical implementation of the device, since the Ising interaction matrix J i j contains highly non-local interaction terms that are difficult to realize with only natural interactions. We note that the analysis presented in this paper was restricted entirely to the semiclassical case where there was no coherence between the sites, which can potentially positively or negatively impact the performance of the device. We leave extensions involving coherence to future work.
Although the device discussed in this paper is a computational device that uses quantum effects, it is not a quantum computer in the conventional sense, since the off-diagonal density matrix elements of the state of the device are explicitly zero at all times. The scheme is more closely related to a physical implementation of a simulated annealing algorithm, rather than quantum annealing or quantum adiabatic computation. Considering the wide applications of simulated annealing to optimization problems, we expect that the proposed device can also be applied in similar contexts. The resources required for realizing the feedback circuitry only grow as ∝ M 2 , which follows from the number of connections between BEC sites. The speedup of the device can be increased by simply increasing the number of bosons N in a site, with no additional resources in terms of apparatus. This should be compared to a completely classical approach to the problem, where either the computational circuitry must be made smaller or massive parallelization is required in order to reduce the simulation time.
One key advantage that the device has over the standard quantum computer is that the 'computation' is performed in an explicitly open dissipative system. In spirit it is similar to the method of quantum simulation, where a physical implementation yields information about a given Hamiltonian [21]. In contrast to a quantum computer where decoherence is treated as unwanted noise, in this device the dissipation drives the computation towards the desired state [22,23]. Since there is no coherence between the sites the experimental demands are far less stringent than for making a full quantum computer. Considering devices that are able to harness quantum effects in such open dissipative systems may lead to novel types of information processor that outperform classical computation devices.
14 Writing S k = |S k |σ k , where σ k = sgn(S k ), and using the fact that if S is the ground state, we have Now consider another configuration where the magnitude of the kth spin is changed, but its sign is not changed: S = {S 1 , S 2 , . . . , |S k |σ k , . . . , S M }. The energy difference is where we have again used the fact that S is the ground state. , there is a degeneracy in the ground state and the spin on site k can take an arbitrary value. A simple example is a three-site system with J 12 = J 23 = J 13 > 0. In this case, one of the spins can take an arbitrary value, with the other two taking maximally antiparallel spin configurations. Now we would like to connect the above result with the ground state of the original Ising Hamiltonian (1). Let the energy spectrum of (1) be E n (n = 0, 2 M − 1), with E 0 E 1 · · · E 2 M −1 . For maximized spins, the Hamiltonian (3) has exactly the same energy spectrum, except that it is multiplied by a factor N 2 . For sites corresponding to the case (A), from the above result, the spins are maximized S k = N σ k . For case (B) sites, due to the energy degeneracy, we may assume that one particular ground state configuration can be written with maximized spins. It then follows that this state must have an energy N 2 E 0 and has the same σ k configuration as (1).

Appendix B. Average single-site spin of an M -site Hamiltonian
We consider the general case of a Hamiltonian of the form It is convenient to rescale the total spin by a factor of two S i = S i /2 so that the spin runs between S i = [−N /2, N /2] in steps of one. The partition function for both bosonic and distinguishable particles is Looking at the distinguishable case, we may absorb the combinatorial factor using Stirling's formula Dropping the constant term and absorbing the first term into the exponential of equation (B.2), we obtain an extra term in the Hamiltonian This gives us a direct understanding of why the distinguishable particle case has a lower proportion of particles in the ground state. The i S 2 i term favors spin configurations that are near S i = 0 rather than the extrema S i = ±N /2. At low temperatures β → ∞, the effect of this term is diminished, and one recovers the same result as the bosonic case.
Turning the sums into integrals (which is valid for N 1), the partition function is The quadratic form of the spin variables in the exponent suggests that we change variables by diagonalizing the matrix J = U U T , where is a diagonal matrix with elements ξ i i = 1, M.
Defining new spin variablesS = U S, the partition function is where R S is the region of integration under the transformed coordinates. To approximate this integration, we take a hypercube of the same volume as the original integration, along theS i axes. Using these variables, the integration variables decouple, giving where and we have introduced variables The average spin on a given site can be derived using the same approximations as that used for the partition function, which yields where and σ k is the ground state spin configuration for the site k. Completing the square in the exponential, we calculate . (B.14) For large N the second term is small; thus the temperature dependence (if any) will be dominated by the first term. Substituting the first term into equation (B.12), we obtain Thus for distinguishable particles ( = 1), the temperature only appears as the product β N . This implies that as N is increased, the average spin S k can be kept constant if k B T /N ∝ 1, as seen in figure 2(b). For bosonic particles ( = 0), the temperature dependence cancels out, so that we need to examine the second term in equation (B.14). Examining the temperature and N dependence of the terms within the exponents, it is easily verified that all terms have a ∝ N 2 β dependence, while the terms in the error function have an N √ β dependence. Similarly, the factor of 1/ √ A i , combined with the 1/N factor in equation (B.12), has a dependence 1/N √ A i ∝ N √ β. Thus, for the bosonic case, the temperature has a dependence k B T /N ∝ N , as seen in figure 2(a).

Appendix C. Time evolution of the bosonic Ising model
We first derive the transition rate for the process |k → |k + δk i . To derive the transition rate for |δk i | > 1, we require a generalization of Fermi's golden rule beyond first order in perturbation theory. Such a generalization is provided in [19], considering the case of a multi-level atom transitioning to a continuum. In our case, the bosonic site plays the role of the multi-level atom and the continuum is the phonon reservoir to which the system loses energy. We repeat the final result for an N th order transition here for convenience: |V n,n−1 | 2 h 2 n (E κ ) + s 2 n (E κ ) , , s n−2 (E κ ) = |V n−2,n−1 | 2 s n−1 (E κ ) h 2 n−1 (E κ ) + s 2 n−1 (E κ ) , and E κ is the energy of the continuum state κ, E n is the energy of the nth atomic energy level, ρ(E κ ) is the density of states for the continuum state κ and V nm is the transition matrix element between the nth and mth atomic levels. For N = 1 the transition rate reduces to that of Fermi's golden rule.
Applying this formula to our case, we assume a perturbative transition on the system and energy levels in the bosonic site that are equally spaced The matrix elements are | k i + 1|V |k i | 2 = g 2 (k i + 1)(N − k i ).
We take the parameters α and ξ as phenomenological parameters that may be determined by experiment. For exciton-polaritons the timescale for equilibration is typically of the order of ps [12], while magnon condensation has been reported to occur with thermalization times of 50 ns [24]. For atomic BECs, times of thermalization of the order of tens of ms are typical [25]. For our numerical simulations, we measure all timescales relative to α (or equivalently set α = 1) and choose a small value for ξ , i.e. ξ 1. At thermal equilibrium, the transition rates are equal between the states |k ↔ |k + δk , which ensures that d p k dt = 0. Following [15], we demand that w(k, δk i ) w(k + δk i , −δk i ) = p k+δk i p k . (C.9) The transition rates calculated in (C.7) cancel, and on the rhs we can calculate from the probability distribution at thermal equilibrium Equations (9) are evolved numerically. We use the standard numerical differential equation solver supplied by Mathematica to obtain the results for small boson numbers. For the results presented in figure 3, we used a kinetic Monte Carlo method, which yielded identical results with the differential equation approach for small boson numbers.