Optimization using Bose-Einstein condensation and measurement-feedback circuits

We investigate a computational device that harnesses the effects of Bose-Einstein condensation (BEC) to accelerate the speed of finding the solution of a given optimization problem. Many computationally difficult problems, including NP-complete problems, can be formulated as a ground state search problem. In a BEC, 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. We find that the speed of convergence to the ground state can be sped up by a factor of $ N$ at a given error, where N is the boson number per site.

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]. The 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 phenomena such as Bose-Einstein condensation (BEC) [3]. We show that by using bosonic particles it is possible to speed up the computation of a given optimization problem. The method takes advantage of the fact that bosonic particles tend to concentrate in the minimal energy state at low temperatures. Since many difficult computational problems can be reformulated as an energy minimization problem [4], this is attractive for such computational purposes that a large number of bosons lie in the ground state configuration. The origin of the speedup is due to bosonic final state stimulation, an effect that is familiar from stimulated emission of photons in lasers [5]. This allows the system to move towards the ground state at an accelerated rate.
We formulate the computational problem to be solved as an energy minimization problem of an Ising Hamiltonian [4]. For example, the NP-complete MAX-CUT problem [6], where the task is to group M vertices into two groups A and B such as to maximize the number of connections between the groups, is known to be equivalent to the Hamiltonian H P = ij J ij σ i σ j , where J ij 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 }. In simulated annealing [7], very long annealing times are necessary to ensure that the system does not get caught in local minima. Quantum annealing [8] overcomes such problems due to local minima by introducing a quantum tunneling term but requires a slow adiabatic evolution to prevent leaks into excited states.
The computational device we have in mind is 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 by σ = ±1. Any particle that displays bosonic statistics with an internal spin state may be used, such as exciton-polaritons in semiconductor microcavities, which have recently observed to undergo BEC [9][10][11] or neutral atoms with an unpaired electron in atom chips [12]. Systems that undergo BEC are natural choices for implementation of such a device, since similar principles to the formation of a BEC are required in order for the rapid cooling to the solution of the computational problem.
where S i = N k=1 σ k i is the total spin on each site i, and J ij is the same matrix as in H P which specifies the computational problem. The ground state spin configuration of (1) is equivalent to the original Ising model Hamiltonian H P [21]. This can be seen by noting that the same spectrum as H P is obtained when the site spin is maximized |S i | = N . The energy between these levels connect linearly as the spin on a particular site is changed from S i = −N to N or vice versa.
The interaction Hamiltonian (1) may be produced by measuring the total spin S i on each site, processing those measurement results and feeding an appropriate control signal back into the system by applying a local dc field on site i. For example, say at a particular instant a spin measurement of all the sites are made, giving the result {S j }. Then at that moment a local field B i = j J ij S j is applied on site i, yielding the effective Hamiltonian H = i B i S i . The measurement and the feedback process are continuous. Although J ij has a large connectivity and is long-ranged, by using such a feedback method to induce the interactions there is no restriction to the kind of interactions J ij that can be produced in principle. The above argument can be formulated in the framework of quantum feedback control. We start with the Wiseman-Milburn feedback master equation [15] dρ c /dt = L 0 ρ c + 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 [16].
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 Fig. 1. First consider one particular site i. The meter measures the z-component of the spin, thus we have 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. In order that the system can dissipate energy out of the system we have a dissipation term on each site L 0 ρ c = αD[S − i ]ρ c , where α is a rate constant determining the time scale of the dissipation (cooling), S − i = a † i− a i+ , and a iσ is the annihilation operator for a boson on site i in the state σ = ±1. The first two terms of the master equation thus describe a cooling process with 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 is the current resulting from the measurement of site j, J ij is the same matrix specifying the problem Hamiltonian (1), and Γ is a overall constant. Inserting these expressions into the feedback master equation . Due to the symmetric nature of the J ij matrix, the last term in the above equation can be written [19] as −iΓ [H, ρ], where H is given in equation (1). This gives the time evolution The first term is an evolution of the system according to the Hamiltonian (1), which shows that the feedback Hamiltonian H fb (t) indeed reproduces the desired Hamiltonian (1). The second term is a cooling of the system as before, and the third is a dephasing term originating from the measurement on each site, as well as a contribution from the feedback circuit noise.
Initially each site is prepared with equal populations of σ = ±1 spins, which can be achieved by using a linearly polarized pump laser, in the case of exciton-polaritons. The system is cooled in the presence of the interactions between the sites, by immersing the system in an external heat bath. The readout of the computation is simply performed by measuring the total spin on each site after the system cools down by dissipating heat into the environment. The sign of the total spin gives the information of σ i = ±1 for the original spin model. Since the "computation" here is the cooling process itself, no complicated gate sequence needs to be employed to obtain the ground state.
To understand the effect of using bosons, first compare the thermal equilibrium configuration of a system described above with an equivalent system that uses classical, distinguishable particles. As a simple example, consider the two site Hamiltonian H = −JS 1 S 2 − λN (S 1 + S 2 ), where the second term is included such as there is a unique ground state in spite of the S i ↔ −S i symmetry of the first term in the Hamiltonian. For a single spin on each site and J, λ > 0, the ground state configuration is σ 1 = 1, σ 2 = 1, which we regard as the "solution" of the computational problem. We neglect the presence of an on-site particle interaction ∝ S 2 i here since we assume that the strength of the interactions J produced by the induced feedback method can be made much larger than such a term which may occur naturally.
In Figure 2 we show the average spin on a single site of the two site Hamiltonian, which can be calculated from standard partition function methods accounting for bosonic counting factors. Comparing bosonic particles and classical distinguishable particles, 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 results in an improved signal to noise ratio for the bosons in comparison to distinguishable particles. 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.
We now turn to the time taken to reach thermal equilibrium, after initially preparing the system with equal populations of σ = ±1 particles on each site. We generalize the methods of Ref. [13] to our bosonic Ising model, also accounting for transitions beyond first order in perturbation theory. Given the M -site Hamiltonian where the k i range from 0 to N , a † iσ is the creation operator for a boson on site i in the state σ, and we have defined the vector k = (k 1 , k 2 , . . . , k M ). 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 Mdimensional 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 (1 ± γ) factor to ensure that the system evolves to the correct thermal equilibrium distribution, in a similar way to that discussed in Ref. [13]. We calculate the weight factors to have the form [21] w(k, δk i ) = (1 + γ i (δk)) αξ δk−1 where α is a rate constant determining the overall timescale, γ i (δk) = The F(k, δk) factors are final state stimulation factors due to bosonic statistics, from matrix elements | k + δk|V δk |k | 2 in Fermi's golden rule, where the perturbation causing the transition is V ∝ a † + a − + a † − a + [20]. Transition beyond order one are suppressed by the coefficient ξ 1.
We use the standard numerical differential equation solver supplied by Mathematica to evolve p k for small boson numbers. Figure 3a shows the cooling of the system for N = 1 and N = 50 particles. As the number of particles is increased, we see that the time taken to reach equilibrium is considerably reduced, as well as a high proportion of particles occupying the ground state. For low temperatures, the time dependence of the single site case can be approximated by the rate equations dn 1 dt = − dn 2 dt = α(n 1 + 1)n 2 , where n i are the populations on levels i = 1, 2. Analytically solving this gives a equilibration time of τ ∼ 1/αN for large  Figure   3b, where the time rapidly increases as → 0 (i.e. T → 0) unlike the single site case. In our simulations, we assume that the Hamiltonian (1) is correctly implemented by the feedback scheme, and use a kinetic Monte Carlo method [22] to numerically calculate the cooling time starting from a T = ∞ configuration. A final thermal equilibrium temperature is set, which determines the error probability. Fig. 3b shows that as the boson number is increased, there is a significant speedup at constant error of several orders of magnitude. There is a small odd/even effect due to the definition of the error . The curves approach zero equilibriation time as ∝ 1/N for large N (Fig. 3c). In 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 scheme is also compatible with a thermal annealing procedure, where the temperature is gradually reduced to zero starting from a high temperature configuration. We calculate the residual energy, defined as the average energy above the ground state of the system following the annealing procedure. An exponential annealing schedule with time constant τ 0 is used, starting from a temperature corresponding to an error of = 0.7. Times up to 4τ 0 are annealed where the system no longer responds to the cooling. Fig. 3d shows that the residual energy is suppressed for all N > 1, thus again displaying an improvement due to bosonic final state stimulation.
We conclude that the scheme as shown in Figure 1 has a systematic way of improving the standard Ising model, in terms of a speedup proportional to the number of bosons N per site. 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 ij to ∼ N J ij . 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 emphasize that although the device discussed in this letter is a computational device that uses quantum effects, it is rather different to a quantum computer, since the off-diagonal density matrix elements of the state of the device are explicitly zero at all times. For these reasons we expect the scaling of the equilibration time with the site number M is not faster than exponential, in analogy to the classical case.
The speedup then manifests itself as a suppressed prefactor of this exponential function, which can be accelerated by a factor of N . In its present form, the device can simulate any kind of optimization problem that can be written as an Ising model involving two spins, such as the graph partitioning problem, 2SAT, MAX-2SAT, and others. Extension of the device to involve k-body interactions give a natural implementation of problems such as  We derive an effective Bose-Hubbard model that predicts a phase transition from Bose-Einstein condensate to Mott insulator in two different systems subject to applied periodic potentials: microcavity exciton polaritons and indirect excitons. Starting from a microscopic Hamiltonian of electrons and holes, we derive an effective Bose-Hubbard model for both systems and evaluate the on-site Coulomb interaction U and hopping transition amplitudes t. Experimental parameters required for observing a phase transition between a Bose-Einstein condensate and a Mott insulator are discussed. Our results suggest that strong periodic potentials and polaritons with a very large excitonic component are required for observing the phase transition. The form of the indirect exciton interaction is derived including direct and exchange components of the Coulomb interaction. For indirect excitons, the system crosses over from a Bose-Hubbard model into a double layer Fermi-Hubbard model as a function of increasing bilayer separation. The Fermi-Hubbard model parameters are calculated, and the criteria for the location of this crossover are derived. We conjecture that a crossover between a Bose Mott insulator to a Fermi Mott insulator should occur with increasing bilayer separation.

I. INTRODUCTION
The observation of the Bose-Einstein condensation ͑BEC͒ of exciton polaritons has generated a large amount of interest in recent years. [1][2][3] The focus has now turned to examining various properties of the condensate, such as thermal equilibration, 4 superfluidity, 5 vortex formation, 6,7 and elementary excitations. 8 If the trend followed by atom optics physics community holds for the exciton-polariton community, one important branch of study of exciton-polariton BECs will be the application of periodic potentials on the BEC system. Optical lattices have attracted much attention, spurred on by the experiments demonstrating the phase transition between a BEC and a Mott-insulating state in a Bose-Hubbard model. 9 The application of the periodic potential simultaneously increases the particle-particle interaction, as well as decreasing the kinetic energy. This allows the ratio of the Hubbard on-site interaction to the hopping amplitude U / t to be varied at will. The experiment has been of particular interest in the quantum information community since the experiment realizes a nearly ideal quantum simulator. 10,11 A quantum simulator is a device that directly recreates a quantum many-body problem in the laboratory. By experimentally modifying physical parameters, such as the periodic potential amplitude, temperature, and density, one may explore the phase diagram of the system.
The formation of a Bose-Hubbard model using polaritonic systems was first proposed in Refs. 12 and 13. In this paper we develop the theory for exciton-polaritons subject to a periodic potential ͓see Fig. 1͑a͔͒. In contrast to the works of Refs. 12-17, where the polariton interaction originates from an effective nonlinearity due to a coupling to atomic sites, our interaction originates from the excitonic components of the polaritons, which ultimately originates from a Coulomb interaction. Starting from a microscopic Hamiltonian for electrons and holes and their Coulomb interaction, we derive the origin of the Bose-Hubbard model that is assumed in Ref. 18, allowing an accurate determination of the Bose-Hubbard parameters U and t. From an experimental point of view, steps toward a similar experimental configuration as the optical lattice have been realized already by modifying the semiconductor microcavity system. In Ref. 19, it was shown that a band structure was successfully formed using a metal deposition technique. The periodic metal structure on the surface changes the boundary conditions of the photon field, thus creating a static periodic potential for the polaritons. 20 Other methods for trapping polaritons have been proposed by etching the microcavity. 21 Such etching techniques are anticipated to produce stronger trapping potentials and access a more strongly correlated regime. 14 The formal similarity in the treatment of indirect excitons allows us to write general formulas that capture both the polariton and indirect exciton interaction ͓see Fig. 1͑b͔͒. Polaritons are described in the d = 0 limit of the formulas, where d is the bilayer separation of the indirect exciton system. Indirect excitons have a nonzero d, but a vanishing photon component. Although only exciton polaritons have currently been observed to undergo Bose-Einstein condensation so far, 22,23 there is a large amount of interest in BECs of indirect excitons, as well as works of indirect excitons in periodic lattices, 24 motivating us to write the generalized formulas for both cases. We place particular interest on what parameters are required for observing a Bose-Hubbard Mott transition. For indirect excitons, as the bilayer separation d is increased, the bosonic nature of the excitons gradually diminishes due to the reduced electron-hole interaction. The system is more appropriately described as a Fermi-Hubbard model in this limit. We discuss the criterion for this crossover to occur for our model. We also conjecture that a crossover between a Bose Mott insulator ͑BMI͒ to a Fermi Mott insulator ͑FMI͒ should take place with increasing bilayer separation, in analogy to the more commonly known BCS-BEC crossover. 25,26 SI units are used throughout this paper.

II. BOSE-HUBBARD MODEL
We assume that a periodic potential of the form is applied on the photonic ͑ph͒ and excitonic ͑exc͒ components of the exciton polaritons respectively, where k 0 =2 / and is the wavelength of the periodic potential created. As mentioned in the introduction, a variety of experimental methods exist to create such a potential on the photon field. [19][20][21]27 For the excitonic part, metal gates may be applied to the surface, trapping the excitons under the gates. 28 The potential W exc ͑r͒ is an effective potential for the center of mass motion of the exciton, obtained after integrating over the relative motion of the exciton. The potential can in principle be either type I ͑where the electron and holes share the same potential minimum locations͒ or type II ͑electrons and holes have minima on alternate sublattices͒. For example, deformation potentials induced by surface acoustic waves 29 and dipolar traps 28 are type I potentials.
Meanwhile, the piezoelectric trapping technique 30 is an example of a type II trapping potential. Type II potentials need rather strong trapping potentials for each individual component of the exciton ͑electron and hole͒ since their effective amplitude is suppressed by a factor of ͑k 0 a B ͒ 2 ͓see Eq. ͑30͒ in Ref. 30͔. However, since for strong potentials type II potentials tend to ionize the excitons, we believe that a type I potential is more promising in order to avoid these undesired effects.
The total Hamiltonian of the system is then are the annihilation operators for the quantum well ͑QW͒ excitons and microcavity photons respectively, M = m e + m h is the exciton mass, g is the exciton-photon coupling, and U exc ͑r , rЈ͒ is the effective interaction between two excitons. The photon acquires an effective mass m ph through the dispersion in a two-dimensional ͑2D͒ microcavity, where the photon energy is E ph = m ph c 2 , and c is the speed of light in GaAs. We do not consider the spin of the excitons explicitly because we assume that the polaritons are injected with a linear polarization such that only one spin species is present. The form of the exciton interaction is discussed in detail in Sec. III A. H sat is the nonlinear interaction due to the excitonphoton coupling 31 and is discussed in Sec. III B. Substituting Eqs. ͑8͒ and ͑9͒ into Eq. ͑2͒ and defining the upper ͑ = ↑͒ and lower ͑ = ↓͒ polariton operators, where ͉u ͉ 2 + ͉v ͉ 2 = 1. We obtain the polariton Hamiltonian where we have only included terms where lower polariton operators appear, and compacted the notation such that p q ϵ p q ↓ . Physically, disregarding the upper polariton operators corresponds to a low-temperature regime where there is negligible upper polariton population, which is routinely achieved experimentally. The Hopfield coefficients u ϵ u ↓ and v ϵ v ↓ are taken around q = 0, again assuming that only the low energy states are excited. The lower polariton dispersion is obtained by expanding around q = 0 giving where the lower polariton mass is given by Reverting to real space makes it clear that we have polaritons in a periodic potential For sufficiently low temperatures, we may retain the lowest-energy band of Hamiltonian ͑18͒ to a good approximation. A necessary temperature criterion is that the thermal energy k B T is less than the band gap ⌬. In one dimension, any nonzero potential W 0 will open a band gap. In two dimensions however, for small potentials the lowest-energy band overlaps in energy with the second lowest-energy band. To ensure that these bands are separated, a potential of approximately W 0 Ϸ ប 2 k 0 2 / 2m pol is needed, 32 where is the total potential amplitude due to exciton and photon parts. Under these circumstances, we may make a Wannier transformation and retain only states in the lowest-energy band. This yields H = ͚ n,nЈ t͑n,nЈ͒p n † p n Ј + 1 2 ͚ n 1 ,n 2 ,n 3 ,n 4 U͑n 1 ,n 2 ,n 3 ,n 4 ͒p n 1 † p n 2 † p n 3 p n 4 , ͑20͒ where

͑23͒
w͑r − n͒ is the Wannier function centered around the lattice point n = ͑n x , n y ͒. The Hamiltonian ͑20͒ is a Bose-Hubbard Hamiltonian. We shall be only concerned with the nearestneighbor tunneling matrix elements t and the on-site Coulomb interaction U in this paper: and U ϵ U͑n,n,n,n͒ = ͵ d 2 rd 2 rЈ͉w͑r͉͒ 2 U pol ͑r,rЈ͉͒w͑rЈ͉͒ 2 .

͑25͒
In writing Eq. ͑20͒ we assume that the polariton lifetime is sufficiently long such that the extended "superfluid" and Mott states can occur. For example, for an extended superfluid state we require that there is enough time for the polariton to hop several times before decaying Similarly for the Mott-insulating state, the Coulomb energy should obey In addition to the band-gap criterion k B T Ͻ⌬, it is also necessary to have k B T Ͻ U , t, such that only the low-energy physics of the Bose-Hubbard model is probed. To evaluate Eq. ͑25͒ we require a form for the polariton-polariton interaction which we evaluate in the next section.

A. Exciton-exciton interaction contribution
To obtain the effective interaction between excitons, we use the methods of de-Leon and Laikhtman. 33 There are many approaches to find the effective interaction between excitons in the literature, such as usage of the Usui transformation, 34 variational wave function methods, 35 and operator methods. 36 We find that the wave function methods of Ref. 33 are most transparent and systematically give the exciton interactions to order a B 2 / A, where a B =4⑀ប 2 / 2e 2 is the 2D Bohr radius, A is the trapping area of the excitons, and is the reduced mass = m e m h / ͑m e + m h ͒. We note that a similar method was used by Ciuti et al. 37 to find the same result for the Coulomb interactions, but it is unclear how to treat corrections to the kinetic-energy operator ͑"kinematic corrections"͒ based solely on their method. Ref. 33 makes it clear that such terms cancel in the end and do not give rise to a physical interaction. Here, we generalize the results to indirect excitons in a periodic potential.
Following Ref. 33, the effective Hamiltonian for the excitons can be decomposed into the following terms ͑omitting function labels for brevity͒: where the terms are the direct exciton scattering, the exciton exchange scattering, the electron exchange scattering, the hole exchange scattering, the correction due to nonorthonormality of the wave functions, and the contribution due to excited states of the excitons. Since only exciton-exciton interactions to order a B 2 / A are kept, only two-body exciton interactions need to be considered. Explicit expressions for the above terms are as follows. The direct term is where Hamiltonian for the two exciton system is with V͑r͒ = e 2 / 4⑀r ͑⑀ Ϸ 13⑀ 0 is the permittivity in GaAs͒. The exciton exchange term is The electron ͑hole͒ exchange terms are obtained by multiplying Eq. The evaluations of the direct and exchange terms are deferred to Appendix A. We find the direct term to be where E 1s is the binding energy of a 1s exciton. The function I dir ͑q , d͒ is plotted for various d in Fig. 2. The exciton exchange term is where ⌬Q = ͉QЈ − Q͉ and is the angle between QЈ − Q and q.

͑39͒
where in the square brackets we evaluated half of the operator on the initial states and half on the final states. Doing this we see that these terms exactly cancel with the corrections due to nonorthonomality ͓i.e., the fifth term in Eq. ͑28͔͒. Numerical evaluations of the exchange integral I exch are shown in Fig. 3. Substituting Eqs. ͑36͒-͑39͒ into Eq. ͑28͒ we obtain the final effective Hamiltonian for the two-exciton system. Subtracting the kinetic-energy and binding-energy terms, we obtain an expression for exciton-exciton interaction The dependence on d for small q may be evaluated exactly for the direct term by expanding the term in the square brackets in Eq. ͑A1͒. We obtain This is the zero momentum limit of the Fourier transform of the Coulomb interaction for oriented dipoles U͑q͒ = e 2 ͓1 − exp͑−dq͔͒ / ⑀q. The d dependence of the exchange term must be evaluated numerically, and our results are shown in Fig. 4. We find an approximately linear dependence of the exchange term with d, which changes sign at d / a B Ϸ 0.66. In Fig. 4 we also plot the combined contributions of the direct and exchange integrals. We see that the total interaction remains repulsive for all d despite the exchange term changing sign.

B. Saturation contribution
The saturation contribution to the polariton interaction comes from the coupling of the electron and holes to the electromagnetic field The last two terms in Eq. ͑2͒ may be found by considering the matrix element between the states Starting from the two exciton wave function, H EM can either destroy one of the excitons and produce a photon or take an electron and hole from each of the excitons and produce a photon. These two processes give rise to the matrix element and we have set QЉ = Q + QЈ − q due to momentum conservation. The first two terms in Eq. ͑47͒ correspond to the destruction of an exciton to create a photon, with another exciton acting as a bystander. This corresponds to the second last term in Eq. ͑2͒. The last two terms correspond to an electron and hole being taken from each exciton, resulting in a new exciton being formed from the remaining electron and hole. This process clearly requires two excitons in the initial state, giving the nonlinear last term in Eq. ͑2͒. We only consider the case of zero bilayer separation ͑d =0͒ here since the Hamiltonian ͑42͒ requires that the electron and hole recombine into a photon at the same spatial position. In our approximation ͓Eq. ͑34͔͒ where the electron and hole wave functions are perfectly confined to their respective quantum wells for nonzero bilayer separation there is zero overlap of the electron and hole wave function, which gives a zero matrix element for Eq. ͑47͒. We thus find  where I sat ͑Q , QЈ , q͒ is evaluated in Appendix A numerical evaluation of the integral as a function of the photon momentum is shown in Fig. 3.

IV. MOTT-HUBBARD TRANSITION
Returning to the effective Bose-Hubbard Hamiltonian of Eq. ͑20͒, we may now estimate the size of the tunneling and Coulomb energies from Eqs. ͑24͒ and ͑25͒. Starting from Eq. ͑25͒ we make a change in variables R CM = R + RЈ and = ͑R − RЈ͒ / 2, giving

͑51͒
From Figs. 2 and 3 we see that interaction is large up to a momentum of the order of ϳ1 / a B . Thus the largest contributions to the above integral occur when the variable is of the order of ϳa B , which is much smaller than the length scale of the Wannier functions ϳ. We may thus write Figure 5 shows the two contributions to U as well as the hopping amplitude t. The results are normalized to the characteristic energies where we have used Eq. ͑49͒ to convert G into g. Here, 2បg is the Rabi splitting of the polaritons. For polaritons, the bilayer separation is d = 0, while for indirect excitons the exciton component is ͉u͉ 2 = 1. We see that increasing the potential strength W 0 decreases the hopping t while increasing U as expected. Increasing d enhances U coul as the dipole moment of the excitons are enhanced with an increasing d.
We now derive a criterion for a quantum phase transition from a BEC state into a Mott insulator state. In two dimensions, the phase transition is expected to occur at approximately 40 By turning up the potential W 0 it is clear that at some point U / t will reach this critical amplitude. The other variable that may be changed to reach the phase transition is the detuning of the polaritons which changes the polariton mass. For a potential of size W 0 Ϸ ប 2 k 0 2 / 2m pol , we may derive a criterion for the polariton mass necessary to reach the phase transition using Eq. ͑57͒ and the ratio of the dimensionless parameters in Fig. 5. For GaAs, this gives where in the second relation we assumed that ͉u͉ 2 is of the order of unity and a B Ϸ 10 nm. This corresponds to extremely far blue-detuned polaritons, i.e., very excitonlike polaritons. The lack of the dependence on the wavelength is due to the cancellation of the dependence of the Coulomb interaction energy U 0 and the kinetic energy t 0 . The wavelength is still important however as it sets the energy scale of the whole Bose-Hubbard Hamiltonian. The energy scale should be set such that the parameters U and t are larger than the temperature of the experiment, such that only the lowestenergy band is occupied. Furthermore, semiconductor systems possess an inevitable disorder potential due to reasons such as crystal imperfection and damage during fabrication, thus should be set small enough such that the hopping energy t overcomes this disorder potential strength. For larger potentials than W 0 Ϸ ប 2 k 0 2 / 2m pol , a lighter polariton mass is allowable for reaching the phase transition. Thus in practice a combination of blue-detuned polaritons and large potentials is probably the most favorable experimental configuration. For example, using typical experimental parameters for polaritons ͑d =0͒ in GaAs using the criterion ͓Eq. ͑58͔͒, corresponding to u Ϸ 0.999, with = 0.5 m, a B = 10 nm, 2បg = 15 meV, and an applied potential of W 0 = 6 meV, we obtain U = 0.24 meV and t =9 eV. This corresponds to temperatures in the vicinity of T Ϸ 0.1 K, which are reachable using today's refrigeration methods. In order that the system is stable in the Mott-insulating state, the lifetime of the polaritons should be longer than the time scale set by Eq. ͑27͒. Assuming that the lifetime of the very excitonlike polaritons can be approximated by typical exciton lifetimes Ϸ 1 ns, 41 this corresponds to an energy scale U Ͼប/ Ϸ 1 eV. We thus see that for the above parameters the lifetime requirement is satisfied.
The Coulomb interaction is increased for indirect excitons as shown in Fig. 5. However, the increase is fairly modest for bilayer separations of the order of the Bohr radius. Thus considering that indirect excitons have not been observed to undergo BEC yet, the moderate advantage of increased interactions ͑at the sacrifice of a lighter polariton mass͒ is outweighed by the difficulty of cooling the system into the ground state.
The state of the polaritons may be measured using standard photoluminescence measurements that measure the coherence across the condensate. 42 In analogy to the experiments of Greiner et al., 9 the transition to the Mott insulator state should lead to a lack of spatial coherence across the sample, resulting in the destruction of the far-field interference pattern. However, the disappearance of interference fringes does not unambiguously demonstrate the presence of a Mott insulator since an uncondensed state will also have the same interference characteristics. Therefore, a second order coherence Hanbury-Brown-Twiss measurement is also necessary to determine the correlations between the photons emanating from the sample. At unit filling in a Mott insulator state, the conditional probability of detection of a photon originating from a particular site following detection from the same site vanishes. A similar diminished probability is present at higher filling factors. This is the identical technique used to observe anti-bunching behavior in single photon generation.

V. FERMI-BOSE CROSSOVER BOUNDARY
We now turn to the effect on the particle statistics of indirect excitons as the bilayer separation is increased. For d = 0, excitons are well approximated by bosons for sufficiently low density, which motivated us to describe the system as a Bose-Hubbard model in Sec. II. In the limit d = ϱ, excitons cannot be described by bosons and are more properly described as a double Fermi-Hubbard model, with electron-hole interactions between them. Since the two descriptions are rather different with differing phase transition criteria, it is of interest to know at what d this crossover occurs. The bosonic nature of ͑or the lack of͒ the excitons may be seen by examining the commutation relation We have omitted the label n and written b n † → b † to simplify the notation. Only the commutation relation in the same potential minima of the periodic potential is considered here ͑i.e., the same Wannier function͒, as deviations from bosonic behavior should be most apparent in this case. The operator K may be interpreted as the correction operator to the commutation relation ͑59͒ as it contains the nonbosonic component of the exciton operators.
Due to the presence of the K operator, n-particle states defined using the b † operators defined in Eq. ͑60͒ do not have the simple 1 / ͱ n! normalization of bosonic states. We must instead define such states according to where N͑n͒ is present for proper normalization. A derivation of this normalization factor is given in Appendix B, up to powers linear in the correction operator K. We obtain where I BF ͑W 0 , d͒ is the integral expression given by Eq. ͑B8͒ and has an order of magnitude ϳa B 2 / 2 . All terms neglected in Eq. ͑63͒ have higher powers of a B 2 / 2 , which have a small contribution for the typical periodic potential dimensions that are possible using current fabrication methods ͑ӷa B ͒.
Defined in the way Eq. ͑62͒, the states ͉n͘ provide an orthonormal basis set. Deviations from bosonic behavior occur due to the operator b † not providing the correct mapping between these states ͗n͉ b † ͱ n ͉n −1͘ 1. Using the definition ͓Eq. ͑62͔͒, we have Deviations from unity of the right-hand side represents nonbosonic behavior. Substituting Eq. ͑63͒, this factor is to lowest order in I BF ͑W 0 , d͒ ͱ N͑n͒ Since the Bose to Fermi transition is a smooth crossover, 25 strictly speaking it is arbitrary where to mark the boundary.
However, a reasonable criterion for the location of the crossover from bosonic to fermionic behavior may be defined as when the second term in the above expression becomes of the order of unity: The solution to the above criterion is plotted in Fig. 6 in the space ͑d , W 0 ͒. We see that with decreasing n, d, and W 0 , the excitons become more bosonlike. This dependence on n and d is a restatement of the well-known result that the excitons become nonbosonic when their wave functions start overlapping, i.e., na B 2 / A ϳ 1. The dependence on W 0 may be understood by considering the spread of the Wannier functions with W 0 . As W 0 is increased, the Wannier functions become more localized, effectively reducing the area that the excitons are confined in. This enhances the overlap between the excitons, thus pushing the boundary toward the fermion side of the crossover. For n Յ 1, there is no solution to Eq. ͑66͒, meaning that the boundary for bosonic behavior extends all the way to infinity in W 0 and d. The fact that solutions for bosonic behavior exist with n Ͼ 1 means that in reality one cannot treat the excitons completely as hard-or soft-core bosons, and their true nature lies in between these two limits.

VI. FERMIONIC DESCRIPTION OF THE BILAYER SYSTEM
For parameter regions where the bosonic approximation is invalid, we must write the Hamiltonian in its full form involving both electron e ͑r͒ and hole h ͑r͒ operators: with i = e , h and ␦ ih is a Kronecker delta. The Wannier functions w i ͑r͒ differ for electrons and holes due to their different masses. There is a lattice offset of ͑1,1͒/2 since we assume a type II potential, i.e., the potential minima locations for electrons and holes differ by half a lattice unit. A minimal approximation to Eq. ͑69͒ is to retain the nearestneighbor terms in Eq. ͑70͒ and on-site terms in Eq. ͑71͒. Figure 7 shows the results of our numerical evaluations of t i = t i (͑n x , n y ͒ , ͑n x +1,n y ͒) = t i (͑n x , n y ͒ , ͑n x , n y +1͒) and U ij = U ij ͑n , n , n , n͒. In a similar way to Bose-Hubbard parameters of Fig. 5, the application of the periodic potential W e ͑r͒ acts to increase the electron-electron and hole-hole interaction and decrease the hopping amplitude. The electron-hole interaction plateaus off since potential minima of the two particle species sit on two spatially separate sublattices. Comparison with Ref. 43 reveals that for U ij / ͑t e + t h ͒ ӷ 1 and at a density of one exciton per site ͑half-filling in the terminology of Ref. 43͒, the excitons will be in a Mott-insulating regime in both the electron and hole layers. Thus for a large enough potential W 0 e the system will lie in such a Mottinsulating phase. We again assume that the lifetimes of the indirect excitons ͑which can exceed ϳs according to Ref. 45͒ should exceed the requirements given in Eqs. ͑26͒ and ͑27͒ in the respective phases.
Examining various limits leads us to draw a qualitative phase diagram as shown in Fig. 8. First consider traveling up the d axis, with W 0 = 0. Assuming a periodic potential with = 0.1 m, one exciton per potential minima corresponds to a density of n exc =10 10 cm −2 . Monte Carlo calculations have predicted that Wigner crystallization should occur at r s Ϸ 37, 46 corresponding to a density of n WC = 2.3 ϫ 10 8 cm −2 in GaAs. As n exc is above the Wigner crystal melting density n WC , we expect that the bilayer system should be conducting ͑i.e., a metallic phase͒ for d → ϱ. For d = 0, the system is still in a fairly low-density regime ͑a B 2 / 2 Ӷ 1͒, and thus we expect that the ground state may be described by a BEC ͑i.e., a nonlocalized metallic state͒ for sufficiently low temperatures. As d is increased, the Bohr radius of the excitons increase, until the exciton wave functions start to overlap. Beyond this point, the excitons cannot be described as bosons anymore, and the system enters a BCS phase. 25,26 Moving in the direction of increasing W 0 for small d, as discussed in Sec. IV, we expect a Bose-Hubbard transition into a Mott-insulating phase. From the considerations of Ref. 43, at unit exciton density we expect the system to be in a Mott-insulating phase for U ee Ͼ U eh , U hh Ͼ U eh , and U ij / ͑t e + t h ͒ ӷ 1. We thus expect that a transition should occur from the electron-hole plasma phase to a Mottinsulating phase for large d. Connecting the two boundaries for small and large d leads to the phase diagram Fig. 8. It is plausible to expect that the Bose and Fermi Mott-insulating states can be smoothly connected, in a similar way to a BEC-BCS crossover. 25,26 We thus conjecture that the first order transition line between the metallic and Mott-insulating states can also be smoothly connected throughout the phase diagram. The repulsion between the particle species generally increases with increasing d, as can be seen in Fig. 4. Thus qualitatively the transition should shift to smaller values of W 0 for the fermionic limit, as shown in Fig. 8.

VII. SUMMARY AND CONCLUSIONS
We have considered the effect of applying a periodic potential on interacting exciton polaritons and indirect excitons. Our main result is shown in Fig. 5, where the Bose-Hubbard parameters for the on-site interaction U and the tunneling amplitude t was calculated. We also derived a guideline ͓Eq. ͑58͔͒ for the range of parameters necessary to realize a phase transition from a BEC phase into a Mott-insulating phase. The results suggest that very excitonlike polaritons are required to observe the transition. Loosely speaking, the reason is that for the typical experimental parameters, the tunneling amplitude t is far greater than the interaction energy U. Thus in order to make these parameters on the same order, the polariton mass needs to be increased to reduce t. This results in the necessity of polaritons with a large exciton component. Alternatively, a very large potential amplitude W 0 can be applied. The experimental challenge in this case is to maintain U and t greater than the experimental temperature and system disorder. Since the energy scale of the Hubbard parameters are set by the applied potential period, this favors small in order to increase the energy scale. Although we focused mainly on parameters for GaAs, we note our formulas are general enough such that a simple substitution of material and geometrical parameters in Fig. 5 should be enough to find the Hubbard parameters for any semiconductor system.
We have also considered the effect of increasing the bilayer separation for indirect excitons, where there is a crossover from a Bose-Hubbard model to a double Fermi-Hubbard model. The Hubbard parameters for the fermionic limit were derived ͑Fig. 7͒. A Mott transition should be present for both limits, thus we argue that there should be a transition for all intermediate d.
In an analogous way that there is a BEC-BCS crossover for zero potential, 25,26 the Mott-insulating limit should also crossover from a Bose Mott insulator to a double Fermi Mott insulator for large potentials. Our argument is based on connecting the various limits of the system and requires a more rigorous numerical investigation to confirm our conjecture. A more detailed investigation of the various phases would require an extensive numeric survey of the parameter space, which we leave as future work.

ACKNOWLEDGMENTS
This work is supported by the Special Coordination Funds for Promoting Science and Technology, Navy/SPAWAR Grant No. N66001-09-1-2024 and MEXT. P.R. would like to acknowledge financial support from the German Research Foundation ͑DFG͒ via Tr950/1-1 and Re2978/1-1. T.B. and P.R. thank Hui Deng, David Press, and Sven Höfling for valuable comments regarding the paper.

͑A3͒
and J 0 ͑x͒ is the Bessel function of the first kind. The normalization assuming the electrons and holes are confined as delta-functions in the z direction is The electron and hole exchange terms may be obtained by following the derivation given in the Appendix B of Ref. 37. We obtain Eqs. ͑38͒ and ͑39͒, where

͑A7͒
The dimensionless integral appearing in Eq. ͑50͒ is

͑A8͒
where ͑Ј͒ is the angle between q and Q͑QЈ͒. where we have used the fact that the relative wave function ͑r͒ extends out to a distance of the order of ϳa B , while the Wannier function extends out to at a distance ϳ, with ӷa B . The Wannier function has dimensions of the inverse length ͑in 2D͒, hence the order of magnitude of the first integral is ϳ1 / 2 . The order of magnitude of the second integral is ϳa B 2 , making the whole integral of the order of ϳa B 2 / 2 . As can be shown by direct calcu-lation, integrals involving higher powers in the operator K involve higher powers of a B 2 / 2 . Therefore, the approximation made in Eq. ͑B5͒ is thus reasonable as long as a B 2 Ӷ 2 .