Phase ordering kinetics of a nonequilibrium exciton-polariton condensate

TL;DR

This paper studies the universal phase ordering dynamics in a nonequilibrium exciton-polariton condensate, revealing different scaling behaviors of vortex-antivortex annihilation depending on system parameters.

Contribution

It demonstrates that the late-time coarsening process follows universal kinetics with distinct scaling laws, extending understanding of nonequilibrium condensate dynamics.

Findings

Late-time dynamics exhibit universal phase ordering.

Scaling of correlation functions follows known models.

Different regimes show distinct dynamical exponents.

Abstract

We investigate the process of coarsening via annihilation of vortex-antivortex pairs, following the quench to the condensate phase in a nonresonantly pumped polariton system. We find that the late-time dynamics is an example of universal phase ordering kinetics, characterized by scaling of correlation functions in time. Depending on the parameters of the system, the evolution of the characteristic length scale L(t) can be the same as for the two-dimensional XY model, described by a power law with the dynamical exponent z \aprox 2 and a logarithmic correction, or z \aprox 1 which agrees with previous studies of conservative superfluids.

Tables

Table 1. Table 1: Parameters used in simulations to obtain the data presented in the Figures.

Figures	1(a,b), 2(a) and 4(a)	2(b), 3(a,b) and 4(b)
A	0	0.1
$g_C$	$3.4\mu eV\mu m^2$	$3\mu eV\mu m^2$
$g_R$	$7.2\mu eV\mu m^2$	$6\mu eV\mu m^2$
${\gamma }_C^{-1}$	$50ps$	$3.3ps$
${\gamma }_R^{-1}$	$8ps$	$3ps$
R	$5.5\times {10}^{-3}\mu m^2$	$2.3\times {10}^{-4}\mu m^2$
P	$40\mu m^{-2}p{s}^{-1}$	$3\times {10}^3\mu m^{-2}p{s}^{-1}$
$\xi$	$2.9\mu m$	$1.6\mu m$

Equations

$$g^{(1)}(\mathbf{d},t)$$

$$=f(d/L(t)).$$

$$i\text{d}\psi =$$

$$+\frac{i}{2}\left(R{n}_R-{\gamma }_C\right)]\psi \text{d}t+\text{d}W,$$

$$\frac{\partial n_R}{\partial t}=$$

$$⟨\text{d}W(\mathbf{x})\text{d}{W}^{\ast }({\mathbf{x}}^{\prime })⟩$$

$$⟨\text{d}W(\mathbf{x})\text{d}W({\mathbf{x}}^{\prime })⟩$$

$$L(t)\sim {\left(\frac{t}{\ln (t/t_0)}\right)}^{1/z},$$

$$g^{(1)}(x,t)=\frac{1}{N}⟨{\psi }^{\ast }(-x,0,t)\psi (x,0,t)⟩$$

Full text

I Introduction

One of the major achievements of statistical physics is the ability to describe complex systems of many particles with a limited set of variables describing their collective behavior. The universality of phase transitions is a particularly striking example of such reduction, where the multitude of physical models is divided into a finite number of universality classes characterized by certain symmetry properties and critical scaling laws. Phase transitions from disordered to ordered states are often accompanied by the creation of defects, such as domain walls, vortices or strings KibbleZurek . In most realistic situations, these defects subsequently decay in time, and the system undergoes gradual phase ordering or coarsening Bray_PhaseOrderingKinetics . In this nonequilibrium, late-time stage of dynamics, physical systems frequently exhibit universality characterized by a single length scale $L(t)$ that dictates the temporal evolution of all relevant quantities, such as correlation functions. The knowledge of symmetries and the character of the dominant coarsening process is sufficient to determine the evolution of this length scale. The theory of universal coarsening has been successfully applied to a wide variety of systems, from metallurgy and phase separation of fluids Cahn ; Furukawa_Scaling to biological systems and opinion dynamics Castellano_StatisticalSocial .

Recently, the classical concept of phase ordering kinetics was extended to the quantum realm in studies of atomic Bose-Einstein condensates Sachdev_PhaseOrdering ; Moore_Spinor ; Sarma_Coarsening ; Svistunov . In both the spinless Sachdev_PhaseOrdering and spinor superfluid gases Moore_Spinor ; Sarma_Coarsening ; Kawaguchi_DomainGrowth ; Kawaguchi_FerromagneticCoarsening ; Blakie_CoarseningSpin1 links with the corresponding classical systems were established. In the area of quantum fluids of light Carusotto_QuantumFluids ; Polaritons ; Kasprzak_BEC , spontaneous creation of vortices during nonadiabatic exciton-polariton condensation was observed in Deveaud_VortexDynamics and investigated theoretically in the context of the Kibble-Żurek mechanism Matuszewski_UniversalityPolaritons ; Liew_InstabilityInduced . Vortex dynamics was a topic of many studies, eg. Vina_VorticesCoherently ; Bramati_AngularVortexChain ; Deveaud_HydrodynamicVortices ; Deveaud_VortexStreets ; Fraser_VortexAntivortex ; Rubo_WarpingVortices ; Skolnick_InteractionsOnVortices ; Szymanska_TOPOVortices ; Yamamoto_VortexPair ; Whittaker_TriggeredVortices ; Deveaud_SelectivePhotoexcitation ; Baumberg_GeometricallyLocked , and the process of vortex-antivortex annihilation was observed experimentally in Vina_VortexAntivortex ; Sanvitto_AllOpticalControl ; Bramati_MergingVortices . However, to date the universal coarsening dynamics has not been investigated.

Here, we verify the scaling hypothesis in the model of a nonresonantly pumped polariton condensate Wouters_ExcitationSpectrum ; Bobrovska_Adiabatic ; Szymanska_NonequilibriumBKT ; Sanvitto_TopologicalOrder . We consider two sets of parameters such as material constants, pumping power etc. We find examples of universal phase ordering with complete collapse of correlation functions after rescaling spatial coordinates by the length scale $L(t)$. The length scale evolves according to a power law with the exponent $z$ depending on the parameters. In the first case, $z\approx 2$ with a logarithmic correction, as predicted previously for two-dimensional vector or complex fields with purely diffusive dynamics Bray_PhaseOrderingKinetics ; Pargellis_DefectDynamics . In the second case, the dynamical scaling of $L(t)$ is found to be the same as determined previously for conservative superfluids Sachdev_PhaseOrdering ; Blakie_CoarseningSpin1 , with $z\approx 1$. This shows that polariton systems can display various types of universal dynamics, which can be achieved by modifying the material parameters of the sample.

It is important to note that in this work we consider the time-evolving properties of a system that has suddenly crossed a phase transition, and not the critical properties of the phase transition itself. The latter have been intensively investigated in both recent experimental Yamamoto_PNAS ; Yamamoto_AlgebraicBKT ; Sanvitto_TopologicalOrder theoretical Sieberer_DynamicalCritical ; Szymanska_NonequilibriumBKT works. It has been claimed that polariton systems display a kind of dissipative Berezinskii-Kosterlitz-Thouless (BKT) transition, while critical exponents may differ from the ones obtained in thermal equilibrium Yamamoto_PNAS ; Yamamoto_AlgebraicBKT ; Szymanska_NonequilibriumBKT ; Sanvitto_TopologicalOrder . Here, we assume that the system is sufficiently far away from the critical point on the ordered side of it, so the system converges to an approximately defect-free phase.

The idea of universal coarsening dynamics is grounded on the scaling hypothesis, which states that at late times there is a single characteristic length scale describing the large-scale features of the system Furukawa_Scaling ; Bray_PhaseOrderingKinetics . The configuration of defects remains unchanged in time, in the statistical sense, if the spatial coordinates are scaled by this length scale which usually grows according to a power law $L(t)\sim t^{1/z}$. Here $z$ is the nonequilibrium dynamical exponent, which is in general different from the dynamical critical exponent of the phase transition which may have produced the defects in the first place KibbleZurek ; Biroli_IsingAnnealing . Consider, for instance, the first order, equal-time correlation function. It follows that the following scaling holds

[TABLE]

with $d=|{\bf d}|$ and $f(0)=1$, where $N=\langle\int{|\psi|^{2}}{\rm d}\mathbf{r}\rangle$. While the scaling hypothesis has been rigorously proven only in several cases, numerical studies indicate its validity in many physical systems Bray_PhaseOrderingKinetics .

The value of the exponent $z$ depends in general on the dimensionality of the system, the character of the coarsening processes (diffusive, inertial, etc.), symmetries, and conservation laws Bray_PhaseOrderingKinetics . For non-conserved scalar fields, such as the Ising model or model A of diffusion-reaction Hohenberg_DynamicalCriticalPhenomena , it takes the value $z=2$. In the case of conserved scalar fields the coarsening is slower with $z=3$, which can be understood as the effect of the reduced number of accessible intermediate states Moore_Spinor . When the transport is inertial rather than diffusive, faster scaling with $z=3/2$ is predicted Furukawa_Droplet ; Kawaguchi_DomainGrowth .

In the case of vector or complex fields, the existence of topological defects often dominates the phase ordering dynamics. This leads to different values of $z$, and in some cases, to logarithmic corrections, with the notable case of the two-dimensional XY model displaying the $L(t)\sim(t/\ln t)^{1/2}$ dependence Pargellis_DefectDynamics ; Yurke_CoarseningXY ; Bray_XYModel ; BrayBlundell_OnModel ; Jelic_QuenchXYModel . For the subclass of diffusive models, arguments based on the comparison of the local and global energy change allow for the prediction of the growth laws in the general case Bray_EnergyScalingApproach ; Bray_GrowthLaws . Renormalization group methods are also useful, however they fail to predict the logarithmic corrections Bray_PhaseOrderingKinetics ; Bray_EnergyScalingApproach , which may turn out to be significant Bray_XYModel .

II Model

We consider the wave function of nonresonantly pumped polariton condensate $\psi({\bf x},t)$ coupled to the reservoir described by a density field $n_{\rm R}({\bf x},t)$  Wouters_ExcitationSpectrum ; Wouters_ClassicalFields

[TABLE]

where $P$ is the exciton creation rate determined by the external optical or electrical pumping, $m^{*}$ is the effective mass of lower polaritons, $\gamma_{\rm C}$ and $\gamma_{\rm R}$ are the polariton and reservoir loss rates, and $R$ is the rate of stimulated scattering from the reservoir to the condensate, and $g_{\rm C}$, $g_{\rm R}$ are the rates of repulsive polariton-polariton and reservoir-polariton interactions, respectively. We also introduced $D=1-iA$ with $A$ being a small constant accounting for the energy relaxation in the condensate Bobrovska_Stability ; Sieberer_DynamicalCritical ; Bloch_ExtendedCondensates ; Bloch_GapStates . Alternatively, one may introduce a complex coefficient in front of the time derivative term Wouters_Superfluidity ; Wouters_EnergyRelaxation . The above phenomenological model has been successful in describing a number of different experimental situations in exciton-polariton condensates Deveaud_VortexDynamics ; Yamamoto_VortexPair . The complex stochastic noise $\textrm{d}W$, corresponding to disturbance associated with particles incoming and leaving the condensate, can be obtained within the truncated Wigner approximation Wouters_ClassicalFields

[TABLE]

In the absence of noise, a spatially uniform solution is given by $\psi({\bf x},t)=\psi_{0}e^{-i\mu_{0}t}$ and $n_{R}({\bf x},t)=n_{R}^{0}$. Above the threshold pumping $P>P_{\rm th}=\gamma_{C}\gamma_{R}/R$ a stable condensate exists with the condensate density $|\psi_{0}|^{2}=(P/\gamma_{C})-(\gamma_{R}/R)$ and $\mu_{0}=g_{C}|\psi_{0}|^{2}+g_{R}n_{R}^{0}$. We define the healing length, which corresponds to the size of the vortex core as $\xi=2\pi\hbar/\sqrt{mg_{\rm C}|\psi_{0}|^{2}}$. The scaling hypothesis may be valid only when typical distance between vortices is much larger than their size, i.e. the condition $L(t)\gg\xi$ must be fulfilled.

III Results

To investigate the process of phase ordering, we solved Eqs. (2) numerically on a rectangular mesh with size $l=150\,\mu$m with periodic boundary conditions. Starting form an empty initial condition, $\psi,n_{\rm R}=0$, the emergence of a polariton condensate from Wigner noise for $P>P_{\rm th}$ is accompanied by the spontaneous creation of phase vortices in the process analogous to the Kibble-Żurek mechanism. The detailed description of this process was presented in Matuszewski_UniversalityPolaritons ; Liew_InstabilityInduced . Here, however, we are not interested in the process of defect creation, but rather in the long-time dynamics of coarsening, which occurs when the defects are already established. We find that, analogously as in the case of the two-dimensional XY model Pargellis_DefectDynamics ; Yurke_CoarseningXY ; Bray_XYModel ; BrayBlundell_OnModel ; Jelic_QuenchXYModel ; Kawaguchi_FerromagneticCoarsening , it occurs predominantly via annihilation of vortex-antivortex pairs.

To illustrate this, in Fig. 1 we present snapshots of the amplitude and phase of the condensate wave function. The number of vortex-antivortex pairs decreases monotonously. In Figure 2, we display the evolution of the condensate density and the number of vortices in function of time. The number of vortices is estimated from the number of points on the numerical mesh where the wave function is approximately zero and phase winding occurs, at a specific time. Early in the evolution this number is very large as we start from a disordered state at low density, but actual vortices become well established only when the density becomes large.

The results are displayed for two sets of parameters, corresponding to two situations that can occur in the system. In Fig. 2(a), the pair annihilation is effective already at the stage of the dynamics when the condensate density is not yet fully established due to slow saturation of the density. This means that the defect creation and phase ordering overlap temporally, and a clear distinction between the two processes is not possible, similar as in Zurek_WhenSymmetryBreaks2D ; Biroli_IsingAnnealing ; Jelic_QuenchXYModel ; Deutschlander_KZMColloids . On the other hand, if parameters of the system are chosen such that they correspond to a lower quality sample, with a shorter polariton lifetime, the stationary density is established more quickly, and the dynamics for $t\gtrsim 20\,$ps is practically purely due to phase ordering, as shown in Fig. 2(b). The parameters used throughout the paper are summarized in Table 1.

A more complete information about the statistical properties of the system is given by the correlation functions. We confirm the scaling hypothesis by directly verifying the scaling property of the first order correlation function (1). In Figure 3(a) we show $g^{(1)}(d)$ plotted at several instants of time during the pure phase ordering stage, averaged over 16 realizations of the truncated Wigner simulations. As an estimate the length scale $L(t)$ in (1) we choose the value of $d$ for which the condition $g^{(1)}(d)=0.25$ is fulfilled. We note that in contrast to the case of scalar fields, where the correlation function often exhibits an oscillatory tail, in the present case there are no oscillations, which is generally the case if sharp domains walls are absent Bray_XYModel ; BrayBlundell_OnModel ; Blakie_CoarseningSpin1 . We obtain a perfect collapse for the scaled correlation function $f(d/L(t))$, which confirms that the scaling hypothesis is valid in this case, see Fig. 1(b).

The verification of the scaling hypothesis allows one to expect a particular form of the scaling law for the time dependent length scale $L(t)$. We find that in the “clean” case of Fig. 2(b) it follows closely the scaling law predicted for two-dimensional systems with a vector order parameter

[TABLE]

with $z\approx 2$, as shown in Fig. 4(b). In particular, it is the same as in the case of the XY model in two dimensions Pargellis_DefectDynamics ; Yurke_CoarseningXY ; Bray_XYModel ; BrayBlundell_OnModel ; Jelic_QuenchXYModel . We note that the logarithmic correction stems from the existence of a second relevant length scale and is absent if the initial conditions contain no free vortices Bray_XYModel . This scaling law is different from the one predicted for conservative atomic condensates both in the spinless Sachdev_PhaseOrdering and spinor cases Blakie_CoarseningSpin1 , where $z\approx 1$ was obtained. Nevertheless, an atomic condensate model including the effects of dissipation Kawaguchi_FerromagneticCoarsening predicted $z=2$. This highlights the crucial difference between the conservative and dissipative systems from the point of view of coarsening and demonstrates that in the case of polariton condensates dissipation is essential.

The above scaling law with $z=2$ can be explained by the balance between the vortex-antivortex attractive force and the effective friction Bray_EnergyScalingApproach ; Kawaguchi_FerromagneticCoarsening ; Bray_XYModel . Consider an isolated (anti)vortex of the form $\psi=A(r,t){\rm e}^{\pm i\phi(r,t)-i\mu_{0}t}$. When the dynamics is diffusive Wouters_ExcitationSpectrum ; Szymanska_NonequilibriumCondensation , far from the vortex core $A\approx|\psi_{0}|$ and the evolution of phase is given by the Kardar-Parisi-Zhang equation Diehl_TwoDimensionalSuperfluidity ; Bobrovska_Adiabatic . In the $L(t)\rightarrow+\infty$ limit this equation reduces to $\partial\phi/\partial t\approx-(1/\Gamma)\delta H/\delta\phi$, where $H$ is the nonlinear sigma model Hamiltonian corresponding to the kinetic part of (2). The energy of the vortex is divergent as $E_{v}\sim\ln(l/a)$, where $l$ is the system size and $a\approx\xi$ is the “microscopic” cutoff Fetter_vortices ; Bray_XYModel . For a vortex-antivortex configuration, $l$ is replaced by $R$, the distance between the vortices. From the pair energy we obtain the attractive force $F=-dE/dR\sim 1/R$. The energy dissipation for a vortex moving with velocity $v$ can be calculated as $dE/dt=\int d^{2}r(\delta H/\delta\phi)(\partial\phi/\partial t)\sim-\int d^{2}r(\partial\phi/\partial t)^{2}=-v^{2}\int d^{2}r(\partial\phi/\partial x)^{2}\sim-v^{2}E_{v}$, with the friction constant $\gamma\sim E_{v}\sim\ln(R/a)$. The evolution of the average distance between pairs is $dR/dt\sim F/\gamma\sim 1/R\ln(R/a)$, which results in $R(t)\sim(t/\ln(t/t_{0}))^{1/2}$.

We also investigate in more detail the case of the second set of parameters from Fig. 2(a). As we mentioned before, here the phase ordering dynamics is not pure, but takes place when the stationary state density is yet to be established. We find that although in this case the scaling hypothesis does not precisely describe the system dynamics, some quantitative predictions can still be formulated about the phase ordering process. Indeed, we find only slight deviations from the collapse of the $g^{(1)}$ correlation functions (not shown), which is related to the existence of a second time scale corresponding to the slow saturation of density and not to the coarsening. In Fig. 4(a) we show the time evolution of the length scale in this case. The data fits well to the theoretical prediction (4) with $z\approx 1$, in agreement with the results on vortex driven coarsening in ferromagnetic atomic condensates in the easy-axis configuration Blakie_CoarseningSpin1 . This result is further supported by the value of exponent $z=1$ for model E of superfluid helium Hohenberg_DynamicalCriticalPhenomena , and the estimated value $z\approx 1.1$ for conservative spinless condensate Sachdev_PhaseOrdering . We note that in this regime the ratio of the average distance between the vortices to the vortex size is larger, $L(t)/\xi=17.2$ (at $t=150$ ps) as compared to the previous case where $L(t)/\xi=4.5$ (at $t=90$ ps). The possible transition between the two scalings in function of this parameter would require more detailed numerical study with longer evolution times and computational box sizes.

We note that at large distances, the force between vortex and antivortex may become repulsive Diehl_Duality , which could lead to slowing down of the annihilation, and the saturation of the length scale L(t) at late times. However, we did not observe such behaviour for the parameters that were considered.

In conclusion, we confirmed that universal phase ordering can occur in exciton-polariton condensates. We found a scaling regime corresponding to purely diffusive dynamics and the one which is similar as in systems with conservative dynamics, while in the latter case the precise physical interpretation is not clear. We note that this is not the only system to display various universal behavior in different parameter regimes; in binary liquids diffusive, viscous hydrodynamic and inertial hydrodynamic regimes exist with different values of critical exponents Furukawa_Droplet ; Bray_PhaseOrderingKinetics . To our best knowledge polariton condensates are unique in that the transition is between the scaling laws determined by the dynamics of topological defects.

We thank Alejandro Zamora, Marzena Szymańska, and Nikolaos Proukakis for stimulating and valuable discussions. We acknowledge support from National Science Center, Poland grants DEC-2011/01/D/ST3/00482 and 2015/17/B/ST3/02273.

Appendix A Phase ordering with hard-wall boundary conditions

In Fig. 5 we show the evolution of the length scale in the case with Dirichlet boundary conditions for the condensate wave function $\psi(\mathbf{r},t)$. The obtained fit with dynamical exponent $1/z=1.12\pm 0.10$ agrees with the one obtained in the case of periodic boundary conditions, see Fig. 4 (b) in the main text. Note that error bars are larger than in the periodic boundary case solely due to the way we perform averaging of the correlation function in this case. As the condensate density always tends to zero close to the Dirichlet boundaries, we are no longer able to average over spatial coordinates r as in Eq. (1) in the main text. Instead we calculate correlation function from $x$ and $-x$ points on the sample

[TABLE]

The absence of averaging over spatial coordinates leads to larger variations in the estimated correlation functions and the length scale L(t).

Appendix B Decay of the number of vortices

We used two automated methods of vortex counting, one based on the counting of local density dips, and the other based on calculation of local phase winding around a particular point on the grid. After the initial stage of evolution, when the condensate density is relatively high, the vortices are well defined and the two methods give the same results. The evolution of the vortex number is presented in Fig. 6. It is in excellent agreement with the scaling of correlation function, taking into account that the number of vortices scales as $N\sim t^{-d/z}$ with $d=2$ being the number of dimensions and $z=2$.

Appendix C The effect of the logarithmic correction

In Figure 7 we show the result of fitting the numerical data to the pure algebraic function without the logarithmic correction. The fit is also very good, but the value of the scaling exponents do not agree with theoretical predictions, Eq. (4). This is very similar to the situation described in the XY model Bray_XYModel , where this issue was discussed at length. The logarithmic correction effectively changes the slope of the fit, but does not result in significant bending in the bi-logarithmic scale. Such bending could be observable at early times t, but in this limit the universal scaling is not valid.