Externally Controlled Lotka-Volterra Dynamics in a Linearly Polarized Polariton Fluid

TL;DR

This paper derives and analyzes a simplified Lotka-Volterra model to understand externally controlled pattern switching in a polarized polariton fluid, linking nonlinear optical dynamics with population competition theory.

Contribution

It introduces a generalized Lotka-Volterra model with external control for polariton systems, providing analytical insights into pattern switching mechanisms.

Findings

The model captures key switching behaviors observed in simulations.

Analysis of stability and bifurcations in the parameter space.

External control enables manipulation of pattern states.

Abstract

Spontaneous formation of transverse patterns is ubiquitous in nonlinear dynamical systems of all kinds. An aspect of particular interest is the active control of such patterns. In nonlinear optical systems this can be used for all-optical switching with transistor-like performance, for example realized with polaritons in a planar quantum-well semiconductor microcavity. Here we focus on a specific configuration which takes advantage of the intricate polarization dependencies in the interacting optically driven polariton system. Besides detailed numerical simulations of the coupled light-field exciton dynamics, in the present paper we focus on the derivation of a simplified population competition model giving detailed insight into the underlying mechanisms from a nonlinear dynamical systems perspective. We show that such a model takes the form of a generalized Lotka-Volterra system for two competing populations explicitly including a source term that enables external control. We present a comprehensive analysis both of the existence and stability of stationary states in the parameter space spanned by spatial anisotropy and external control strength. We also construct phase boundaries in non-trivial regions and characterize emerging bifurcations. The population competition model reproduces all key features of the switching observed in full numerical simulations of the rather complex semiconductor system and at the same time is simple enough for a fully analytical understanding of the system dynamics.

Equations

$$\mathrm{i}\hslash {\partial }_t{E}^{\pm }=$$

$$\mathrm{i}\hslash {\partial }_t{p}^{\pm }=$$

$$+T^{++}|p^{\pm }{|}^2{p}^{\pm }+T^{+-}|p^{\mp }{|}^2{p}^{\pm }.$$

$${\dot{A}}_1$$

$${\dot{A}}_2$$

$${\partial }_t{\tilde{A}}_i={\tilde{A}}_i\left(r_i-\sum _{j=1}^2{c}_{ij}{\tilde{A}}_j\right)$$

$$C=2\left(\begin{array}{cc}{\beta }_1& {\theta }_1\\ {\theta }_2& {\beta }_2\end{array}\right),$$

$$S_1^{\ast }(\delta \alpha )={\beta }_2{\sqrt{\frac{1+\delta \alpha }{{\theta }_1}}}^3-\sqrt{\frac{1+\delta \alpha }{{\theta }_1}},$$

$$S_2^{\ast }(\delta \alpha )=6\sqrt{\frac{3{\left[(1+\delta \alpha ){\theta }_2-{\beta }_1\right]}^3}{{\beta }_1^2({\theta }_1{\theta }_2-{\beta }_1{\beta }_2)}}.$$

$$V=\sum _{i=1}^2\frac{{\alpha }_i}{2}{A}_i^2-\frac{{\beta }_i}{4}{A}_i^4-\frac{\theta }{4}{A}_i^2{A}_{i+1}^2+S_i{A}_i$$

$$(j,0,0),(0,j,0),(0,0,j),3\times (j,j,j),2\times (j,j+1,j+1).$$

$$\mathrm{i}\hslash {\partial }_t{p}_j^y=({ϵ}_j-\hslash \omega -\mathrm{i}{\gamma }_e)p_j^y+$$

$$+$$

$$-$$

$$+$$

$$-$$

$$\mathrm{i}\hslash {\partial }_t{E}_j^y=(\hslash {\omega }_j^y-\hslash \omega -\mathrm{i}{\gamma }_c)$$

$$E_j=\frac{\Omega }{\hslash {\omega }_j-{\delta }_j-\hslash \omega -\mathrm{i}{\gamma }_c}{p}_j\equiv \frac{2{\lambda }_j}{{\alpha }_{\mathrm{PSF}}\Omega }{\mathrm{e}}^{\mathrm{i}{\theta }_j}{p}_j$$

$${\lambda }_j{\mathrm{e}}^{\mathrm{i}{\theta }_j}=\frac{{\alpha }_{\mathrm{PSF}}{\Omega }^2}{2(\hslash {\omega }_j^y-{\delta }_j-\hslash \omega -\mathrm{i}{\gamma }_c)}\mathrm{for}j=1,2$$

$${\lambda }_0{\mathrm{e}}^{\mathrm{i}{\theta }_0}=\frac{{\alpha }_{\mathrm{PSF}}{\Omega }^2}{2(\hslash {\omega }_0^x-\hslash \omega -\mathrm{i}{\gamma }_c)}\mathrm{for}j=0.$$

$${\partial }_t{\tilde{p}}_j^y=L_j{\tilde{p}}_j^y+\sum _{k=1}^2{C}_{jk}{\tilde{p}}_k^{y2}{\tilde{p}}_j^y$$

$${\partial }_t{\phi }_j=K_j+\sum _{k=1}^2{D}_{jk}{\tilde{p}}_k^{y2}.$$

$$\hslash L_j=$$

$$-\frac{1}{2}(T^{++}-T^{+-}){\tilde{p}}_0^{x2}\sin (-2{\varphi }_j)$$

$$\hslash C_{jk}=$$

$$\hslash K_j=$$

$$-\frac{1}{2}(T^{++}-T^{+-}){\tilde{p}}_0^{x2}\cos (-2{\varphi }_j)$$

$$\hslash D_{jk}=$$

$$\partial {\tilde{p}}_j^y=\left[L_j({\varphi }_j^{(0)})+\sum _{k=1}^2\left(C_{jk}({\varphi }_j^{(0)},{\varphi }_k^{(0)})-\frac{D_{jk}({\varphi }_j^{(0)},{\varphi }_k^{(0)})L_j^{\prime }({\varphi }_j^{(0)})}{K_j^{\prime }({\varphi }_j^{(0)})}\right){\tilde{p}}_k^{y2}\right]{\tilde{p}}_j^y.$$

$${\partial }_t{\tilde{p}}_1^y$$

$${\partial }_t{\tilde{p}}_2^y$$

$${\partial }_t{A}_1$$

$${\partial }_t{A}_2$$

Full text

I Introduction

The demand for integrated optoelectronic devices in optical communication networks has resulted in an increase of research activities targeted at functional all-optical components. For example, a wide range of different approaches has been proposed to realize efficient all-optical switches exploiting the nonlinear optical properties of different material systems, including organic photonic crystals Hu2008 , rubidium atomic-vapor cells Dawes672 , or GaAs semiconductor microcavities as in Refs. ballarini2013all ; amo2010exciton ; sanvitto2016road and Refs. Kheradmand2008 ; Schumacher2009 ; PhysRevB.87.205307 . The latter utilize the optical control of transverse optical patterns to achieve transistor-like switching performance Dawes672 ; Lewandowski17 . Spontaneous formation of spatially extended patterns has been intensively studied during the past decades cross1993 with applications to different areas of science including formation of sand ripples and desert sand dunes bowmannewell98 , animal coat patterns such as zebra stripes Meinhardt1982 ; murray03 , or geographic patterns in parasitic insect populations hassel . With applications to optical switching, however, the quest for efficient external control of these patterns naturally arises, but in many cases has not been explored in detail. In the present work we investigate a specific example of all-optical switching of polariton patterns in a semiconductor quantum-well microcavity system. In contrast to our previous work Lewandowski17 we give a detailed analysis of the optical switching dynamics from a nonlinear dynamical systems perspective. To this end we derive a simplified mode competition model that governs the essentials of the system dynamics but is simple enough to fully characterize the possible stationary solutions and phase-space singularities in the relevant parameter space. This simplified model has the very general mathematical form of a generalized Lotka-Volterra system, with the addition of an inhomogeneity for external control. In the main part of the paper (cf. Sec. III), we present a complete steady state analysis of this nonlinear system. The solution space we obtain in dependence of spatial anisotropy and external control strength is of very general nature and may be similarly realized in other systems where external control of population competition is studied such as chemical reactions or in the life sciences.

As our specific example, here we study planar semiconductor microcavities with a strong coupling between the cavity field and the exciton polarization that gives rise to the formation of exciton polaritons Kavokin . These quasiparticles consist of a photonic and an excitonic part and are characterized by a normal-mode splitting in the dispersion relation, on the lower branch leading to long coherence times and strong nonlinear interactions. The latter are driven by four-wave mixing processes of coherent polariton fields which can also be interpreted as polariton-polariton scattering, mediated through exciton-exciton scattering. For certain excitation conditions, small spatially varying density fluctuations can experience huge growth in particular modes due to the intrinsic feedback mechanisms driven by four-wave mixing. This causes spatially homogeneous density distributions to become unstable such that the system’s symmetry is spontaneously broken. This results in the formation of stationary patterns, directly observable in the far field emission from the microcavity. Figure 1(a) shows the excitation geometry used with the pump at normal incidence (zero in-plane momentum) and finite off-axis ($k{\neq}0$) signals due to amplified fluctuations. Figure 1(b) shows the normal-mode splitting of the dispersion relation into a lower (LPB) and upper polariton branch (UPB) alongside the bare exciton and cavity dispersions. Phase-matching conditions determine the efficiencies of the different scattering processes, and therefore determine the possible pattern geometry. For scalar (only one circular polarization state) polariton fields hexagon patterns are favored ardizzone ; egorov2014motion ; saito2013order . Extending the model with polarization dependence a complex interplay of the TE/TM cavity-mode splitting in the linear regime and the spin-dependent exciton-exciton interactions in the nonlinear regime arises. The resulting polarization-induced spatial anisotropy PhysRevB.94.045308 determines the possible unstable modes and stable patterns formed. For linearly polarized pump excitation slightly above instability threshold, cross-linearly polarized resonant modes parallel and perpendicular to the pump’s polarization plane constitute the basic instabilities PhysRevB.77.073302 , resulting in two-spot or four-spot patterns. Making use of this spatially anisotropic polariton-polariton scattering a reversible optical switching for orthogonal two-spot patterns can be realized that is triggered by a weak external control beam Lewandowski17 . The advance of taking the polarization dependence into account, compared to scalar switching schemes Kheradmand2008 ; Schumacher2009 ; PhysRevB.87.205307 , is the explicit breaking of the rotational symmetry, offering to control the orientation of the initial pattern via a linearly polarized pump and to provide an automatic back-switching mechanism. It also makes background-free detection possible, since the pattern switching takes place in the cross-polarized channel.

II Orthogonal Switching of Two-Spot Patterns

A detailed numerical investigation of the present system for the following excitation conditions was already discussed in Ref. Lewandowski17 . The system is excited (driven) by an $x$-linearly polarized continuous wave pump at normal incidence with flat-top Gaussian shape in the QW plane and an intensity slightly above the off-axis instability threshold. In this case spontaneous breaking of spatial and polarization symmetry is observed. Signals are formed by resonant polariton scattering onto the elastic circle defined by the polariton dispersion (cf. Fig. 1). As can also be understood based on a linear stability analysis schumacher2007influence ; PhysRevB.77.073302 , the scattering occurs predominantly in four spatial directions oriented orthogonal and parallel to the pump polarization plane, respectively. In general, in any spatial direction signals can form either in the TE or in the TM mode as illustrated in Fig. 2. However, for the pump polaritons scattered to finite $k$ through Coulomb interaction, due to the underlying spin-dependent exciton-exciton interactions, the scattering probability is higher for a polarization state orthogonal to the pump polarization state schumacher2007influence . Therefore, in the spatial direction parallel to the pump polarization state, the scattered signals preferentially form in the TE mode, and in the TM mode for scattering orthogonal to the pump polarization state (cf. Fig. 2). Out of these four signals, a stationary two-spot pattern (only one mode pair with opposite in-plane momenta) can be prepared by introducing some anisotropy favouring one direction over the other. Alongside a slight polarization induced spatial anisotropy that is due to a slightly higher density of states for the TE modes PhysRevB.77.073302 , an additional anisotropy can be introduced by tilting the pump beam slightly away from normal incidence (cf. Fig. 2). In the nonlinear system studied, at sufficiently high densities of the favored two-spot pattern, cross-saturation processes will lead to the extinction of the signals in the other orthogonal direction. This allows us to single out a two-spot pattern $\mathrm{T}_{1}$ which will serve as the initial state in the switching process. Here, we start with the two-spot pattern oriented in the $x$-direction (direction of intrinsic anisotropy). This selection process is schematically visualized in Fig. 2.

If now a weak (compared to the pump intensity) $y$-polarized control beam with an in-plane momentum on the elastic circle and spatially orthogonal to the initial pattern is applied, stimulated scattering leads to population revival of the corresponding two-spot pattern. Due to the cross-saturation effect, the initial pattern is destabilized and finally switched off completely while the orthogonal pattern reaches a steady state $\mathrm{T}_{2}$ (cf. Fig. 3). After the control beam has been switched off again, the anisotropy leads to re-emergence of the initial $\mathrm{T}_{1}$ pattern while the $\mathrm{T}_{2}$ state vanishes again. If the control beam is too weak so that complete switching may not be possible, the system will remain in a stationary four-spot pattern state $\mathrm{F}$. Based on this scheme, in Ref. Lewandowski17 transistor-like reversible switching was demonstrated including a systematic study of switching times, minimum control power needed, and achievable gain. This was done by numerical simulations of the nonlinear set of equations of motion governing the coherent coupled light field and exciton dynamics in the microcavity system in the two-dimensional QW plane in real space,

[TABLE]

Here, the index $\pm$ denotes the different components in the circular polarization basis, $H_{c}{=}\hbar\omega_{0}{-}\tfrac{\hbar^{2}}{4}(\tfrac{1}{m_{\mathrm{TM}}}{+}\tfrac{1}{m_{\mathrm{TE}}})(\partial_{x}^{2}{+}\partial_{y}^{2})$ is the cavity Hamiltonian and $\epsilon_{0}^{e}$ the flat exciton dispersion. $\gamma_{c}$ and $\gamma_{e}$ represent the photon and exciton loss rates, $\Omega$ describes the photon-exciton coupling, and $H^{\pm}{=}{-}\tfrac{\hbar^{2}}{4}(\tfrac{1}{m_{\mathrm{TM}}}{-}\tfrac{1}{m_{\mathrm{TE}}})(\partial_{x}{\mp}\mathrm{i}\partial_{y})^{2}$ couples the two polarization components due to TE/TM-splitting. The cubic nonlinearities consist of a phase-space filling term $\alpha_{\mathrm{psf}}$, repulsive interaction $T^{++}$ for excitons with parallel spins, and attractive interaction $T^{+-}$ for excitons with opposite spins.

Based on Eqs. (1) with an $x$-linearly polarized continuous wave pump beam at $k\approx 0$ we perform a numerical simulation to demonstrate the basic switching process in the two-dimensional plane of in-plane momenta. Figure 3 shows the photonic component, $|E|^{2}$, in $k$-space without control beam such that the stationary pattern $\mathrm{T}_{1}$ forms in panel (a) and with the control beam on, switching to the stationary pattern $\mathrm{T}_{2}$ in panel (b). The control beam has the same frequency as the pump beam (cl. Fig. 1). System parameters used and details of the calculations are given in Appendix A. Upon switching off the control beam the switching action is reversed and the pattern returns to its original state $T_{1}$. Switching times in the range of ${\sim}100~{}\mathrm{ps}$ are achievable Lewandowski17 .

The focus of the present paper is not on the full numerical simulations of Eqs. (1) and resulting switching performance. Rather, starting from Eqs. (1) we will derive a simplified population competition model for selected modes in $k$-space (details of the derivation are given in Appendix B) to provide further insight into the underlying phase space singularities that dictate the global behavior and solution space of the nonlinear dynamical system studied. In a similar fashion this approach was previously applied to the switching between subsets of a hexagonal pattern PCmodel . We will systematically analyze the existence and stability properties of possible steady states in dependence of the strength of the different involved physical processes for a typical orthogonal switching setup comparable to the one introduced above. To this end, we will construct phase boundaries in representative regions of parameter space and characterize relevant bifurcations. This will lead to a general understanding of crucial parameter dependencies such as the ratio between the control beam strength and the anisotropy, and also of the coexistence of solutions in certain parameter regions and hysteresis behavior. Based on the simplified model, we will also be able to show that the polariton dynamics and optical switching phenomena discussed in the present paper can mathematically be understood based on an generalized Lotka-Volterra model including an external control parameter.

III Population Competition Model

The simplified population competition model discussed in the remainder of the present paper can be derived from the set of equations of motion in Eq. (1) as detailed in Appendix B. It reads

[TABLE]

Here $A_{1}$ and $A_{2}$ are real-valued amplitudes of the two elementary states of the system, i.e., the two orthogonal polariton mode pairs in $k$-space. In case of stationary solutions the phases in the full model (1) become “locked”, allowing us to remove them as dynamical variables from the equations which is suitable for the following steady state analysis. The six dimensionless real-valued positive parameters, $\alpha_{i}$, $\beta_{i}$, $\theta_{i}$ are directly related to the main (up to third order) polariton-polariton scattering processes of the original system as illustrated in Fig. 4. They can be calculated from the physical parameters of the full model (cf. Appendix C). These parameters are intrinsically different for $A_{1}$ and $A_{2}$ due to the polarization dependence and anisotropy. The linear process representing growth of the resonant modes is described by $\alpha_{i}$. Saturation processes are represented by the cubic terms which can be divided into self-saturation ($\beta_{i}$) and cross-saturation ($\theta_{i}$). The external control is described by the inhomogeneity $S$. If we rewrite Eqs. (2) as $(\dot{A}_{1},\dot{A}_{2})^{\mathrm{T}}\equiv(f_{1},f_{2})^{\mathrm{T}}=\mathbf{f}$, steady states are characterized by $\mathbf{f}{=}\mathbf{0}$. Four qualitatively different stationary solutions are possible: i) two-spot pattern $\mathrm{T}_{1}$ with $A_{1}{\neq}0$ and $A_{2}{=}0$, ii) two-spot pattern $\mathrm{T}_{2}$ with $A_{1}{=}0$ and $A_{2}{\neq}0$, iii) four-spot pattern $\mathrm{F}$ with $A_{1}{\neq}0$ and $A_{2}{\neq}0$, and iv) the trivial solution with $A_{1}{=}0$ and $A_{2}{=}0$. Here we are only interested in physical solutions $A_{i}\geq 0$, and therefore the state space is $\mathbb{R}^{2}_{\geq 0}$, i.e. the first quadrant of the $(A_{1},A_{2})$-plane. The linear term ($\alpha$) leads to exponential growth of the corresponding mode pair and the self-saturation ($\beta$) has a stabilizing effect. For these two processes, the equations (2) are automatically decoupled and may be solved separately, resulting in a stable four-spot pattern. However, the cross-saturation terms ($\theta$) couple the two modes and tend to suppress a particular mode pair, favoring the other mode pair, resulting in a two-spot pattern. Additionally, the external control acts as a source term for $A_{2}$. The PC model thus describes the dynamical competition between the three types of possible stationary patterns: $\mathrm{T}_{1}$, $\mathrm{T}_{2}$ and F.

A phase in the PC model is defined as a set consisting of the number of steady states and their stability properties. These sets, i.e. the phase, are functions of the seven parameters and a phase boundary in parameter space indicates the change in the number of steady states and/or their stability. This can only happen at points where at least one eigenvalue of the corresponding Jacobian matrix $J{\equiv}\left(\partial_{A_{j}}f_{i}\right)$ is zero wiggins2006introduction , which is equivalent to the condition $\det{J}{=}0$. We thus find phase boundaries in parameter space at points satisfying $\{\mathbf{f},\det{J}\}{=}\mathbf{0}$.

III.1 Homogeneous Case $S=0$

The homogeneous case of Eq. (2) ($S=0$) has the form of a generalized Lotka-Volterra (GLV) model brenig1988complete with cubic nonlinearities. The transformation $A_{i}\rightarrow A_{i}^{2}\equiv\tilde{A}_{i}$ yields the usual Lotka-Volterra (LV) equations (hofbauer1988theory, ) with quadratic nonlinearities while conserving the $(A_{1},A_{2})$ phase space structure in the positive quadrant hernandez1997lotka . Hence, for $S=0$ the phase portraits of (2) are topologically equivalent to those of the following well-known system

[TABLE]

with $i=1,2$, growth rate vector $\mathbf{r}=2(\alpha_{1},\alpha_{2})$, and community matrix

[TABLE]

which describes self- and cross-interaction. This LV model (3) for interspecific competition has been studied in many different contexts, e.g. ecology wangersky1978lotka , chemistry Hering1990 , economics goodwin1982growth , physics lamb1964theory and it was shown roughgarden1979theory ; hofbauer1988theory that its dynamical behavior is limited to three cases depending on the parameters: i) Coexistence regime for $\mathrm{sgn}(\det{C}){=}{+}1$ (larger self-saturation), ii) bistability regime for $\mathrm{sgn}(\det{C}){=}{-}1$ (larger cross-saturation), iii) dominance regime for sufficiently large/small growth rate ratios. In this regime long-time dynamics are independent of the interaction parameters and always result in extinction of one population. Although the LV system, Eq. (3), is well-known, we present a short discussion of the homogeneous case of the PC model, Eq. (2), and point out the importance for application to the polariton switching dynamics. Since for $S=0$ the system is solvable analytically we obtain explicit expressions for all steady states and phase boundaries: i) $\mathrm{T}_{\mathrm{1}}$: $A_{1}{=}\sqrt{\tfrac{\alpha_{1}}{\beta_{1}}}$, $A_{2}{=}0$, stable for $\tfrac{\alpha_{1}}{\alpha_{2}}>\tfrac{\beta_{1}}{\theta_{2}}$, ii) $\mathrm{T}_{\mathrm{2}}$: $A_{1}{=}0$, $A_{2}{=}\sqrt{\tfrac{\alpha_{2}}{\beta_{2}}}$, stable for $\tfrac{\alpha_{1}}{\alpha_{2}}<\tfrac{\theta_{1}}{\beta_{2}}$, iii) $\mathrm{F}$: $A_{1}{=}\sqrt{\tfrac{\alpha_{2}\theta_{1}-\alpha_{1}\beta_{2}}{\theta_{1}\theta_{2}-\beta_{1}\beta_{2}}}$, $A_{2}{=}\sqrt{\tfrac{\alpha_{1}\theta_{2}-\alpha_{2}\beta_{1}}{\theta_{1}\theta_{2}-\beta_{1}\beta_{2}}}$, only exists for $\tfrac{\beta_{1}}{\theta_{2}}{\lessgtr}\tfrac{\alpha_{1}}{\alpha_{2}}{\lessgtr}\tfrac{\theta_{1}}{\beta_{2}}$ and is stable for $\theta_{1}\theta_{2}{<}\beta_{1}\beta_{2}$, iv) trivial solution: $A_{1}{=}A_{2}{=}0$. We already see that the four-spot solution does not exist in the entire $S{=}0$ plane in contrast to the two-spot pattern. The trivial solution also exists everywhere but is always unstable since the eigenvalues of its Jacobian $\alpha_{1},\alpha_{2}$ are positive. This solution will not be listed hereinafter. For a systematic discussion we introduce a variable anisotropy parameter $\delta\alpha$ in the first equation of (2), $\alpha_{1}\rightarrow\alpha_{1}+\delta\alpha$, and set $\alpha_{1}{=}\alpha_{2}{=}1$. It favors $A_{1/2}$ for $\delta\alpha\gtrless 0$, respectively.

The resulting phase boundaries are shown in Fig. 5 where the homogeneous case ($S{=}0$) is included in the extended region at the bottom of each plot. It shows the structure of the usual LV model. Green (dark gray) / red (light gray) letters mark stable / unstable steady states. For large anisotropy $|\delta\alpha|$ only stable and unstable two-spot patterns are possible. This dominance regime does no longer depend on the interaction parameters $\beta_{i}$ and $\theta_{i}$. The middle region on the other hand is divided into two cases: A stable four-spot pattern for larger self-saturation ($\beta_{1}\beta_{2}{>}\theta_{1}\theta_{2}$, coexistence, upper panel) and two simultaneously stable two-spot patterns for larger cross-saturation ($\beta_{1}\beta_{2}{<}\theta_{1}\theta_{2}$, bistability, lower panel).

The first four cases in Fig. 6 show the corresponding flow given by representative trajectories in state space for the homogeneous case. The system’s dynamics are unambiguous in cases 1, 2, and 3a where only one attractor exists which determines the long time behavior. However, in case 3b two attractors exist and the dynamics now depend on the system’s history leading to a hysteresis effect. The two basins of attraction are marked with different colors blue (dark gray) and orange (light gray). They are separated by the stable manifold of the saddle point F (often called separatrix), defined by the set of points $(A_{1},A_{2})$ which satisfy $(A_{1},A_{2})\rightarrow\mathrm{F}$ for $t\rightarrow\infty$. In contrast to this, the line connecting all three steady states is the unstable manifold of F, consisting of points which satisfy $(A_{1},A_{2})\rightarrow\mathrm{F}$ for $t\rightarrow-\infty$. Any initial point will stay in its region (blue or orange) and end up at the corresponding attractor ($\mathrm{T}_{1}$ or $\mathrm{T}_{2}$). The system thus shows hysteresis behavior which might prevent complete back-switching, and therefore should be avoided for switching purposes, for example by increasing the anisotropy.

The Hartman Grobman theorem wiggins2006introduction ensures that the behavior near hyperbolic equilibrium points is completely determined by its linearization. This no longer holds for non-hyperbolic fixed points which are characterized by the existence of at least one zero real-part eigenvalue of the corresponding Jacobian. Therefore, the behavior at the phase boundaries in parameter space can not be analyzed via linear stability analysis. Instead one can use center manifold theory and normal forms wiggins2006introduction in order to determine the occurring bifurcations. They are transcritical for the usual LV model but different in the case of cubic nonlinearities due to the additional $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetry. For both two-spot solutions symmetric pitchfork (PF) bifurcations occur at the phase boundaries when the anisotropy parameter is changed. They are supercritical for $\beta_{1}\beta_{2}{>}\theta_{1}\theta_{2}$ (larger self-saturation) leading to the coexistence regime with a stable F pattern and subcritical for $\beta_{1}\beta_{2}{<}\theta_{1}\theta_{2}$ (larger cross-saturation) leading to the bistability regime with stable $\mathrm{T}_{1}$ and $\mathrm{T}_{2}$ pattern. Pitchfork bifurcations are typical for dynamical systems with inversion symmetry, here $A_{i}\rightarrow{-}A_{i}$. Since we are only interested in positive solutions in the first quadrant of the phase plane, we either observe a transition from a stable T solution to an unstable T and a stable F solution (supercritical) or the same transition with reversed stability (subcritical).

In conclusion, all phase portraits of the homogeneous case can be completely reduced to the results of the well-known Lotka-Volterra model for two competitive species. Due to the additional inversion symmetry we obtain pitchfork bifurcations at the phase boundaries. The homogeneous case is crucial for the polariton switch because it describes the initial pattern formation and back-switching process in absence of the control beam. Both can only work reliably in regions where only the $\mathrm{T}_{1}$ pattern exists as a single attractor. Therefore, a sufficient minimum anisotropy in favor of $A_{1}$ is needed. We have also seen that if the cross-saturation dominates over the self-saturation, the coexistence of $A_{1}$ and $A_{2}$ is destabilized resulting in the extinction of $A_{i}$ and survival of $A_{i+1}$. A strong interspecific competition thus prevents the coexistence. Otherwise, if the self-saturation is stronger than the cross-saturation, the dominating intraspecific competition promotes coexistence. The parameters used for the full numerical simulations in Sec. II lead to the larger self-saturation case and sufficient anisotropy (cf. Appendix C), providing reliable initial pattern formation and back-switching, since the F pattern is unstable.

III.2 Inhomogeneous Case $S>0$

The inhomogeneous ($S{>}0$) case of the PC model (2) describes the actual switching process induced by an external control beam. This term can also be motivated for other systems where the GLV model is commonly used, e.g. to include constant migration/harvesting in the description of ecological systems or constant influx in chemical reactions. Therefore, it can be interpreted as an extension of the generalized Lotka-Volterra dynamics and the following analysis is of very general interest but has not been investigated before. The influence of constant terms in the usual LV model with quadratic nonlinearities was investigated in Ref. Saputra2010 from a purely mathematical perspective, but in general inhomogeneous population competition models did not receive much attention in the past. Here, we analyze and apply the GLV model with a constant inhomogeneity to the case of the orthogonal switching of two-spot polariton patterns.

Also including inhomogeneity, all steady states and phase boundaries can still be determined analytically due to the system’s simplicity. However, we note that in general solving multivariate polynomial equation systems is a difficult task which can be simplified using algebraic methods (e.g. Gröbner Basis buchberger1 ; buchberger2 ), as was shown in Ref. PCmodel . The non-zero source term $S$ in (2) breaks the inversion symmetry for $A_{2}$ and prevents a $\mathrm{T}_{1}$ solution. This leaves us with only two competing patterns (called $\mathrm{T}_{2}$ & $\mathrm{F}_{2}$) this time around. However, as we will see below it includes the possibility of an additional qualitatively different four-spot pattern solution $\mathrm{F}_{1}$. The control parameter $S$ acts as a constant source term for the $A_{2}$ mode pair, leading to linear growth. Supporting $A_{2}$ means simultaneously suppressing $A_{1}$ due to the cross-saturation effect, leading to a more asymmetric $\mathrm{F}_{2}$ solution in favor of $A_{2}$. A second four-spot pattern, $\mathrm{F}_{1}$, which is in favor of $A_{1}$ replaces $\mathrm{T}_{1}$. The phase boundaries are given explicitly by

[TABLE]

and additionally, for the case of larger cross-saturation ($\beta_{1}\beta_{2}{<}\theta_{1}\theta_{2}$)

[TABLE]

The explicit expression for the steady states $\mathrm{T}_{2}$, $\mathrm{F}_{1}$, and $\mathrm{F}_{2}$ are not given here due to their excessive length. The phase boundaries in the $(\delta\alpha,S)$ parameter space can be seen in Fig. 5 for both the case of higher self- and higher cross-saturation. Starting from the threshold values $\delta\alpha^{*}_{1,2}$ of the homogeneous case, we find continuously changing phase boundaries ($S^{*}_{1,2}$) to higher anisotropy values for increasing control $S$. In the case of higher self-saturation there is only one phase boundary $S^{*}_{1}$ due to the absence of a $\mathrm{T}_{1}$ state. For sufficiently strong control $S$ there is only the stable steady state $\mathrm{T}_{2}$. This region represents the desired outcome of a successful switching process. Its phase boundary determines the minimum control strength required to achieve switching for a given anisotropy value. For lower values of $S$, a stable four-spot pattern arises and depending on the ratio $\tfrac{\theta_{1}\theta_{2}}{\beta_{1}\beta_{2}}$ a second unstable four-spot pattern occurs together with a stability transition of $\mathrm{T}_{2}$.

Again, we can draw the flow in the $(A_{1},A_{2})$ phase plane and mark the basins of attraction as shown in Fig. 6, cases 4, 5, and 6. For $\theta_{1}\theta_{2}{>}\beta_{1}\beta_{2}$ a region with two attractors occurs similar to the $S{=}0$ case, but this time with a stable $\mathrm{F}_{1}$ state instead of $\mathrm{T}_{1}$, which again implies hysteresis behavior. For increasing anisotropy, the unstable $\mathrm{F}_{2}$ state will approach $\mathrm{T}_{2}$ along the unstable manifold until they meet and an unstable $\mathrm{T}_{2}$ and stable $\mathrm{F}_{1}$ state emerge. This corresponds to the approach of the two phase boundaries $S^{*}_{1}$ and $S^{*}_{2}$ and the annihilation of $S^{*}_{2}$, which can only happen for parameter values $1{<}\tfrac{\theta_{1}\theta_{2}}{\beta_{1}\beta_{2}}{<}\tfrac{3}{2}$. Otherwise the two boundaries diverge. Thus, the case of higher cross saturation is divided into two subcases, defined by either the survival or the vanishing of the middle region with two attractors.

Another important effect of the inhomogeneity $S$ is the explicit breaking of the $\mathbb{Z}_{2}$ symmetry in the time evolution of $A_{2}$. In the case of larger cross-saturation a cusp point (C) arises in the two-dimensional bifurcation diagram in Fig. 5(b). The pitchfork bifurcation at $S{=}0$ is replaced by a saddle-node bifurcation for $S{\neq}0$, resulting in the creation/destruction of a stable-unstable pair of F-solutions while varying one of the parameters. In the cusp point two saddle-node phase boundaries, $S_{2}^{*}$ and $-S_{2}^{*}$ (not shown), meet tangentially. The other pitchfork bifurcation ($\delta\alpha_{1}^{*}$) remains ($S_{1}^{*}$) in presence of the inhomogeneity since the $\mathbb{Z}_{2}$ symmetry of $A_{1}$ is still conserved. This is the only bifurcation in the case of larger self-saturation in Fig. 5(a) for $S{\neq}0$. Figure 7 shows the evolution of the equilibrium surfaces for $A_{1}$ and $A_{2}$ in the $(\delta\alpha,S)$ parameter space and the corresponding bifurcations. The pitchfork bifurcations are visible on the $S{=}0$ line. In the case of higher cross-saturation (lower panel) one of the pitchfork bifurcations unfolds into a saddle-node bifurcation with increasing $S$. In the case of higher self-saturation (upper panel) the supercritical pitchfork bifurcation remains stable with increasing $S$.

In conclusion, the main difference for the inhomogeneous case is the absence of a steady $\mathrm{T}_{1}$ state and the existence of an additional stable $\mathrm{F}_{1}$ state in the case of higher cross-saturation. This corresponds to the unfolding of one of the subcritical pitchfork bifurcations into a saddle-node bifurcation. The control term $S$ thus prevents extinction of $A_{2}$ and promotes coexistence respectively. Furthermore, the $\mathrm{F}_{1}$-$\mathrm{T}_{2}$-bistability region vanishes for high anisotropy values, if the parameter condition $1{<}\tfrac{\theta_{1}\theta_{2}}{\beta_{1}\beta_{2}}{<}\tfrac{3}{2}$ is satisfied, resulting in a single remaining phase boundary $S_{1}^{*}$ similar to the case of higher self-saturation. The occurring saddle-node and subcritical pitchfork bifurcations are problematic since they both imply sudden vanishing of a stable fixed point, meaning the system undergoes an abrupt transition to another stable fixed point, i.e. hysteresis is possible. This can only happen in the bistability regime (larger cross-saturation). The numerical simulation of the switching process presented in Section II takes places in the larger self-saturation regime (model parameters calculated in Appendix C), which is advantageous for the switching purpose, since no hysteresis can occur. However, if the external control beam is not sufficiently strong, we observe a stable $\mathrm{F}_{1}$ pattern in the numerical simulations, as predicted by the PC model.

III.3 Remarks

So far, we have discussed steady states, their stability properties, and bifurcations occurring in the solution space of the population competition model in Eq. (2). From a dynamical perspective, we observe critical slowing down scheffer2009early ; tredicce2004critical near the phase boundaries due to the continuously vanishing real part of the Jacobian’s eigenvalue responsible for the bifurcation. This corresponds to the divergence of switching times observed in the numerical simulations in Ref. Lewandowski17 . For example, approaching $S^{*}_{1,2}$ for $\beta_{1}\beta_{2}{\gtrless}\theta_{1}\theta_{2}$ from above results in divergence of the switching time. Similarly, approaching $\delta\alpha^{*}_{2,1}$ for $\beta_{1}\beta_{2}{\gtrless}\theta_{1}\theta_{2}$ from the left side results in divergence of the back-switching time. Hence, to achieve favorable performance, switching should be done for parameters sufficiently far away from the phase boundaries.

We note that in general a nonlinear dynamical system can have periodic solutions which are characterized by closed orbits in the phase plane. Here, we use Dulac’s Criterion (strogatznonlinear, , p. 202) to rule out any periodic solutions. A simplified version reads: The existence of a function $g(A_{1},A_{2})$ with the property that $\nabla{\cdot}(g\mathbf{f})$ is sign definite in the entire considered state space, rules out any closed orbits in this area. If we choose $g{=}\tfrac{1}{A_{1}A_{2}}$, we obtain $\nabla\cdot(g\mathbf{f}){=}{-}2(\beta_{1}\tfrac{A_{1}}{A_{2}}+\beta_{2}\tfrac{A_{2}}{A_{1}}+\tfrac{S}{A_{1}A_{2}^{2}})$, and therefore closed orbits in the positive quadrant are impossible. Another observation is that for $\theta_{1}{=}\theta_{2}{\equiv}\theta$ (symmetric coupling) we can write the system in Eq. (2) as a gradient field $\mathbf{f}{=}\nabla V$ with the potential function

[TABLE]

with $A_{i+2}{=}A_{i}$, $S_{1}{=}0$, and $S_{2}{=}S$. Closed orbits are impossible in gradient systems (strogatznonlinear, , p. 199), however, for the general case $\theta_{1}{\neq}\theta_{2}$ (asymmetric coupling) this argument is no longer applicable. Also, assuming the gradient system defined by the potential function (7) allows us to use the language of catastrophe theory and observe that two of Thom’s seven elementary catastrophes occur, namely folds (denoted as $A_{2}$ in Arnold’s notation) and cusps ($A_{3}$). We further note that our detailed investigation above is limited to destabilizing linearity, stabilizing nonlinearities, and positive control parameter. These conditions match the numerical and experimental observations for the physical system under investigation here. For this case we are able to completely characterize all possible steady states, their corresponding bifurcations, and rule out any other bifurcations involving double zero and pure imaginary eigenvalues of the Jacobian.

IV Conclusions

We have presented a detailed analysis of an all-optical switching concept based on transverse patterns in an interacting polariton system in a planar quantum-well semiconductor microcavity. From the relevant equations based on a microscopic semiconductor theory here we derived a simplified population competition (PC) model describing the system dynamics restricted to selected modes in $k$-space. Interestingly, the resulting rather simple PC model shows all key features of the system dynamics also observed in the full numerical simulations in the parameter range of interest here. In addition to what can be learned from the numerical simulations, the PC model enables us to systematically identify phase boundaries in parameter space and singularities governing the global dynamics of the nonlinear system. Interestingly the rather complicated original system of interacting microcavity polaritons, for the switching phenomenon studied here can be completely characterized by only seven flow portraits in a two-dimensional state space. The model derived has the very generic form of a generalized Lotka-Volterra (GLV) system extended with an inhomogeneity term to achieve external control. Such a system has not been investigated before. Considering the wide-spread use of GLV systems the understanding obtained in the present work is of very general nature and will be similarly applicable also to other fields where external control of population competition is studied such as chemical reactions or population competition in life-sciences.

V Acknowledgments

We gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft through project SCHU 1980/5-2 and through the Heisenberg program (No. 270619725). We further acknowledge valuable discussions with Rolf Binder and a grant for computing time at the Paderborn Center for Parallel Computing ($\mathrm{PC^{2}}$).

Appendix A Numerical Details

We used the following parameters (appropriate for GaAs systems) to numerically simulate the switching presented in section II: $m_{\mathrm{TE}}{=}1.05{\cdot}m_{\mathrm{TM}}{=}0.215~{}\mathrm{meV}\mathrm{ps}^{2}\mu\mathrm{m}^{-2}$, $\Omega{=}6.5~{}\mathrm{meV}$, $\gamma_{c}{=}0.8~{}\mathrm{meV}$, $\gamma_{e}{=}0.2~{}\mathrm{meV}$, $\alpha_{\mathrm{psf}}{=}5.188\cdot 10^{-4}\mu\mathrm{m}^{2};$ $T^{++}{=}-5T^{+-}{=}5.69\cdot 10^{-3}\mathrm{meV}\mu\mathrm{m}^{2}$, $\epsilon_{0}^{e}=1.497~{}\mathrm{eV}$. Coherent cw excitation $2.5~{}\mathrm{meV}$ above the lower polariton branch was used with flat-top Gaussian profiles with intensities $I_{\mathrm{pump}}=2.1\cdot 10^{5}\cdot I_{\mathrm{probe}}=54~{}\mathrm{kW}\mathrm{cm}^{-2}$. The equations (1) were solved on a finite 2D grid in real space using a 4th order Runge-Kutta method with variable step size. In Section III we use parameters $\beta_{1}\beta_{2}/\theta_{1}\theta_{2}{=}1.14$ for the case of larger self-saturation and $\beta_{1}\beta_{2}/\theta_{1}\theta_{2}{=}0.69$ for the case of larger cross-saturation. All bifurcations were determined analytically for the homogeneous case and numerically with the help of MATCONT, a matlab software package, for the inhomogeneous case.

Appendix B Derivation of the PC Model

We follow qualitatively the derivation of the hexagon PC model PCmodel , but here we are including polarization effects for linearly polarized excitation, and therefore considering a different reduced $k$-space. We start with the coupled equations of motion (1) for the exciton polarization and the cavity field and transform them into $k$-space and in the linear polarization basis. Introducing the reduced $k$-space consisting of modes $\{k_{0}{\equiv}0,k_{1},k_{2},k_{3},k_{4}\}$ with relations $k_{3}{=}-k_{1}$ and $k_{4}{=}-k_{2}$ (cf. Figure 8) results in $20$ equations.

We consider an $x$-polarized pump. This leads to the following selection rules PhysRevB.94.045308 :

• •

$\begin{subarray}{c}\mathrm{x}\\ \mathrm{y}\end{subarray}$-polarized probe with $\mathbf{k}\parallel\mathbf{e}_{x}$ excites $\begin{subarray}{c}\mathrm{TM}\\ \mathrm{TE}\end{subarray}$-mode

• •

$\begin{subarray}{c}\mathrm{x}\\ \mathrm{y}\end{subarray}$-polarized probe with $\mathbf{k}\perp\mathbf{e}_{x}$ excites $\begin{subarray}{c}\mathrm{TE}\\ \mathrm{TM}\end{subarray}$-mode

We assume that all dynamical quantities oscillate with pump frequency $\omega$. Removing the phase factor $\mathrm{e}^{-\mathrm{i}\omega t}$ yields a shift of the dispersion by $-\hbar\omega$. We consider all phase-matched scattering processes uo to third order within the reduced k-space. In agreement with the linear stability analysis of the corresponding system reporting a $D_{2}\cong\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetry PhysRevB.77.073302 , we assume an equal excitation of opposite modes in the reduced $k$-space, i.e. $p_{\mathbf{k}_{1}}{=}p_{\mathbf{k}_{3}}{\equiv}p_{1}$ and $p_{\mathbf{k}_{2}}{=}p_{\mathbf{k}_{4}}{\equiv}p_{2}$ (analogously for $E$). We further assume the pump induced densities $E_{0}^{x/y}$ and $p_{0}^{x/y}$ at $k=0$ to be constant. Using a cyclic definition for the two mode pairs $p_{j}{=}p_{j+2}$ with $j{=}1,2$ results in the following phase-matched ($\mathbf{q}{=}\mathbf{k}-\mathbf{k}^{\prime}-\mathbf{k}^{\prime\prime}$) scattering processes $(\mathbf{q},\mathbf{k}^{\prime},\mathbf{k}^{\prime\prime})$ for each mode pair in Eq. (1):

[TABLE]

In our setup the pump is $x$-polarized. Therefore, terms $\propto p_{0}^{y}$ or $\propto E_{0}^{y}$ are omitted. Furthermore, we assume the co-linearly polarized off-axis modes to be very small as for these modes the instability threshold is not reached, i.e. $p_{j}^{x}\ll p_{j}^{y}$ and $E_{j}^{x}\ll E_{j}^{y}$. Hence, we also neglect terms $\propto p_{j}^{x}$ or $\propto E_{j}^{x}$. This leaves us with the following four equations for the $y$-polarized mode pairs:

[TABLE]

If we further assume that the time evolution of $E$ follows adiabatically the evolution of $p$ and set $\partial_{t}E_{j}{\approx}0$, we can write

[TABLE]

with $j{=}0,1,2$. The pump field was also set to zero $E_{\mathrm{pump},j}{=}0$ and will be added later manually. We define $\theta_{j}$ as the phase between $E_{j}^{y}$ and $p_{j}^{y}$ and the ratio of their amplitudes is given by $\tfrac{2\lambda_{j}}{\alpha_{\mathrm{PSF}}\Omega}$. Here, $\delta_{j}$ includes anisotropy effects due to higher density of states for the TE modes and tilting of the pump. So the parameters are given by

[TABLE]

[TABLE]

Finally, we obtain two equations for the two elementary states of the system. We also factorize the exciton field into phase and magnitude, i.e. $p_{j}^{y}{=}\tilde{p}_{j}^{y}\mathrm{e}^{\mathrm{i}\varphi_{j}}$, and split the results into separate equations of motion for magnitude and phase:

[TABLE]

If we define $\phi_{j}{\equiv}\varphi_{j}-\varphi_{0}$, the coefficients are then given by

[TABLE]

Analogously to Ref. PCmodel , we remove the phases as dynamical variables by assuming locked phases and linearization. We define the time-dependent phase as $\phi_{j}(t){\equiv}\delta\phi_{j}(t)+\phi_{j}^{(0)}$ where the locked phases satisfy $K_{j}(\phi_{j}^{(0)}){=}0$ and $\delta\phi_{j}(t)$ is a small deviation. Expanding equations (14) and (15) up to first order in $\delta\phi_{j}$ around $\phi_{j}^{(0)}$ and neglecting terms $\propto\delta\phi_{j}\tilde{p}_{k}^{y2}$ leads to an explicit expression for $\delta\phi_{j}$ which can be substituted back to obtain

[TABLE]

We rewrite these two equations in a shorter form:

[TABLE]

We replace $\tilde{p}_{j}^{y}$ with a product of a dimensionless quantity $\hat{\tilde{p}}_{j}^{y}$ and a characteristic quantity $\tilde{p}_{j,c}^{y}$ which carries the original dimension, i.e. $\tilde{p}_{j}^{y}{=}\hat{\tilde{p}}_{j}^{y}\cdot\tilde{p}_{j,c}^{y}$. We do the same for the independent time variable $t{=}\hat{t}\cdot t_{c}$, so that the derivative changes to $\tfrac{\partial\tilde{p}_{j}^{y}}{\partial t}{=}\tfrac{\tilde{p}_{j,c}^{y}}{t_{c}}\tfrac{\partial\hat{\tilde{p}}_{j}^{y}}{\partial\hat{t}}$. The characteristic values are chosen in a way that the corresponding dimensionless quantities are of magnitude $1$. With the definitions $\hat{t}\equiv t$ and $\hat{\tilde{p}}_{j}^{y}\equiv A_{j}$, we can finally write down the population competition model as

[TABLE]

where we have manually added a control parameter $S$ for the $A_{2}$ mode pair.

Appendix C Calculation of Model Parameters

The model parameters can be calculated from the physical parameters via equation (17). Their values, especially the signs, depend on the specific choice of the locked phases. The system’s physical behavior observed in the full numerical simulations suggest destabilizing linear terms and stabilizing nonlinear terms in the PC model. Chosing phases satisfying this condition, characteristic values $\tilde{p}_{j,c}^{y}{=}1~{}\mu\mathrm{m}^{-1}$ and $t_{c}{=}1~{}\mathrm{ps}$, and anisotropy effects $\delta_{1}{=}0.2~{}\mathrm{meV}$ and $\delta_{2}{=}0~{}\mathrm{meV}$, leads to the following model parameters for the switching simulation presented in section II: $\alpha_{1}{=}0.49$, $\alpha_{2}{=}0.43$, $\beta_{1}{=}0.007$, $\beta_{2}{=}0.01$, $\theta_{1}{=}0.006$, $\theta_{2}{=}0.005$. This corresponds to the case of larger self-saturation. The anisotropy is sufficiently high for the initial pattern-formation and back-switching, and also no hysteresis can occur.