A Polariton Graph Simulator
AA Polariton Graph Simulator
Pavlos G. Lagoudakis , and Natalia G. Berloff , , ∗ Skolkovo Institute of Science and Technology Novaya St., 100, Skolkovo 143025,Russian Federation Department of Physics and Astronomy, University of Southampton, Southampton,SO17 1BJ, United Kingdom Department of Applied Mathematics and Theoretical Physics, University ofCambridge, Cambridge CB3 0WA, United KingdomE-mail:
Abstract.
We discuss polariton graphs as a new platform for simulating the classicalXY and Kuramoto models. Polariton condensates can be imprinted into any two-dimensional graph by spatial modulation of the pumping laser. Polariton simulatorshave the potential to reach the global minimum of the XY Hamiltonian in a bottom-upapproach by gradually increasing excitation density to threshold or to study large scalesynchronisation phenomena and dynamical phase transitions when operating above thethreshold. We consider the modelling of polariton graphs using the complex Ginzburg-Landau model and derive analytical solutions for a single condensate, the XY model,two-mode model and the Kuramoto model establishing the relationships between them.
1. Introduction
Engineering a physical system to reproduce a many-body Hamiltonian has been at theheart of Richard P. Feynman’s idea of an analogue Hamiltonian simulator [1]. Suchsimulators could address problems that cannot be solved by a conventional Turingclassical computer, which would allow us to test various models of lattice systems ordiscover new states of matter. Analogue Hamiltonian simulations led to the observationof a superfluid-insulator phase transition in ultracold Bose gases that is closely related tometal-insulator transition in condensed-matter materials [2]. In the last decade variousother systems have been proposed as classical or quantum simulators: ultracold bosonicand fermionic atoms and molecular gases in optical lattices [3, 4, 5, 6], photons [7],trapped ions [8, 9], superconducting q-bits [10], network of optical parametric oscillators(OPOs) [11, 12], and coupled lasers [13] among other systems.The design of an analogue Hamiltonian simulator consists of several importantingredients [14]: (i) mapping of the Hamiltonian of the system to be simulated into theelements of the simulator and the interactions between them; (ii) preparation of thesimulator in a state that is relevant to the physical problem of interest: one could be a r X i v : . [ c ond - m a t . o t h e r] S e p Polariton Graph Simulator i start to condense with thewavefunction ψ i = (cid:113) ρ i ( r ) exp[ iS i ( r )], characterized by the number density ρ i ( r ) anda phase S i ( r ) with the relative phase-configuration that carries the highest polaritonoccupation due to the bosonic stimulation during the condensate formation [28]. Bycontrolling the pumping intensity and profile, the graph geometry and the separationdistance between the lattice sites one can control the couplings between the sites andrealise various phase configurations that minimize the XY model as was shown in[29]. This gives rise to the use of the polariton graph as an analogue XY Hamiltoniansimulator. The search for the global minimum of the XY Hamiltonian is via a bottom-upapproach which has an advantage over classical or quantum annealing techniques, wherethe global ground state is reached through either a transition over metastable excitedstates or via tunnelling between the states in time that depends on the size of the system.Figure 1(a) shows the schematics of the classical thermal annealing, quantum annealingvia tunnelling between the metastable states and the bosonic stimulation.The XY model has been previously simulated by other physical systems: ultra coldatomic optical lattices [30] and coupled photon lasers network [13]. Polariton graphs Polariton Graph Simulator
Energy'landscape Energy'landscape(a) (b)
Classical'Annealing Quantum''tunnelling
Bosonicstimulation Bosonicstimulation above'threshold'
Figure 1. (a) Schematics of different types of annealing for finding the global minimumof the energy landscape of the simulated XY Hamiltonian: classical annealing via thethermal activation, quantum annealing via quantum tunnelling during the adiabaticadjustment of the simulated Hamiltonian and the bosonic stimulation that leads tothe system condensation at the threshold. (b) Schematics of the comparison betweenphoton laser and polariton graph operation.
The interest in using an analogue simulator for finding the global minimum of theXY Hamiltonian is motivated by the recent result in the theory of quantum complexity:there exist universal Hamiltonians such that all other classical spin models can bereproduced within such a model, and certain simple Hamiltonians such as the 2DIsing model with fields on a square grid with only nearest neighbour interactions arealready universal [31]. This suggests that hard computational problems that can beformulated as a universal Hamiltonian can be solved by a simulator that is designedfor finding the global minimum of such Hamiltonian. In particular, the XY model is aquadratic constrained optimization model, which is an NP-hard problem for non-convexand sufficiently dense matrices [32, 29].A more traditional use of analogue Hamiltonian simulators has been based onmodelling electrons as they move on a lattice generated by the periodic array of atomiccores. To elucidate such a behaviour ultra cold atomic condensates were loaded intooptical lattices such as cubic-type lattices [33], superlattice structures [34, 35], triangular Polariton Graph Simulator
2. An approximate analytical solution for a single condensate
The mean field of polariton condensates can be modelled [41, 42] in association withatomic lasers by writing a driven-dissipative Gross-Pitaevskii equation (aka the complexGinzburg-Landau equation (cGLE)) for the condensates wavefunction ψ ( r , t ) : i ¯ h ∂ψ∂t = − ¯ h m (1 − iη d R ) ∇ ψ + U | ψ | ψ + ¯ hg R R ψ + i ¯ h R R R − γ C ) ψ (1)coupled to the rate equation for the density of the hot exciton reservoir, R ( r , t ): ∂ R ∂t = ( γ R + R R | ψ | ) R + P ( r ) . (2)In these equations m is polariton effective mass, U and g R are the strengths of effectivepolariton-polariton and polariton-exciton interactions, respectively, η d is the energyrelaxation coefficient specifying the rate at which gain decreases with increasing energy, R R is the rate at which the reservoir excitons enter the condensate, γ C is the rates of thecondensate losses, γ R is the redistribution rate of reservoir excitons between the differentenergy levels, and P is the pumping into the reservoir. In the limit γ R (cid:29) γ C one canreplace Eq. 2 with the stationary state for the reservoir R = P ( r ) / ( γ R + R R | ψ | ).To non-dimensionalize the cGLE we use ψ → (cid:113) ¯ h / mU (cid:96) Ψ , r → (cid:96) r , t → mt(cid:96) / ¯ h ,where we choose (cid:96) = 1 µ m and introduce the notations g = 2 g R /R R , γ = mγ C (cid:96) / ¯ h , p = m(cid:96) R R P ( r ) / ¯ hγ R , η = η d ¯ h/mR R (cid:96) , and b = R R ¯ h / m(cid:96) γ R U . The dimensionless form of Eqs. (1)-(2) becomes the cGLE with a saturablenonlinearity
Polariton Graph Simulator i ∂ Ψ ∂t = − (1 − iηn R ) ∇ Ψ + | Ψ | Ψ + gn R Ψ + i ( n R − γ ) Ψ , (3) n R = p ( r ) / (1 + b | Ψ | ) . (4)By taking the Taylor expansion for small | Ψ | in the expression for the reservoir n R we arrive at the more standard cGLE i ∂ Ψ ∂t = − (1 − iηp ) ∇ Ψ + (1 − pb ) | Ψ | Ψ + gp Ψ + i (cid:16) p − γ − pb | Ψ | (cid:17) Ψ , (5)where, in the view of smallness of η we dropped gη | Ψ | term. We can compare therelative strength of nonlinearities in Eqs. (3-4) and (5) depending of the physicalquantities that define b. By taking the values the system parameters typically acceptedfor GaAs microcavities [43, 44, 26] ¯ hR R = 0 . · µ m , U = 0 . − . · µ m weobtain b = 2 − /γ R ps . With γ R on the order of 1ps − we have b of the order of the realnonlinearity.We consider the asymptotics and approximations of the steady state solutions for asingle radially symmetric Gaussian pumping profile p ( r ) = p exp( − σr ), where p is themaximum pumping intensity and σ characterises the inverse width of the Gaussian. TheMadelung transformation Ψ = √ ρ exp[ iS ] relates the wavefunction to density ρ = | Ψ | and velocity u = ∇ S . Separating the real and imaginary parts of Eqs. (3-4) we arriveat the mass continuity and the Bernoulli equations:1 rρ d ( rρu ) dr = p ( r )1 + bρ (cid:18) η ( r ( √ ρ ) (cid:48) ) (cid:48) r √ ρ − ηu (cid:19) − γ, (6) µ = − ( √ ρ ) (cid:48)(cid:48) √ ρ − ( √ ρ ) (cid:48) r √ ρ + u + ρ + p ( r )1 + bρ (cid:18) g − ηrρ d ( rρu ) dr (cid:19) . (7)Away from the pumping spot, where p ( r ) = 0, the velocity u = | u | is given by theoutflow wavenumber k c = const with rρ r + ρ = − γρr/k c , which can be integrated toyield ρ ∼ exp[ − rγ/k c ] /r. From Eq. (7) at infinity we obtain µ = k c − γ / k c . In the viewof their asymptotic behaviour the condensate density and velocity can be approximatedby ρ ( r ) = a γr exp( γrk − c ) k − c + ξ − γr/k c + a r , (8)and u ( r ) = k c tanh( lr/k c ) , (9)utilizing their behaviour at the origin and infinity and introducing parameters ξ, a , a ,and l that define the parametric family of solutions. Their values should be found fromthe governing equations via matching asymptotics, as shown below.We neglect η in the view of its smallness and substitute the expressions for thedensity, velocity and the pumping profile into Eqs. (6-7). By expanding the resultingexpressions about r = 0 and setting the term to the order O ( r ) to zero we obtain theequations that define the unknown parameters ξ, a , a , l and k c in terms of the system Polariton Graph Simulator g, b, γ, p , σ . The leading order expansion of Eq. (6) and the first orderexpansion of Eq. (7) fix ξ and a as ξ = a b (2 l + γ ) p − l − γ , a = − γ k c . (10)The expansion to O ( r ) of Eq. (6) determines k c as k c = [4 a bl p (2 l + γ ) + 3 γ (8 l − l ( p − γ ) + ( p − γ ) γ (11)+ 2 l (2 p − p γ + 3 γ )] / (3 a bp σ (2 l + γ ) )]Finally, the expansions to O ( r ) of Eq. (7) define the remaining parameters a and l through two nonlinear equations k c = a ξ + γ ( ξ + 8)4 k c ξ + gp ξa b + ξ , (12) l = a γ k c ξ + 5 γ k c ξ − γ k c ξ − a bgp γ k c ( a b + ξ ) + gp σξa b + ξ . (13)Equations (8) and (9) with the parameters defined by Eqs. (10-13) for the given systemparameters g, b, γ and the pumping parameters p and σ fully specify the approximateanalytical solution of Eqs. (6-7). Figure 2 shows the comparison of the approximateanalytical solutions (solid lines) given by Eqs. (8-9) and the numerical solutions (dashedlines) of Eqs. (6-7) for b = 1 . , γ = 0 . , g = 0 . p = 5 , σ = 0 . p = 10 , σ = 0 .
4. The values specifyingvelocity are k c = 1 . , l = 0 . k c = 1 . , l = 0 . - - - - - - - Distance)( µ m) Distance)( µ m) u u ! ! P P (a) (b)
Figure 2.
Approximate analytical (blue lines) and numerical (black lines) solutionsfor density (solid lines) and velocity (dashed lines) of Eq. (5) for the pumpingprofile given by p ( r ) = p exp( − σr ) (green shaded area). The system parametersare b = 1 . , γ = 0 . , g = 0 . σ = 0 . p = 5; (b) σ = 0 . p = 10. Polariton Graph Simulator
3. Mapping of phases into the classical XY Model
In the previous section we obtained solutions of the governing equation Eq. (5) for asingle pumping Gaussian spot. Spatial light modulator can be used to pump condensatesat the vertices of a distributed graph via p ( r ) = N (cid:88) i =1 p i exp[ − σ i | r − r i | ] , (14)where p i stands for the pumping intensity at the center of the spot at position r = r i .In what follows we shall assume that all vertices are pumped identically, so that p i = p , σ i = σ for all i = 1 , · · · , N . To the leading order and assuming that all condensatesare well-separated we can approximate the resulting condensate wave function, ψ N ,as Ψ N ( r , t ) ≈ (cid:80) Ni =1 Ψ ( | r − r i | ) exp( iθ i ), where Ψ = Ψ ( r ) is the solution of thestationary Eq. (5) for a single localized radially symmetric condensate pumped by p ( r ) = p exp( − σr ) , found in the previous sectionΨ ( r ) = (cid:113) ρ ( r ) exp (cid:34) i k c l log cosh (cid:18) lk c r (cid:19)(cid:35) , (15)where ρ ( r ) is given by Eq. (8). To find the total amount of matter M we write M = (cid:90) | Ψ N | d r = 1(2 π ) (cid:90) | ˜Ψ N ( k ) | d k , (16)˜Ψ N ( k ) = (cid:90) exp( − i k · r )Ψ N ( r ) d r = ˜Ψ ( k ) N (cid:88) i =1 exp( i k · x i + iθ i ) , (17)where ˜Ψ ( k ) = 2 π (cid:82) ∞ Ψ ( r ) J ( kr ) rdr and J is the Bessel function. The total massbecomes M = 2 πN (cid:90) ∞ | Ψ | rdr + (cid:88) i
0, solid lines) or antiferromagnetic ( J ij <
0, dashed lines) by either varying thedistance between the pumping spots, as illustrated in Figs. 3(a-e,g) or by changing thepumping parameters of the individual spots as shown in Fig. 3(f). Spin configuration,such as ferromagnetic, where all J ij > J ij > J ij < J ij < Polariton Graph Simulator Figure 3.
Theoretical prediction of classical magnetic spin configurations showingcontour plots of the normalised density and spin orientations (arrows): (a) ferromagnet,with all J ij >
0, (b) rhomb, with all horizontal J ij >
0, (c) spiral 1, with all J ij < J ij <
0, (e) defect, with a single ferromagnetic coupling(top edge), (f) star, with J ij > | Ψ | , with Ψ = (cid:80) Ψ i . Theindividual wave functions Ψ i are analytic approximations as described in the text with b = 1 . , g = 0 . , p = 3 , σ = 0 . Polariton Graph Simulator J ij , by tuning the distancebetween polariton sites, or the characteristics of the pumping spots (the intensity, p , orthe inverse width of the Gaussian, σ ) leading to more exotic phases. We illustratethe control over an individual J ij on the seven-vertex graphs of Figs. 3(e,f). InFigure 3(e) we utilise control over an individual J ij by tuning the distance between twovertices and introduce a single ferromagnetic coupling (defect edge) into an otherwiseantiferromagnetic configuration. In Figure 3(f) we control the coupling of a singlepolariton site (central vertex) to all its neighbours by tuning the intensity of its pumpingspot and switch them to ferromagnetic in an otherwise antiferromagnetically coupledgraph (star configuration). Finally, polariton graphs allow for fully disordered systemsto be addressed as shown in Fig. 3(g).
4. Kuramoto Model
The cGLE can be reformulated as the Kuramoto model which is a paradigm for aspontaneous emergence of collective synchronization and that has been widely used tounderstand the topological organization of real complex systems from neural networksto power grids [45]. In this context the polariton lattice describes collective dynamics of N coupled phase oscillators with phases θ i ( t ), characterized by the natural frequencies ω i which are associated with the chemical potential of individual condensates.To show this we adopt a two-mode model that neglects the spacial variations andrepresents the network of interacting polariton condensates via the radiative couplings J ij between i th and j th condensates i Ψ it = | Ψ i | Ψ i + (cid:18) ( g + i ) p b | Ψ i | − iγ (cid:19) Ψ i + i (cid:88) j J ij Ψ j . (20)We neglect the blueshift g and without loss of generality let γ = 1. For the densitiesand phases of the individual condensates we obtain12 ˙ ρ i ( t ) = pρ i bρ i − ρ i + (cid:88) J ij √ ρ i ρ j cos θ ij ˙ θ i ( t ) = − ρ i − (cid:88) J ij √ ρ j √ ρ i sin θ ij . (21)The radiative coupling J ij is due to the interference of the condensates from differentpumping spots [46], and are such that J ij (cid:28) ρ i for any j . First, we shall assume thatthe density number dynamics is faster than the phase dynamics, so that the densitiesacquire the instantaneous steady state values ρ i = ρ = ( p − /b to the leading orderin J ij . In this case we get the equations of the phase dynamics represented by theKuramoto model:˙ θ i ( t ) = − ρ − (cid:88) J ij sin θ ij . (22) Polariton Graph Simulator V ( θ ) = ρ N (cid:88) i =1 θ i − N (cid:88) i,j =1 J ij cos θ ij , (23)so that the Kuramoto model (22) describes the gradient flow to the minima of V ( θ )and, therefore, minimizes the XY Hamiltonian.Next, we will allow for the density variations and consider two spots only. For thesystem with just two condensates we introduce the average density R = ( ρ + ρ ) / z = ( ρ − ρ ) / θ = − z − J R √ R − z sin θ , ˙ R = p [( R + z ) Q + + ( R − z ) Q − ] − R + 2 J √ R − z cos θ , ˙ z = p [( R + z ) Q + − ( R − z ) Q − ] − z, (24)where we defined Q ± = (1 + b ( R ± z )) − . Assuming that z and J are small comparedto R we can expand these equations in small parameters z and J and consider thesteady state for R , which to the leading order is R = ( p − /b . Eliminating z fromEq. (24) we see that θ satisfies the second-order differential equation¨ θ + 2 (cid:18) − p + J cos θ (cid:19) ˙ θ = − (cid:18) − p (cid:19) J sin θ . (25)As the pumping increases from the threshold value of p th the oscillations between twocondensates become damped with the rate proportional to 2 (cid:18) − p − + J ij cos θ (cid:19) . Therelative phases lock to 0 or π difference depending on whether J > J < , respectively. This again agrees with the minimization of the XY Model, which for twopumping spots is H XY = −J cos θ with minimum at θ = 0 if J > π if J <
5. Conclusions
In conclusion, we discussed polariton graphs as an analog platform for minimizing theXY Hamiltonian and therefore emulating classical spin model with a potential of solvingcomputationally hard problems. We demonstrated that the search for the global groundstate of a polariton graph is equivalent to the minimisation of the XY Hamiltonian H XY = − (cid:80) J ij cos( θ ij ). The theoretically predicted phase transitions explained the Polariton Graph Simulator
6. References [1] Feynman R R 1982 Simulating physics with computers
Int. J. Theor. Phys. Nature
Advances in Physics Rev.Mod. Phys. Nature
Annu. Rev. Condens.Matter Phys. Nature Photonics Nature
Science
Nature Commun. Opt. Express Nature Photonics Phys. Rev. Lett.
Polariton Graph Simulator [14] Bloch I, Dalibard J and Nascimb´ene S 2012 Quantum simulations with ultracold quantum gases Nature Physics Phys. Rev. Lett. et al. Nature
Science
Contemporary Physics Rev. Mod. Phys. Phys. Rev. B et al Phys. Rev. Lett. et al
Nature
Science
Phys. Rev. Lett. Nature Physics Nature Comm Phys. Rev. X Nature Materials , also arXiv:1607.06065[30] Struck J et al
Science
Science
J.Global Optim. Phys. Rev. Lett. Phys. Rev. A et al Nature et al
New J. Phys. Nature et al
Phys. Rev. Lett.
Phys. Rev. Letts. (2017)[40] Exotic states of matter with polariton chains, in review by Phys. Rev. Letts. (2017)[41] Wouters M and Carusotto I 2007 Excitations in a nonequilibrium Bose-Einstein condensate ofexciton polaritons,
Phys. Rev. Lett. Polariton Graph Simulator [42] Keeling J and Berloff N G 2008 Spontaneous rotating vortex lattices in a pumped decayingcondensate Phys. Rev. Lett.
Phys. Rev. Lett. et al
Nat. Phys. Physics Reports , Phys. Rev. B 85
Phys. Rev. Letts. Physica D: Nonlinear Phenomena R135 (2001)[47] Plumhof J D, St¨oferle T, Mai L, Scherf U and Mahrt R F 2014 Room-temperature Bose-Einsteincondensation of cavity exciton-polaritons in a polymer
Nat. Mat. et al Nature497