Effects of the non-parabolic kinetic energy on non-equilibrium polariton condensates

In the study of non-equilibrium polariton condensates it is usually assumed that the dispersion relation of polaritons is parabolic in nature. We show that considering the true non-parabolic kinetic energy of polaritons leads to significant changes in the behaviour of the condensate due to the curvature of the dispersion relation and the possibility of transfer of energy to high wavenumber components in the condensate spatial profile. We present explicit solutions for plane waves and linear excitations, and identify the differences in the theoretical predictions between the parabolic and non-parabolic mean-field models, showing the possibility of symmetry breaking in the latter. We then consider the evolution of wavepackets and show that self-localisation effects may be observed due to the curvature of the dispersion relation. Finally, we revisit the dynamics of dark soliton trains and show that additional localized density excitations may emerge in the dynamics due to the excitation of high frequency components, mimicking the appearance of near-bright solitary waves over short timescales.

saturation interaction is strongest when the polariton consists equally of photon and exciton components. We assume that the scattering length is much shorter than the de Broglie wavelength, and we do not preserve the effect of Hopfield coefficients on the interaction as it would take into account the decrease of both Coulomb and saturation effects for large exchanged momenta 4 . However we do include the polarisation of polaritons, so the polariton condensate spinor wave function ψ = (ψ + , ψ − ) for the lowest energy mode is governed by two coupled complex Ginzburg-Landau-type partial differential equations 2,3,5,18 , which take into account the effects of weakly short-range and polariton-reservoir interactions, as well as non-equilibrium properties such as: incoherent gain and decay of condensate polaritons; energy relaxation 19 ; and spin/polarisation 12,13,20 . Note that whenever we consider the spin coherent case ψ + = ψ − we use the notation ψ for the condensate wave function. Due to a constant loss of information, corresponding to the rapid decay of polaritons, the condensate is effectively semiclassical (as observed experimentally e.g. in refs 12-14). For the more general quantum treatment and implications however we refer to refs 4, 15, 16. In addition we take into account the fact that the kinetic energy is non-parabolic 2 , i.e. it varies due to the k dependence of the Hopfield coefficents. We define the dispersion relation (or equivalently, kinetic energy density), of the lower (upper) polariton branch as 2,4 , R L,U c av exc c av exc 2 2 Note that the lower and upper polariton modes can be distinctively addressed experimentally due to the large energy gap between the two branches (see Fig. 1). Here the dispersion of the cavity photon is ω = + k q k ( ) c n z cav 2 2 0 . We have used the notation k = |k| for  ∈ k D with D = 1, 2, with c denoting the speed of light, = π q z M l 2 z the quantization of the confined photon in the z-direction, M the number of the quantized z-mode orthogonal to the k-plane, n 0 the refractive index between the cavity mirrors, and l z the cavity spacer length. The dispersion of the exciton ω ω ω ≈ + ≈ − k k ( ) 10 exc e xc 0 42 exc 0 can be assumed to be constant close to the centre of the polariton dispersion. The minimum splitting between the two dispersions, which is obtained at ω cav (k) = ω exc (k), is given by 2Ω R . We set ω = .
1 557 exc 0 eV, the mass of the cavity photon m cav ~ 10 −4 -10 −5 m e and the effective exciton mass m exc ~ 0.1 − 1m e with m e the electron mass, reflecting the experimental values used in refs 12, 13. While we reduce our considerations to the mean-field states we note for the sake of completeness that the polariton can be described quantum mechanically by the Hamiltonian 15 where the kinetic energy of the polariton mode generated by â k and the potential energy due to elastic polariton-polariton scattering are given by k k k k k p k k p k k k p k p kin L ,U p p , , , E b and a B are the exciton binding energy and the Bohr radius, and S is the system area. The incoherent processes for interactions of polaritons with a thermal bath of acoustic phonons and the effect of incoherent pumping in the rotating wave approximation is given by ref. 15 H H with b q being the phonon destruction operators and η d the operators corresponding to the bosonic pumping reservoir, while we refer to ref. 15 for further details about the functions L z , G q and η ⁎ g k and limitations of this approach. Now supposing the above assumptions and the quantum state to be in a product state and that a complex mode ψ carries all relevant information of the many-body system we directly obtain a mean-field partial differential equation 2 . Here in particular the kinetic energy becomes Scientific RepoRts | 7: 1891 | DOI:10.1038/s41598-017-01113-8 kin L ,U 2 whose nature and implications are the topic of this paper. The interaction terms are very small under state-of-the-art experimental conditions and the interactions between the condensed phase and the environment and pumping are modelled in polariton mean-field theory phenomenologically 2,18 . While the kinetic energy density is the "complete" factor for calculating the kinetic energy of the polariton condensate, the effective mass has played a significant role in recent publications. Next we clarify its nature.

Results
Approximations to the kinetic energy. By looking at the kinetic energy density, the effective mass m(k) of particles in semiconductors can be seen to depend on the curvature of the particles' dispersion relation ω(k) at a particular k-value 7,21,22 : where f(k)′′ means second derivative of f in k (i.e. curvature with respect to k). At the simplest level of approximation, the curvature of the dispersion relation is taken to be fixed, as is the case for a free particle with its associated parabolic dispersion. For the case of a polariton condensate this approximation is valid for small values of k. This approximation is widely used in polariton condensate modelling, leading to a second derivative/Laplacian for the dispersion in the equations of motion. The next level of approximation allows for the effective mass to vary with k, i.e. a "velocity dependent effective mass", as has been used for instance in ref. 7. This approach effectively reduces the dispersion relation to the quadratic term in the Taylor series expansion of the full dispersion, with some constant offset. Finally, the third approach, and the one we focus on, is to use the full dispersion in the equations of motion. Note that while we consider here a single site, in contrast to lattices of polariton condensates, the semiconductor itself possesses a lattice structure which acts on the electrons and holes, i.e. the effective mass is a product of the polariton mean-field model. Due to this effective mass the kinetic energy of a polariton condensate varies with k non-parabolically. This leads to the approximate approach to the mean-field description considering the kinetic energy to be where ω i is an energy offset at the bottom of the polariton dispersion, i.e. , which corresponds to the dashed lines in Fig. 1. In case this variation of the effective mass is further neglected we arrive at the simplest parabolic form for the kinetic energy [3][4][5]  which corresponds to cGPE theory 2 and here the kinetic energy density is the dotted line in Fig. 1. Next we address the question how these approximations are connected with the general dispersion relation. For the condensate wave function ψ(k) in k-space it is given by kin L ,U 2 The full polariton dispersion corresponds to the continuous line in Fig. 1 associated with the upper and lower branch respectively. Furthermore by applying Taylor's theorem to ω k ( ) corresponds to the inertial mass which determines the wavepacket velocity from de Broglie's relation introduced in ref. 8. It is zero when expanding at a = 0. By dropping the remainder term in (10) we obtain the parabolic approximation (8). In addition by Taylor's theorem we find that the k dependent effective mass becomes a reasonable concept when expanding the dispersion for k close to zero where ω′ = (0) 0 L,U and when ∫ In Fig. 1 we show the deviation between the approximate parabolic kinetic energy density, the kinetic energy density with velocity dependent effective mass and the complete kinetic energy density in terms of their k dependence, indicating their range of validity. Figure 1(i) shows that the effective mass introduced in ref. 8 approximates the kinetic energy density from below and closer than the (parabolic) cGPE for the lower polariton energy branch. It overrates the kinetic energy shift for the upper branch as seen in (ii) with crossing of the curves at µ . However as the full dispersion and the effective mass concept allow larger k's to be occupied at lower energy in the lower branch widening of the wave packet in k space is expected and thus, due to the Fourier transform's duality properties, highly localised features in the wave packet are expected to emerge. We now look at how these different kinetic energy density approximations affect the resultant condensate model.
Spin sensitive state equations. The kinetic energy term of a polariton condensate guiding equation may be defined in terms of a Fourier multiplier outside the Fourier transform of the condensate wavefunction 7 , i.e. in the form where ω L,U (k) is real-valued and is the Fourier transform of q L,U iff the transform exists and ☆ denotes a convolution between two functions. Eq. (11) corresponds to a kinetic energy (9). We consider this form of the kinetic energy for the spinor polariton field ψ componentwise, which is otherwise governed by cGPEs. These cGPEs are often coupled to a rate equation for the excitonic reservoir n R 5, 23, 24 . Thus the two components' coupled wave equations are given by In (12) we consider a number of physical parameters: α 1 > 0 is the repulsive self-interaction strength, α 2 = −0.1α 1 the attractive cross-interaction strength, γ C gives the scattering rate of the reservoir into the condensate, Γ d the decay rate of condensed polaritons and η approximates the additional energy relaxation processes 14 .
In general we also have an external time dependent potential V ext (r, t). The reservoir dynamics are usually given by the rate equation 2 pump pump denotes the pumping distribution associated with the incoherent scattering into the polariton BEC of component±, and Γ R denotes the reservoir decay rate. We note that a formally equivalent equation has been suggested earlier for the atom laser based on atomic BEC 25 . For fast reservoir relaxation Γ R ≫ Γ d the reservoir dynamics are much faster than that of the condensate and the reservoir population can be approximated to leading order as 14, 23 , which for small amplitudes |ψ ± | 2 can be further simplified.
Consequently the growth and decay terms in (12) can be written in the form refs 2, 23 so that 2 and γ = Γ d /2 identified accordingly (while neglecting the relaxation contribution to the density occupation). Now to compare (12) with the corresponding effective mass approximations we generally write 1 where | | q k ( ) corresponds to one of the three cases: as pointed out in the following discussions.
Unpolarised condensate. When the condensate is unpolarised one considers a simple single component PDE to govern the condensate wave function: which is the spin coherent counterpart of ± n R and its approximation physically justified on the same grounds. One can further simplifying the reservoir and decay of polaritons via a Taylor series given by ref. 23 In addition, for small amplitudes of the wave function, one may approximate η η β Unpolarised plane waves. After discussing the explicit guiding equations we now turn to results and predictions which deviate from the simpler approximations to the polariton kinetic energy (9). We begin by considering the stationary (i.e. time-independent) plane wave solutions of polariton condensates in the presence of pumping and decay processes, i.e. solutions of the equation To solve this model we take the stationary ansatz Consequently we get an algebraic equation for the plane wave amplitudes The complete solution is therefore given by . In Fig. 2(a,b) we present the 1 d plane wave density ρ(k i ) = |ψ(r, k i )| 2 as a function of the dispersion relation indicating the inherent differences of the kinetic energy models and the actual/observable plane waves. While approximations yield results/plane waves that are similar to the full dispersion model when k is below the inflection point the large k behaviour differs significantly.
Superposition in periodic potential and the equilibrium no-interaction case. As an example for deriving dynamical equilibrium plane waves we consider the generalised dispersive PDE for a wave moving in the reference frame with velocity v in a periodic potential V per , i.e.
Now let us make the separation of variables ansatz for a superposition of two different k-modes as a toy model of a wave packet and to study the dispersion between the two k− modes: i b First we note the linearity of the ☆ q operator, i.e.
On the other hand we have to satisfy if the external potential is given by Thus we obtain the analytical condensate wave function Assume a > b then it implies for the monotonically increasing polariton dispersion that ω(a) > ω(b). If we assume that the interference of the nonlinearity is vanishing and there is no external potential we would observe the same behaviour. Here note that the "wave packet" consisting of two plane waves is coherent, iff Any localised structure, such as moving bright solitons in BEC, has to satisfy such a condition component-wise, which is in particular satisfied for certain wave-packets in the focusing case α 1 < 0. Here the attractive interactions cancel out the dispersive effects from the kinetic energy due to opposing sign 26 . Note that the theoretical description in the inertial frame moving with v the plane waves does not include a Doppler term iv∇ψ and thus in such a frame the condition applies. In this sense temporal coherence of a wave packet could be feasible, however note that taking into account interference terms due to the nonlinearity will modify this behaviour.
Spin sensitive results for plane and linear waves. The simplest scenario to begin an examination of the effects of the dispersion relation in spin sensitive systems is the case of a homogeneous external potential (absorbed in the chemical potential), with pumping and decay in the simplest approximate form ref. 23, where the incoherent mode equation becomes  iη). We simplify the consideration to (1 + 1) d and note that the higher dimensional case is analogous. Furthermore we use the notation ω =â  We solve the spin sensitive system under the simplifying (but not necessary) assumption Γ ± = Γ by the analytic plane wave solutions for the two polariton spin components ±,  23. We note that for slowly varying ψ ± (x) and P ± (x) we set P ± → P ± (x) in (34).
It is useful to separate the exact motion into bulk motion plus low amplitude acoustic disturbances to understand the elementary excitations 27 . So we consider the linear waves by making an ansatz of the form where ϕ ± represents the unperturbed part solving the mean-field model such as the plane wave solutions presented above and the linear waves . By inserting in the spin sensitive PDEs and dropping terms of order δψ 2 and including a chemical potential μ ± = n ± α 1 , we get the Bogoliubov equations in k-space for the linearised perturbation dynamics, i.e. Scientific RepoRts | 7: 1891 | DOI:10.1038/s41598-017-01113-8 To arrive at those equations we have used the linearizations ψ ψ φ δψ , i.e. dropping all terms of order δψ 2 and higher, |ϕ ± | 2 ϕ ± and |ϕ ± | 2 ϕ ± for the excitations dynamics, i.e. we have separated those equations from the bulk dynamics. So the bulk of the spinor condensate is guided by Assuming the pumping function , and by introducing the abbreviations Here for a given and implicitly defined ±  P we obtain the excitations δψ ± k , which in turn defines implicitly the physical pumping function P ± . Furthermore for Re(c) → 0 and Im(c) → 0 we recover the spin coherent case. We observe a significant modification of the polariton excitation formation by considering the wave packets explicitly given by (37) due to the functional variation of q k ( ) as presented in Fig. 3. Here the direction of propagation of the linear waves is downwards (v x = 0 and v y = −1 with v = (v x , v y )) implying symmetry breaking of the ring structure in the general framework thus showing severe changes in the phenomenology of the polariton linear waves as compared with cGPE-type excitations. More specifically, in Fig. 3(i) we observe a ring structure within the linear wave solution |δψ ± (k x , k y )| when a parabolic approximation is assumed. However, the more general model Fig. 3(ii) shows linear waves that break this ring symmetry, particularly by including a singular line orthogonal to the direction of motion. Finally Fig. 3(iii) represents the most accurate linear wave of the general dispersion, i.e. ω = q k k ( ) ( ) L , and again shows the breaking of the ring's symmetry by a disconnected singular line. This breaking of the ring symmetry is not observed when the linear wave is moving at lower velocities or at rest, as observed as well for the case of its parabolic approximation.
Following our presentation of analytical wave results we now turn to the numerical phenomenology due to the different kinetic energy model predictions.
Conservative localized pulse evolution. First we consider the scenario of a localised wave packet neglecting the non-equilibrium effects and assuming that the condensate evolves without an external potential. Thus the governing equation of motion is Figure 3. Linear wave functions |δψ(k x , k y )| of the lower polariton branch for the parabolic (i), the m(k) model (ii) and the general kinetic energy (iii) showing the structural differences of the predictions of the three kinetic energy models. The direction of motion is in −k y direction and the specific parameters corresponding to the figures are set as follows.
. Further we subtract ω L (0) from the dispersions for better visibility.
where V(x) = V 0 (x) + δP(x) and κ/α = 1.36 and further parameters are specified in the Supplemental Material. As an initial condition we first set  Fig. 4 we show that spatial localisation of the initial data can be preserved to a higher degree above the inflection point at µ −  k m 2 1 despite the dispersive character of the kinetic energy. For k = 0 the initial wave function disperses rapidly. Corresponding results for the parabolic kinetic energy show stronger dispersion for large k, i.e. above the inflection point as presented in Fig. 5. These observations mimic the behaviour of bright-type solitons observed in ref. 28. We note also that the velocity of the wave packet reduces at larger k > k inf , when considering the full dispersion relation, as the velocity of the wavepacket is given by the derivative of the dispersion. At k = 4 the dispersion is flatter than at k = 2.2, so the wavepacket is slower (compare the slope of the trajectory in Fig. 4(b,c)). This stands in stark contrast to the velocity due to the parabolic dispersion, which increases with k (compare Fig. 5(b,c)).
These results are in agreement with the observation by Sich et al. 28 of localisation occuring above the inflection point. We note however one significant difference in the modeling. In ref. 28 the excitonic and photonic modes were explicitly modeled, with the complex dispersion profile of the polariton condensate emerging as a consequence of the coupling between the modes. Here instead we model only the polariton mean-field, and explicitly include the dispersion relation in the dynamics, allowing us to identify directly the contribution of the dispersion on the resulting dynamics.
Dark soliton instability. Following from (12) and (14) the single polariton state equation resembling an incoherent driving scheme can be approximately written as  24 we set the potential in x-space due to this metallic contact to be  and assume a Gaussian pumping spot resembling the spatial form of the incoherent polariton ground state formation, We refer to refs 2, 12, 13 for various applications of this model, when we assume a parabolic dispersion. As shown in ref. 29 a local abrupt change of interaction strength of a condensate establishes a stable and regular dark soliton train within a conservative GP theory. Once the flow in the direction of decreasing interaction due to particle repulsions is locally crossing the speed of sound for a scalar condensate, dark solitons are formed from dispersive shock waves at the point of abrupt change in self-interactions 30 . These solitons then proceed to dissipate the local excess of energy 24 . While in polariton condensates the interaction strength α 1 can be varied by tuning the exciton/photon detuning, and there is an ongoing debate on its experimentally measured value 31 , it is straightforward to apply a tunable potential step V 0 (x, t). The mechanism for soliton generation is again the breaking of the sound-barrier in the region x < 0 in the presence of a perturbation at x = 0 24 . In the regime of soliton-train generation, the frequency ν increases with the magnitude of the potential step as the corresponding increase of mass passing the step at x = 0 allows a more frequent breaking of the local speed of sound. In Figs 6 and 7 we observe a strong modification of the mean-field density dynamics due to the non-parabolic dispersion relation, as compared with the regular dark soliton train patterns of cGP-theory reported in ref. 24. We find that at low k the kinetic energy is in a quasi-parabolic regime supporting dark soliton solutions consistent with the graphs presented in Fig. 1. We see also in Fig. 6 that the effective mass induces additional bright-type high-density waves of varying frequency that fade out for larger times in (i) while they persist with a fixed frequency in the full dispersion regime (ii) on top of the regular dark soliton train arrays -a phenomenon entirely unobserved or neglected in cGP theory. We attribute this localisation phenomenon to the flat part of dispersion relation in the lower polariton branch. Finally we note that larger self-interaction strengths, as reported in ref. 31, would lead to chaotic dark soliton trains Fig. 7(ii) while again stable patterns are reported for the parabolic dispersion approximation for appropriate parameters, thus providing an indirect test for the large interaction hypothesis.

Discussion
The topic of kinetic energy and particularly negative effective mass of polaritons has received increasing attention 7,8,28,32 . Here we find that the effects on the pattern formation can be interpreted in the way suggested so far only for small k below the inflection point, because the kinetic energy of the polariton is composed of an array of terms (which can be stated in terms of Taylor's theorem). To directly and unambiguously observe the switch in . Density plots |ψ(x, t)| 2 of the condensate wave function predictions for the kinetic energy of the effective mass model (i) and the general polariton dispersion (ii) with a stronger self-interaction strength which is 50 times larger than the one used for Fig. 6.
sign of the effective mass, i.e. the second derivative term plus offset within the condensate wave function, as e.g. reported in ref. 28, would require a situation including physical processes that neutralise all the other terms in the expansion (10). Energy shifts are obtained by external potentials, however the variation of the remainder terms at the inflection point k ~ 2 μm −1 of the lower polarition dispersion is larger than that of the effective mass and thus overlaps its sign switch -the dispersion is positive for all k. Furthermore the dispersions of the polariton branches are monotonically increasing in k and thus there is no switch in sign of the kinetic energy. This implies that e.g. the time dependent phase of plane wave solutions due to the kinetic energy cannot be cancelled out by the contribution due to nonlinear interactions. This in turn puts into question the mechanism of bright solitons reported in ref. 28, since those bright solitons are agued by the observation that negative mass and repulsive self-interactions are formally equivalent to positive mass with attractive interactions, i.e. ψ µ ψ ψ ψ , with chemical potential µ = n(0)/2. Such a model implies the well-know bright soliton solutions, because the dispersive effects are compensated by the attractive self-interaction forces. Now when considering the toy model case of a superposition of two plane waves including the polariton dispersion as kinetic energy we have found that they move coherently in time, iff when interference terms are neglected. This illustrates that due to α 1 , ω > 0 only for specific amplitudes coherent motion is feasible. Such condition could be satisfied in ref. 28. Furthermore we have simulated the movement of a gaussian wave packet at different positions of the dispersion and can indeed confirm that the dispersion of the wave packet is reduced above the point of inflection, when considering the full polariton kinetic energy -the spatial localisation of an initial wave is stronger in the more accurate full dispersion model compared to predictions of the parabolic cGP, thus supporting the observation of localisation or suppressed dispersion in ref. 28. However note that the term soliton or bright soliton is defined rigorously 33 and involves properties such as unchanged shape and speed over time and particularly after a collision with another soliton and thus one may more accurately refer to these temporarily localised structures as near-bright-type solitons. We suggest that localisation of a wave packet in x-space is due the flat dispersion relation of the lower polariton branch allowing large k modes to be occupied with less energy than in cGPE models. In addition we note that an alternative route to bright soliton generation in polariton condensates could be along the lines of effective attractive self-interactions as discussed in ref. 34 with results similar to matter-wave bright soliton formation in ultra-cold lithium-7 gases 35,36 . Apart from this the signs of the non-parabolic dispersion relation of polaritons are apparent in many aspects of wave formation starting from bright-type solitons trains on top of dark solitons trains to chaotic dark solitons to the explicit form of linear waves of the polariton spin modes. Thus for accurate description of the polariton mean field mode the full dispersion should be included in future dynamical models for more accurate predictions.

Methods
We used a generalised Gross-Pitaevskii theory to analyse the impact a general form for the kinetic energy has on possible polariton condensate dynamics. The mathematical analysis presented is based on standard analytical tools, complex algebra and integration, and the Bogoliubov approach to linear waves. An efficient numerical method was employed for the simulation of the partial differential equation. In time, the second order time splitting method, i.e. the Strang splitting, is used for the time evolution. In space, the wave function ψ(x, t) is discretized on uniform grids and the spectral method based on the Fourier series is used to deal with the generalized kinetic energy term q*ψ(x, t). For one time step, the algorithm works as follows: • solve i∂ t ψ(x, t) = q*ψ(x, t) for half a time step via the discrete Fourier-Transform (DFT).
• solve the remaining part for one time step via the 4th order Runge-Kutta method.
In the simulations, we choose Δh = 0.005 and Δt = 0.005. A large domain [−50, 50] is used to make sure the solution on the boundary is negligible.