Spin–orbit coupling and nonlinear modes of the polariton condensate in a harmonic trap

The model that we consider in this work consists of coupled GPEs for the two-component polariton spinor wave function $({{\rm{\Psi }}}_{+},{{\rm{\Psi }}}_{-})$ [20, 24]. We assume that the polariton condensate is pumped by means of a nonresonant optical scheme, which creates a density of coherent and incoherent excitons [29, 38]. Adiabatic elimination of the latter yields the following system of coupled equations [38]:

$$\begin{eqnarray}{\rm{i}}{\partial }_{t}{{\rm{\Psi }}}_{+} & = & -\displaystyle \frac{1}{2}(1-{\rm{i}}\eta )({\partial }_{x}^{2}+{\partial }_{y}^{2}){{\rm{\Psi }}}_{+}+(| {{\rm{\Psi }}}_{+}{| }^{2}+\alpha | {{\rm{\Psi }}}_{-}{| }^{2}){{\rm{\Psi }}}_{+}\\ & & +\beta {({\partial }_{x}-{\rm{i}}{\partial }_{y})}^{2}{{\rm{\Psi }}}_{-}+{\rm{i}}(\varepsilon -\sigma | {{\rm{\Psi }}}_{+}{| }^{2}){{\rm{\Psi }}}_{+}+\left[\displaystyle \frac{1}{2}({x}^{2}+{y}^{2})+{\rm{\Omega }}\right]{{\rm{\Psi }}}_{+},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}{\rm{i}}{\partial }_{t}{{\rm{\Psi }}}_{-} & = & -\displaystyle \frac{1}{2}(1-{\rm{i}}\eta )({\partial }_{x}^{2}+{\partial }_{y}^{2}){{\rm{\Psi }}}_{-}+(| {{\rm{\Psi }}}_{-}{| }^{2}+\alpha | {{\rm{\Psi }}}_{+}{| }^{2}){{\rm{\Psi }}}_{-}\\ & & +\beta {({\partial }_{x}+{\rm{i}}{\partial }_{y})}^{2}{{\rm{\Psi }}}_{+}+{\rm{i}}(\varepsilon -\sigma | {{\rm{\Psi }}}_{-}{| }^{2}){{\rm{\Psi }}}_{-}+\left[\displaystyle \frac{1}{2}({x}^{2}+{y}^{2})-{\rm{\Omega }}\right]{{\rm{\Psi }}}_{-}.\end{eqnarray} \tag{ 2 }$$

Here, gain $\varepsilon \gt 0$ is determined by the rate of the pumping of the polariton density, $\sigma \gt 0$ is the gain-saturation coefficient that prevents an unphysical blow-up and allows for the time-independent finite density solutions to exist in the system, and η is an effective polariton diffusion coefficient (viscosity), linked to the carrier diffusion [38, 39]. Further, the strength of the repulsive intra-component interactions is normalized to be 1, while $\alpha =-0.05$ accounts for the inter-component attraction, and β is the strength of the polariton SOC [20], while the strength of the harmonic-oscillator potential is scaled to be 1. The potential corresponds to the shift of the cavity detuning, and originates from the trapping of the photonic component of exciton-polariton condensate. The potential can be created, for example, if either the above-mentioned open cavity with a parabolic top mirror is used, or one of the many available patterning techniques is applied to the top mirror in the traditional closed-resonator setting with Bragg reflectors [23]. Lastly, Ω characterizes the Zeeman-splitting energy, which is proportional to the applied magnetic field. Accessible physical values of the parameters in the linear Hamiltonian part of equations (1), (2), and, in particular, the SOC strength in the presence of the harmonic-oscillator potential have been discussed in detail in [28]. That work demonstrates that SOC energy may reach 1 meV, which is comparable to the characteristic harmonic-oscillator energy, thereby making values of β between 0 up to 1 physically relevant. Achieving large polarization energy splitting has been important, since the beginning of the studies of microcavity polaritons, in the context of the polariton spin Hall effect [33, 34], see also [4] for a detailed review of the previous work. Currently, the design of strong SOC is getting crucially important for the work towards experimental observation and applications of polariton topological insulators [5, 35–37].

Stationary solutions to equations (1) and (2), with chemical potential μ, are sought in the form of

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\pm }={{\rm{e}}}^{-{\rm{i}}\mu t}{\psi }_{\pm }(r,\theta ),\end{eqnarray} \tag{ 3 }$$

where ${\psi }_{\pm }$ are complex functions of polar coordinates $(r,\theta )$. The substitution of expressions (3) in equations (1) and (2) leads to the following stationary equations:

$$\begin{eqnarray}\mu {\psi }_{+} & = & -\displaystyle \frac{1}{2}(1-{\rm{i}}\eta )\left(\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {\theta }^{2}}\right){\psi }_{+}+(| {\psi }_{+}{| }^{2}+\alpha | {\psi }_{-}{| }^{2}){\psi }_{+}\\ & & +\beta {\left[{{\rm{e}}}^{-{\rm{i}}\theta }\left(\displaystyle \frac{\partial }{\partial r}-\displaystyle \frac{{\rm{i}}}{r}\displaystyle \frac{\partial }{\partial \theta }\right)\right]}^{2}{\psi }_{-}+{\rm{i}}(\varepsilon -\sigma | {\psi }_{+}{| }^{2}){\psi }_{+}+\left(\displaystyle \frac{1}{2}{r}^{2}+{\rm{\Omega }}\right){\psi }_{+},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\mu {\psi }_{-} & = & -\displaystyle \frac{1}{2}(1-{\rm{i}}\eta )\left(\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {\theta }^{2}}\right){\psi }_{-}+(| {\psi }_{-}{| }^{2}+\alpha | {\psi }_{+}{| }^{2}){\psi }_{-}\\ & & +\beta {\left[{{\rm{e}}}^{{\rm{i}}\theta }\left(\displaystyle \frac{\partial }{\partial r}+\displaystyle \frac{{\rm{i}}}{r}\displaystyle \frac{\partial }{\partial \theta }\right)\right]}^{2}{\psi }_{+}+{\rm{i}}(\varepsilon -\sigma | {\psi }_{-}{| }^{2}){\psi }_{-}+\left(\displaystyle \frac{1}{2}{r}^{2}+{\rm{\Omega }}\right){\psi }_{-},\end{eqnarray} \tag{ 5 }$$

which are used below as the basic system.

Before tackling the full nonlinear system, it is relevant to address the simplest version of stationary equations (4) and (5), which neglects the nonlinearity, gain, and diffusion, as well as assumes the symmetry between the components, i.e., ${\rm{\Omega }}=0$ (no Zeeman splitting):

$$\begin{eqnarray}\mu {\psi }_{+} & = & -\displaystyle \frac{1}{2}\left(\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {\theta }^{2}}\right){\psi }_{+}+\beta {\left[{{\rm{e}}}^{-{\rm{i}}\theta }\left(\displaystyle \frac{\partial }{\partial r}-\displaystyle \frac{{\rm{i}}}{r}\displaystyle \frac{\partial }{\partial \theta }\right)\right]}^{2}{\psi }_{-}+\displaystyle \frac{1}{2}{r}^{2}{\psi }_{+},\\ \mu {\psi }_{-} & = & -\displaystyle \frac{1}{2}\left(\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {\theta }^{2}}\right){\psi }_{-}+\beta {\left[{{\rm{e}}}^{{\rm{i}}\theta }\left(\displaystyle \frac{\partial }{\partial r}+\displaystyle \frac{{\rm{i}}}{r}\displaystyle \frac{\partial }{\partial \theta }\right)\right]}^{2}{\psi }_{+}+\displaystyle \frac{1}{2}{r}^{2}{\psi }_{-}.\end{eqnarray} \tag{ 6 }$$

It is easy to see that equation (6) gives rise to two exact eigenmodes of the linearized dissipation-free system, corresponding to independent signs ± in equation (7):

$$\begin{eqnarray}&&{\psi }_{+}={A}_{\max }r\exp \left(-{\rm{i}}\theta -\displaystyle \frac{{r}^{2}}{2\sqrt{1\mp 2\beta }}\right),\,{\psi }_{-}=\pm {A}_{\max }r\exp \left({\rm{i}}\theta -\displaystyle \frac{{r}^{2}}{2\sqrt{1\mp 2\beta }}\right),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\mu }_{\pm }^{(0)}=2\sqrt{1\mp 2\beta },\end{eqnarray} \tag{ 8 }$$

where ${A}_{\max }$ is an arbitrary constant, and superscript $(0)$ attached to ${\mu }_{\pm }$ refers to the linear approximation. These exact eigenmodes, which, as a matter of fact, are found following the pattern of standard solutions for the 2D isotropic harmonic oscillator in quantum mechanics, may be naturally called vortex–antivortex states. Note that the symmetric vortex–antivortex state with the lower eigenvalue of the energy, ${\mu }_{+}^{(0)}=2\sqrt{1-2\beta }$, exists if the SOC is not too strong, viz, at $\beta \leqslant 0.5$, while its antisymmetric counterpart with higher energy, ${\mu }_{-}^{(0)}=2\sqrt{1\,+\,2\beta }$, exists at all values of β. In fact, the disappearance of the eigenmode (7) corresponding to the top sign in $\mp$ at $\beta =0.5$, where the radial size of the eigenmode shrinks to zero as ${(1-2\beta )}^{1/4}$, may be considered as a quantum phase transition in the polariton condensate, see [40].

Both types of the vortex–antivortex states, antisymmetric and symmetric ones, were experimentally observed, in a quasi-linear regime, in [28], where they were named, respectively, modes of types (i) and (iii). The respective SOC parameter, estimated in [28], is $\beta \approx 0.04$ in our scaling.

In a sense, the exact vortex–antivortex modes are similar to the above-mentioned semi-vortices found in the nonlinear dissipation-free 2D SOC model of the atomic condensate considered in [7], which feature zero vorticity in one component and vorticity ±1 in the other. Indeed, the difference of the vorticities in the components of the semi-vortex is ${\rm{\Delta }}S=1$ because (as mentioned above too) SOC in the atomic condensates is represented by the linear differential operator of the first order [3], while in the present system the vorticity difference between the components of the vortex–antivortex is ${\rm{\Delta }}S=2$, as the SOC is represented by the second-order operator in equations (1) and (2).

If the nondissipative nonlinear terms are kept in equations (4) and (5) as small perturbations, it is easy to find the first nonlinear correction to the chemical potential, ${\mu }_{\pm }^{(1)}$, using the commonly known method in quantum mechanics, that generates the first correction as the spatial average of the perturbation potential [18, 41] . A simple calculation, which makes use of the unperturbed wave functions (7), yields

$$\begin{eqnarray}&&{\mu }_{\pm }^{(1)}=\displaystyle \frac{1+\alpha }{4}\sqrt{1\mp 2\beta }{A}_{\max }^{2},\end{eqnarray} \tag{ 9 }$$

where ${A}_{\max }^{2}$ should be small enough, to keep correction (9) small in comparison with the result (8), obtained in the linear limit.

Furthermore, if the linear gain, diffusion, and nonlinear loss are also kept as small perturbations in equations (4) and (5), then amplitude ${A}_{\max }$ of the linear mode (7), which is treated in equations (7) and (9) as a free parameter, is uniquely selected by the obvious balance condition for the total norm (it is defined below by equation (14)), written as

$$\begin{eqnarray}&&\varepsilon {\int }_{0}^{\infty }{| {\psi }_{\pm }(r)| }^{2}{r}{\rm{d}}{r}=\sigma {\int }_{0}^{\infty }{| {\psi }_{\pm }(r)| }^{4}{r}{\rm{d}}{r}+\displaystyle \frac{\eta }{2}{\int }_{0}^{\infty }\left[{\left|\displaystyle \frac{{\rm{d}}{\psi }_{\pm }}{{\rm{d}}{r}}\right|}^{2}+\displaystyle \frac{{| {\psi }_{\pm }(r)| }^{2}}{{r}^{2}}\right]{r}{\rm{d}}{r},\end{eqnarray} \tag{ 10 }$$

for the mode $| {\psi }_{+}(r)| =| {\psi }_{-}(r)|$ taken as per equation (7). The result of the simple calculation is

$$\begin{eqnarray}&&{{{({A}_{\max })}_{\mathrm{balance}}^{2}}^{}}_{}^{}=\displaystyle \frac{\varepsilon \sqrt{1\mp 2\beta }}{2\sigma (1\mp 2\beta )+\eta },\end{eqnarray} \tag{ 11 }$$

which makes sense if the pump rate ε is small enough. A more sophisticated analytical approximation for the vortex–antivortex in the full system is developed below, see equations (16)–(18).

In addition to the two species of vortex–antivortex complexes, built as per equations (7), the linearized coupled GPEs (6) give rise to another type of eigenstates, which, following [7], may be called mixed modes, as they contain terms with zero vorticity and vorticities ±2 in each component. Like the symmetric vortex–antivortex states (7), the mixed modes obey constraint ${\psi }_{+}={\psi }_{-}^{* }$, with $\ast$ standing for the complex conjugate. Linear mixed modes cannot be represented by an exact solution (unlike equation (7) for the vortex–antivortex complexes), but they can be easily constructed in an approximate form for $\beta \ll 1$, starting from an obvious ground state of the harmonic-oscillator in each component at $\beta =0$:

$$\begin{eqnarray}&&{\psi }_{\pm }\approx {A}_{\max }\exp \left(-\displaystyle \frac{{r}^{2}}{2}\right)[1\pm {\rm{i}}\beta {r}^{2}\sin (2\theta )],\ \mu \approx 1-\beta .\end{eqnarray} \tag{ 12 }$$

At larger β, the mixed mode eigenmodes were constructed numerically by means of the imaginary-time method, applied to the time-dependent version of linearized equations (6), starting with the ground-state input, ${{\rm{\Psi }}}_{+}={{\rm{\Psi }}}_{-}=\exp (-{r}^{2}/2)$. The numerically found mixed-mode shapes are displayed in figure 1, and the comparison with the analytical approximation (12), for a relevant small value $\beta =0.1$, is presented in figure 1(b).

Zoom In Zoom Out Reset image size

Further, eigenvalues corresponding to the mixed modes were compared to the analytically found eigenvalues for the vortex–antivortex states, given by equation (8), with the aim to identify the ground state, which must have the smallest eigenvalue of energy, μ, for given β. As a result, it was found that, at $\beta \lt 0.46$, the mixed modes realize the ground state (for small β, this is evident from the comparison of equations (8) and (12)), but in a narrow interval of $0.46\lt \beta \lt 0.5$ , the symmetric vortex–antivortex complex, with the upper sign in equation (7), produces a slightly smaller value of μ, which is given by equation (8), as shown in figure 1(c); recall that vortex–antivortex state with the upper sign does not exist at $\mu \gt 0.5$. The same figure 1(c) shows that the mixed mode does not exists either at $\beta \gt 0.5$ (which is another essential manifestation of the above-mentioned quantum phase transition), hence the only eigenmode which survives at $\beta \gt 0.5$ is the antisymmetric vortex–antivortex complex with the lower sign in equation (7). It is worthy to mention here that, as shown below, solely vortex–antivortex modes with opposite signs of their components are found as stable states in the full system of equations (1) and (2), while vortex–antivortex bound states with identical signs of its components cannot be generated by the full system, even if they may realize the ground state of the linearized system at β close to 0.5.

The next step is to consider equations (1) and (2) in which nonlinear terms are kept, but the dissipative ones are still ignored:

$$\begin{eqnarray}{\rm{i}}{\partial }_{t}{{\rm{\Psi }}}_{+} & = & -\displaystyle \frac{1}{2}({\partial }_{x}^{2}+{\partial }_{y}^{2}){{\rm{\Psi }}}_{+}+\displaystyle \frac{1}{2}({x}^{2}+{y}^{2}){{\rm{\Psi }}}_{+}+(| {{\rm{\Psi }}}_{+}{| }^{2}+\alpha | {{\rm{\Psi }}}_{-}{| }^{2}){{\rm{\Psi }}}_{+}+\beta {({\partial }_{x}-{\rm{i}}{\partial }_{y})}^{2}{{\rm{\Psi }}}_{-},\\ {\rm{i}}{\partial }_{t}{{\rm{\Psi }}}_{-} & = & -\displaystyle \frac{1}{2}({\partial }_{x}^{2}+{\partial }_{y}^{2}){{\rm{\Psi }}}_{-}+\displaystyle \frac{1}{2}({x}^{2}+{y}^{2}){{\rm{\Psi }}}_{-}+(| {{\rm{\Psi }}}_{-}{| }^{2}+\alpha | {{\rm{\Psi }}}_{+}{| }^{2}){{\rm{\Psi }}}_{-}+\beta {({\partial }_{x}+{\rm{i}}{\partial }_{y})}^{2}{{\rm{\Psi }}}_{+}.\end{eqnarray} \tag{ 13 }$$

Solving these equations by means of the imaginary-time method makes it possible to produce stationary solutions with a fixed norm,

$$\begin{eqnarray}&&N\equiv \int \int [{| {{\rm{\Psi }}}_{+}(x,y)| }^{2}+{| {{\rm{\Psi }}}_{-}(x,y)| }^{2}]{\rm{d}}{x}{\rm{d}}{y}=1,\end{eqnarray} \tag{ 14 }$$

and eventually identify nonlinear ground states, which realize the minimum of the Hamiltonian,

$$\begin{eqnarray}H & = & \displaystyle \int \displaystyle \int \left\{\displaystyle \frac{1}{2}(| {\rm{\nabla }}{{\rm{\Psi }}}_{+}{| }^{2}+| {\rm{\nabla }}{{\rm{\Psi }}}_{-}{| }^{2})+\displaystyle \frac{1}{2}({x}^{2}+{y}^{2})({| {{\rm{\Psi }}}_{+}| }^{2}+{| {{\rm{\Psi }}}_{-}| }^{2})+\alpha | {{\rm{\Psi }}}_{+}{| }^{2}{| {{\rm{\Psi }}}_{-}| }^{2}\right.\\ & & +\left.\displaystyle \frac{1}{2}(| {{\rm{\Psi }}}_{+}{| }^{4}+{| {{\rm{\Psi }}}_{-}| }^{4})-\beta [{(({\partial }_{x}+{\rm{i}}{\partial }_{y}){{\rm{\Psi }}}_{+})}^{* }({\partial }_{x}-{\rm{i}}{\partial }_{y}){{\rm{\Psi }}}_{-}+{\rm{c}}.{\rm{c}}.]\right\}{\rm{d}}{x}{\rm{d}}{y}\end{eqnarray} \tag{ 15 }$$

(here c.c. stands for the complex-conjugate expression), for the given norm. Figure 2(a) shows the so found profiles of $| {{\rm{\Psi }}}_{+}|$ and $| {{\rm{\Psi }}}_{-}|$ for the ground state at $\beta =0.48,0.45,0.42$. In agreement with the above-mentioned transition of the ground state from the mixed mode to the symmetric vortex–antivortex in the linearized system, as the SOC strength changes from $\beta \lt 0.46$ to $\beta \gt 0.46$, figure 2 demonstrates that, in the dissipation-free nonlinear system, the same transition takes place as the SOC strength increases from $\beta =0.45$ to $\beta =0.48$.

Zoom In Zoom Out Reset image size

The systematic numerical analysis demonstrates that the full system of equations (1) and (2) gives rise to stable stationary states only if the diffusion term is present, $\eta \gt 0$ (at $\eta =0$, chaotic solutions appear, which are not shown here in detail). For β not two large, the system with small values of $\eta \gt 0$ generates stable patterns of the mixed-mode type, which thus plays the role of the fundamental state in the full system, see an example in figure 3(a).

Zoom In Zoom Out Reset image size

At larger values of the SOC strength, β, and $\eta \gt 0$, the full system readily produces the fundamental state in the form of stable vortex–antivortex complexes with the opposite sign of the components, which, in the linear limit, correspond to the lower sign (antisymmetric bound state) in equation (7). A typical example is presented in figure 4, where panel (a) shows real parts of ${{\rm{\Psi }}}_{+}$ and ${{\rm{\Psi }}}_{-}$ (with the opposite signs, as said here), and panel (b) shows $| {{\rm{\Psi }}}_{\pm }(x)|$ (the solid line) in the cross section of y = 0. We stress that the numerical solution of the full system, including the nonlinear and dissipative terms, has not produced any stable vortex–antivortex complex with identical signs of its two component, in spite of the fact that such complexes play the role of the ground state in the dissipation-free system close to $\beta =0.5$, as shown above.

Zoom In Zoom Out Reset image size

The vortex–antivortex states produced by the full system can be quite efficiently approximated in an analytical form. To this end, the following ansatz is adopted for their two components:

$$\begin{eqnarray}&&{\psi }_{\pm }=\pm {A}_{\max }r\exp [-(1/2)\rho {r}^{2}\mp {\rm{i}}\theta ],\end{eqnarray} \tag{ 16 }$$

with free parameters ${A}_{\max }$ and ρ, see equation (7). First, ignoring the dissipative terms, it is possible to predict values of ${A}_{\max }$ and ρ by minimizing the system's energy (15) at a fixed value of norm N (this time, we do not set N = 1). The results is

$$\begin{eqnarray}&&\rho ={[1+2\beta +(1+\alpha )N/8]}^{-1/2},\,{A}_{\max }=\sqrt{N}\rho .\end{eqnarray} \tag{ 17 }$$

Further, restoring the dissipative terms ∼$\varepsilon ,\eta$, and σ, the norm, which appears as a free parameter in equation (17), can be predicted by the balance equation (10). The final result is

$$\begin{eqnarray}&&N=\displaystyle \frac{4}{{\sigma }^{2}}\{-(\sigma \eta -{\epsilon }^{2}(1+\alpha )/4)+\sqrt{{[{\epsilon }^{2}(1+\alpha )/4-\sigma \eta ]}^{2}-{\sigma }^{2}[{\eta }^{2}-{\epsilon }^{2}(1+2\beta )]}\}.\end{eqnarray} \tag{ 18 }$$

The comparison of the vortex–antivortex modes predicted by equations (16)–(18) with their numerically found counterparts is presented in figures 4(b) and (c), showing good accuracy of the analytical approximation. Note that figure 4(c) demonstrates that this family of the vortex–antivortex states does not extend to $\beta \lt 0.4$.

Results of the systematic numerical analysis are summarized in figure 5(a), which displays the existence areas (between the solid and dashed lines) for stable mixed-mode and vortex–antivortex states, in the plane of the most essential control parameters, viz, the SOC strength, β, and diffusion coefficient, η . As mentioned above, the mixed modes and antisymmetric vortex–antivortex complexes tend to be stable, as ground states, at smaller and larger values of $\beta ,$ respectively, with a small bistability region revealed by the figure. Above the dashed lines, the solutions decay to zero, while below the solid lines, the evolution makes them chaotic.

Zoom In Zoom Out Reset image size