Spontaneous topological transitions in a honeycomb lattice of exciton-polariton condensates due to spin bifurcations

TL;DR

This paper theoretically investigates a honeycomb lattice of exciton-polariton condensates, revealing how varying parameters induce topological phase transitions and diverse magnetic-like orderings in the system.

Contribution

It introduces a model of polaritonic graphene showing spontaneous topological transitions driven by spin bifurcations, a novel insight into polariton lattice phase behavior.

Findings

Identification of ferromagnetic, antiferromagnetic, and valence bond phases.

Demonstration of topological order emergence during phase transitions.

Dependence of phases on pump intensity and coupling strength.

Abstract

We theoretically explore nonresonantly pumped polaritonic graphene, a system consisting of a honeycomb lattice of micropillars in the regime of strong light-matter coupling. We demonstrate that, depending on the parameters of the structure, such as intensity of the pump and coupling strength between the pillars, the system shows rich variety of macroscopic ordering, including analogs of ferromagnetic, antiferromagnetic, and resonant valence bond phases. Transitions between these phases are associated with dramatic reshaping of the spectrum of the system connected with spontaneous appearance of topological order.

Equations

$$i\frac{d{\Psi }_n}{dt}=\left[-\frac{i}{2}g(S_n)-\frac{ϵ+i\gamma }{2}{\hat{\sigma }}_x+\frac{1}{2}\left(\bar{\alpha }{S}_n+\alpha S_n^z{\hat{\sigma }}_z\right)\right]{\Psi }_n$$

$$-\frac{1}{2}\sum _{⟨nm⟩}\left[J+\delta J\left(\cos (2{\phi }_m){\hat{\sigma }}_x+\sin (2{\phi }_m){\hat{\sigma }}_y\right)\right]{\Psi }_m.$$

$$S_n\equiv \frac{|{\psi }_{n+}{|}^2+|{\psi }_{n-}{|}^2}{2},S_n^z\equiv \frac{|{\psi }_{n+}{|}^2-|{\psi }_{n-}{|}^2}{2}.$$

$$\Psi_{n}=\left\{\begin{array}[]{l l}\Psi_{n+1},&\quad\text{if}\ \ S_{n}^{z}=S_{n+1}^{z},\\ -\hat{\sigma}_{x}\Psi_{n+1},&\quad\text{if}\ \ S_{n}^{z}=-S_{n+1}^{z}.\\ \end{array}\right.$$

$$W_{\text{bif}}=\Gamma -\gamma +\eta \frac{(ϵ-n_{\uparrow \downarrow }J{)}^2+{\gamma }^2}{\alpha (ϵ-n_{\uparrow \downarrow }J)},$$

$$H_{\mathbf{k}}=-\frac{1}{2}\left(\begin{array}{cc}0& {\hat{\mathcal{J}}}_{\mathbf{k}}\\ {\hat{\mathcal{J}}}_{\mathbf{k}}^{†}& 0\end{array}\right)+\frac{1}{2}\left(\begin{array}{cc}{\mu }_1\Delta {\hat{\sigma }}_z& 0\\ 0& {\mu }_2\Delta {\hat{\sigma }}_z\end{array}\right).$$

$${\hat{\mathcal{J}}}_{\mathbf{k}}={\hat{\mathcal{J}}}_1+{\hat{\mathcal{J}}}_2{\mathrm{e}}^{-i{k}_1}+{\hat{\mathcal{J}}}_3{\mathrm{e}}^{-i{k}_2},$$

$${\hat{\mathcal{J}}}_m=\left(\begin{array}{cc}J& \delta J{\mathrm{e}}^{-2i{\phi }_m}\\ \delta J{\mathrm{e}}^{2i{\phi }_m}& J\end{array}\right),m=1,2,3.$$

$$H_{\mathbf{k}}=\frac{1}{2}\hat{\mathbb{J}}(\mathbf{k})+\frac{\Delta }{2}\cdot \mathrm{diag}({\mu }_1,{\mu }_2,{\mu }_3,{\mu }_4)\otimes {\hat{\sigma }}_z.$$

$$i\frac{\partial \mathbf{\Psi }}{\partial t}=\left[\hat{\mathbb{J}}+i(W-\Gamma )I_{2N}-(ϵ+i\gamma )I_N\otimes {\hat{\sigma }}_x+\left(\begin{array}{cccc}{\hat{S}}_1& & & \\ & {\hat{S}}_2& & \\ & & \ddots & \\ & & & {\hat{S}}_N\end{array}\right)-i\eta \left(\begin{array}{cccc}{S}_1& & & \\ & S_2& & \\ & & \ddots & \\ & & & S_N\end{array}\right)\otimes {\hat{\sigma }}_0\right]\mathbf{\Psi }.$$

$${\hat{S}}_A=\overline{\alpha }{S}_A{\hat{\sigma }}_0+\alpha S_A^z\sigma ,$$

$$S_A\equiv \frac{|{\psi }_{A+}{|}^2+|{\psi }_{A-}{|}^2}{2},S_A^z\equiv \frac{|{\psi }_{A+}{|}^2-|{\psi }_{A-}{|}^2}{2}.$$

$$\left(\begin{array}{cccc}{\hat{S}}_1& & & \\ & {\hat{S}}_2& & \\ & & \ddots & \\ & & & {\hat{S}}_N\end{array}\right)\equiv \underset{A=1}{\overset{N}{⨁}}{\hat{S}}_A\text{or simply}\underset{A}{⨁}{\hat{S}}_A,$$

$$\left(\begin{array}{cccc}{S}_1& & & \\ & S_2& & \\ & & \ddots & \\ & & & S_N\end{array}\right)\equiv \underset{A}{⨁}{S}_A.$$

$$\mathbf{\Psi }=\left(\sqrt{\rho }{\mathtt{e}}_0+{\mathrm{e}}^{-i\omega t}\mathbf{u}+{\mathrm{e}}^{i{\omega }^{\ast }t}{\mathbf{v}}^{\ast }\right){\mathrm{e}}^{-i{\omega }_{0}t},$$

$${\mathtt{e}}_0=\left(\begin{array}{c}{\mathbf{e}}_1\\ {\mathbf{e}}_2\\ ⋮\\ {\mathbf{e}}_N\end{array}\right)=\left(\begin{array}{c}{e}_1^{+}\\ e_1^{-}\\ e_2^{+}\\ e_2^{-}\\ ⋮\\ e_N^{+}\\ e_N^{-}\end{array}\right)\equiv \underset{j=1}{\overset{N}{⨁}}{\mathbf{e}}_j.$$

$$S_A$$

$$S_A^z$$

$${\hat{S}}_A$$

$$+\frac{\sqrt{\rho }}{2}\left(\begin{array}{cc}\boxed{\begin{array}{rl}& {\mathrm{e}}_{+}\left[\begin{array}{c}(\overline{\alpha }+\alpha )\left(e_A^{+}{v}_A^{+}+e_A^{+\ast }{u}_A^{+}\right)\\ +(\overline{\alpha }-\alpha )\left(e_A^{-}{v}_A^{-}+e_A^{-\ast }{u}_A^{-}\right)\end{array}\right]\\ +& {\mathrm{e}}_{-}\left[\begin{array}{c}(\overline{\alpha }+\alpha )\left(e_A^{+}{u}_A^{+\ast }+e_A^{+\ast }{v}_A^{+\ast }\right)\\ +(\overline{\alpha }-\alpha )\left(e_A^{-}{u}_A^{-\ast }+e_A^{-\ast }{v}_A^{-\ast }\right)\end{array}\right]\end{array}}& 0\\ 0& \boxed{\begin{array}{rl}& {\mathrm{e}}_{+}\left[\begin{array}{c}(\overline{\alpha }-\alpha )\left(e_A^{+}{v}_A^{+}+e_A^{+\ast }{u}_A^{+}\right)\\ +(\overline{\alpha }+\alpha )\left(e_A^{-}{v}_A^{-}+e_A^{-\ast }{u}_A^{-}\right)\end{array}\right]\\ +& {\mathrm{e}}_{-}\left[\begin{array}{c}(\overline{\alpha }-\alpha )\left(e_A^{+}{u}_A^{+\ast }+e_A^{+\ast }{v}_A^{+\ast }\right)\\ +(\overline{\alpha }+\alpha )\left(e_A^{-}{u}_A^{-\ast }+e_A^{-\ast }{v}_A^{-\ast }\right)\end{array}\right]\end{array}}\end{array}\right).$$

$${\hat{S}}_{A0}=\left(\begin{array}{cc}\overline{\alpha }{S}_{A0}+\alpha S_{A0}^z& 0\\ 0& \overline{\alpha }{S}_{A0}-\alpha S_{A0}^z\end{array}\right)=\{\begin{array}{rl}& \uparrow ,\left(\begin{array}{cc}{\alpha }_1& 0\\ 0& {\alpha }_2\end{array}\right),\\ & \downarrow ,\left(\begin{array}{cc}{\alpha }_2& 0\\ 0& {\alpha }_1\end{array}\right).\end{array}$$

$${\hat{S}}_A={\hat{S}}_{A0}+\sqrt{\rho }\left(\begin{array}{cc}\boxed{\begin{array}{rl}& {\mathrm{e}}_{+}\left[{\alpha }_1\left(e_A^{+}{v}_A^{+}+e_A^{+\ast }{u}_A^{+}\right)+{\alpha }_2\left(e_A^{-}{v}_A^{-}+e_A^{-\ast }{u}_A^{-}\right)\right]\\ +& {\mathrm{e}}_{-}\left[{\alpha }_1\left(e_A^{+}{u}_A^{+\ast }+e_A^{+\ast }{v}_A^{+\ast }\right)+{\alpha }_2\left(e_A^{-}{u}_A^{-\ast }+e_A^{-\ast }{v}_A^{-\ast }\right)\right]\end{array}}& 0\\ 0& \boxed{\begin{array}{rl}& {\mathrm{e}}_{+}\left[{\alpha }_2\left(e_A^{+}{v}_A^{+}+e_A^{+\ast }{u}_A^{+}\right)+{\alpha }_1\left(e_A^{-}{v}_A^{-}+e_A^{-\ast }{u}_A^{-}\right)\right]\\ +& {\mathrm{e}}_{-}\left[{\alpha }_2\left(e_A^{+}{u}_A^{+\ast }+e_A^{+\ast }{v}_A^{+\ast }\right)+{\alpha }_1\left(e_A^{-}{u}_A^{-\ast }+e_A^{-\ast }{v}_A^{-\ast }\right)\right]\end{array}}\end{array}\right)$$

$${\hat{S}}_A={\hat{S}}_{A0}+$$

$$+$$

$$={\hat{S}}_{A0}+{\hat{S}}_{A+}+{\hat{S}}_{A-}.$$

$$\hat{H}$$

$$+{\mathrm{e}}_{+}\sqrt{\rho }\left[\underset{A}{⨁}{\hat{S}}_{A+}-\frac{i\eta }{2}\underset{A}{⨁}\left(e_A^{+}{v}_A^{+}+e_A^{+\ast }{u}_A^{+}+e_A^{-}{v}_A^{-}+e_A^{-\ast }{u}_A^{-}\right)\otimes {\hat{\sigma }}_0\right]$$

$$+{\mathrm{e}}_{-}\sqrt{\rho }\left[\underset{A}{⨁}{\hat{S}}_{A-}-\frac{i\eta }{2}\underset{A}{⨁}\left(e_A^{+}{u}_A^{+\ast }+e_A^{+\ast }{v}_A^{+\ast }+e_A^{-}{u}_A^{-\ast }+e_A^{-\ast }{v}_A^{-\ast }\right)\otimes {\hat{\sigma }}_0\right]$$

$$={\hat{H}}_0$$

$$+{\mathrm{e}}_{-}\sqrt{\rho }{\hat{\mathbb{S}}}_{+}^{\ast }(-\eta )={\hat{H}}_0+{\hat{\mathbb{S}}}_{+}(\eta )+{\hat{\mathbb{S}}}_{-}(\eta ).$$

$$\{\begin{array}{rl}\omega \mathbf{u}& =\left({\hat{H}}_0-{\omega }_0{I}_{2N}\right)\mathbf{u}+\rho {\hat{\mathbb{S}}}_{+}{\mathtt{e}}_0,\\ \omega \mathbf{v}& =-\left({\hat{H}}_0^{\ast }-{\omega }_0{I}_{2N}\right)\mathbf{v}-\rho {\hat{\mathbb{S}}}_{-}^{\ast }{\mathtt{e}}_0^{\ast }.\end{array}$$

$${\hat{\mathbb{S}}}_{+}{\mathtt{e}}_0$$

$$=\underset{A}{⨁}\left(\begin{array}{c}\left({\alpha }_1-\frac{i\eta }{2}\right)\left(e_A^{+}{v}_A^{+}+e_A^{+\ast }{u}_A^{+}\right)e_A^{+}+\left({\alpha }_2-\frac{i\eta }{2}\right)\text{cancelto}0\left(e_A^{-}{v}_A^{-}+e_A^{-\ast }{u}_A^{-}\right)e_A^{+}\\ \left({\alpha }_2-\frac{i\eta }{2}\right)\text{cancelto}0\left(e_A^{+}{v}_A^{+}+e_A^{+\ast }{u}_A^{+}\right)e_A^{-}+\left({\alpha }_1-\frac{i\eta }{2}\right)\left(e_A^{-}{v}_A^{-}+e_A^{-\ast }{u}_A^{-}\right)e_A^{-}\end{array}\right)$$

$$=\left({\alpha }_1-\frac{i\eta }{2}\right)\cdot \left\{\underset{{\hat{\mathbb{E}}}_0^{\prime }}{\underbrace{\left[\underset{A}{⨁}\left(\begin{array}{cc}{\left|e_A^{+}\right|}^2& 0\\ 0& {\left|e_A^{-}\right|}^2\end{array}\right)\right]}}\mathbf{u}+\underset{{\hat{\mathbb{E}}}_0}{\underbrace{\left[\underset{A}{⨁}\left(\begin{array}{cc}{\left(e_A^{+}\right)}^2& 0\\ 0& {\left(e_A^{-}\right)}^2\end{array}\right)\right]}}\mathbf{v}\right\}.$$

$${\hat{\mathbb{S}}}_{-}^{\ast }{\mathtt{e}}_0^{\ast }={\hat{\mathbb{S}}}_{+}(-\eta ){\mathtt{e}}_0^{\ast }=\left({\alpha }_1+\frac{i\eta }{2}\right)\cdot \left\{{\hat{\mathbb{E}}}_0^{\ast }\mathbf{u}+{\hat{\mathbb{E}}}_0^{\prime }\mathbf{v}\right\}.$$

$$\{\begin{array}{rl}\omega \mathbf{u}& =\left({\hat{H}}_0-{\omega }_0{I}_{2N}\right)\mathbf{u}+\rho \left({\alpha }_1-\frac{i\eta }{2}\right)\left({\hat{\mathbb{E}}}_0^{\prime }\mathbf{u}+{\hat{\mathbb{E}}}_0\mathbf{v}\right),\\ \omega \mathbf{v}& =-\left({\hat{H}}_0^{\ast }-{\omega }_0{I}_{2N}\right)\mathbf{v}-\rho \left({\alpha }_1+\frac{i\eta }{2}\right)\left({\hat{\mathbb{E}}}_0^{\ast }\mathbf{u}+{\hat{\mathbb{E}}}_0^{\prime }\mathbf{v}\right),\end{array}$$

$$\omega \left(\begin{array}{c}\mathbf{u}\\ \mathbf{v}\end{array}\right)=\left(\begin{array}{cc}{\hat{H}}_0-{\omega }_0{I}_{2N}-\rho \left(\frac{i\eta }{2}-{\alpha }_1\right){\hat{\mathbb{E}}}_0^{\prime }& -\rho \left(\frac{i\eta }{2}-{\alpha }_1\right){\hat{\mathbb{E}}}_0\\ -\rho \left(\frac{i\eta }{2}+{\alpha }_1\right){\hat{\mathbb{E}}}_0^{\ast }& -\left({\hat{H}}_0^{\ast }-{\omega }_0{I}_{2N}\right)-\rho \left(\frac{i\eta }{2}+{\alpha }_1\right){\hat{\mathbb{E}}}_0^{\prime }\end{array}\right)\left(\begin{array}{c}\mathbf{u}\\ \mathbf{v}\end{array}\right).$$

$${\omega }_0=\frac{⟨{\mathbf{\Psi }}_0|{\hat{H}}_0|{\mathbf{\Psi }}_0⟩}{⟨{\mathbf{\Psi }}_0|{\mathbf{\Psi }}_0⟩},$$

Full text

††preprint: APS/123-QED††preprint: APS/123-QED

Introduction. Recent decades have witnessed a shift in attention from both the condensed matter and optical communities in investigation of the properties of bulk materials to instead the properties of their interfaces. It is now well known that there exists a particular class of materials with inverted structure of the bands, which possess protected states propagating on the system surface, referred to as topological insulators Hasan and Kane (2010); Qi and Zhang (2011); Bansil et al. (2016). The energy of these edge states lies within the bandgap, and thus they are protected with respect to scattering into the bulk. Depending on the dimensionality of a system, one should distinguish between three-dimensional (3D) topological insulators where topological states appear on the two-dimensional (2D) surface boundary of the bulk Kane and Mele (2005); Fu and Kane (2007); Moore and Balents (2007); Roy (2009); Xia et al. (2009); Zhang et al. (2009); Chen et al. (2009), and 2D topological insulators where chiral 1D channels form on the system boundary Bernevig et al. (2006); König et al. (2007).

It has been recently shown that optical analogs of topologically nontrivial phases, in 2D systems, may arise in purely photonic structures Haldane and Raghu (2008); Wang et al. (2009); Khanikaev et al. (2013) or when driven into the strong light-matter coupling regime, where hybrid quasiparticles known as cavity exciton-polaritons (from here on polaritons) are formed Karzig et al. (2015); Bardyn et al. (2015); Nalitov et al. (2015); Gulevich et al. (2016); Kozin et al. (2018); Klembt et al. (2018). Polaritons combine the advantages of photons, such as extremely low effective mass and long coherence length, with those of excitons, namely the possibility of control by external electric and magnetic fields together with strong nonlinear response stemming from interparticle interactions. The spin structure of the exciton (or rather, the polariton) is then directly related to the circular polarization degree of the cavity photonic mode. Such a union then leads to a rich interplay between nonlinear and topological properties Bleu et al. (2016); Bardyn et al. (2016); Gulevich et al. (2017); Sigurdsson et al. (2017a); Kartashov and Skryabin (2017); Bleu et al. (2017); Mandal et al. (2019) with optical lasing in topologically protected edge modes St-Jean et al. (2017); Klembt et al. (2018).

The vast majority of current proposals on 2D polariton topological insulators are based on $\mathbb{Z}$ (or Chern) insulators with the following requirements: Polaritons should be placed into a 2D lattice of a particular symmetry allowing the appearance of Dirac points in the Brillouin zone where the bands touch each other. Examples are honeycomb Nalitov et al. (2015); Klembt et al. (2018); Kartashov and Skryabin (2019) and Kagome Gulevich et al. (2016, 2017); Sigurdsson et al. (2017a) lattices which can be obtained either by controllable etching of a planar microcavity or by using spatial light modulator to control the profile of the external optical pump. The band inversion and opening of the topological gap is then achieved by cumulative action of the TE-TM splitting of the photonic mode and Zeeman splitting of the excitonic mode induced by the application of an external magnetic field Nalitov et al. (2015). However, in conventional semiconductor materials excitonic $g$-factors are extremely small, and one needs magnetic fields of tens of Tesla to open the topological bandgap of at least several meV. The situation can be potentially improved by using diluted magnetic microcavities Brunetti et al. (2006); Król et al. (2018). However, the technology of producing a high quality patterned semimagnetic cavity is still only in its initial stages.

In the present paper, we develop an alternative approach for the realization of the quantum anomalous Hall effect in a 2D polariton $\mathbb{Z}$-topological insulator without application of any external magnetic fields. Our idea is based on the concept of the spontaneous spin bifurcation in a system of interacting polariton condensates forming a net magnetic polarization, first proposed in Ref. [Ohadi et al., 2015] and developed further in Refs. [Ohadi et al., 2016; Dreismann et al., 2016; Ohadi et al., 2017; Sigurdsson et al., 2017b]. Here we consider a honeycomb polariton condensate lattice (polariton graphene) under nonresonant pumping. We demonstrate that, depending on the pump intensity and coupling strength between the nodes (condensates) of the lattice, the spin bifurcation mechanism can result in spontaneous formation of distinct spin-ordered lattice phases analogous to ferromagnetic (FM), antiferromagnetic (AFM), and resonance valence bond states. Transition between different phases is associated with cardinal reshaping of the spectrum of the system excitations and spontaneous appearance/disappearance of topological order.

Spin bifurcations and phase transitions in polariton graphene. A lattice of driven-dissipative connected polariton condensates is conventionally modeled, in the tight-binding picture, with a set of generalized Gross-Pitaevskii equations for the spinor order parameters $\Psi_{n}=(\psi_{n+},\psi_{n-})^{\text{T}}$, corresponding to spin-up and spin-down polaritons at the $n$-th site,

[TABLE]

Here the summation is taken over the nearest neighbors, $\varphi_{m}$ are the angles of links connecting the neighboring sites $n$ and $m$ of the honeycomb lattice. We define the $n$-th node particle population $S_{n}$ and $z$-component of the condensate pseudospin $S_{n}^{z}$ as:

[TABLE]

We also define an effective decay rate $g(S_{n})=\eta S_{n}+\Gamma-W$ with $\Gamma$ being the polariton decay rate, $W$ the replenishment rate of the condensate non-polarized incoherent pump, and $\eta$ is the gain-saturation nonlinearity. The constants $\epsilon$ and $\gamma$ define the splitting of the $XY$-polarized states in both energy and decay, respectively, due to the inherent cavity birefringence, and $\bar{\alpha}=\alpha_{1}+\alpha_{2}$ and $\alpha=\alpha_{1}-\alpha_{2}$ are spin-anisotropic interaction parameters. Finally, $J>\delta J$ are spin conserving and non-conserving (TE-TM splitting) tunneling rates of polaritons between nodes respectively.

The condensation threshold of the system is defined as the point where an eigenvalue of the linearized Eq. (1) obtains a positive imaginary component due to increase of the laser power $W$, leading to the triggering of the stimulated bosonic scattering into the condensed state at $W_{\text{cond}}=\Gamma-\gamma$. Due to the splitting $\gamma$ in the lifetimes of the linear polarized states the condensate first forms an, in-phase, $Y$-polarized state, i.e., $\Psi_{n}=\Psi_{n+1}\propto(1,-1)^{T}$ (white area in Fig. 1a). This $Y$-polarized state, however, becomes unstable at higher pumping powers and undergoes a bifurcation into a state with a high degree of circular polarization at individual nodes Ohadi et al. (2015, 2017).

We begin our consideration by presenting a class of stationary solutions which minimize the spin bifurcation threshold Sigurdsson et al. (2017b),

[TABLE]

The ansatz above describes in-phase FM bonds and anti-phase AFM bonds between nearest neighbors respectively. Plugging Eq. (3) into Eq. (1), and setting the condition that all nodes have the same number of co- and counterpolarized nearest neighbours (equivalence criteria), the coupled set of the equations of motion reduce to a single equation with a bifurcation threshold,

[TABLE]

where $n_{\uparrow\downarrow}$ denotes the number of AFM neighbors.

In Fig. 1a we plot the minimum of Eq. (4) as a function of coupling strength $J$ (red curve) neglecting TE-TM splitting. The cusps in the red curve indicate that the lowest bifurcation point is shared between two distinct spin phases. The four spin phases of interest, verified by numerical integration of Eq. (1), are shown in Fig. 1b(i-iv) and are labeled AFM, Dipole, Stripe, and FM phases respectively.

We point out that the bifurcation threshold for AFM and FM phases is invariant of $\delta J$ whereas for dipole and stripe phases, strictly speaking, this is not the case as for them the ansatz given by the Eq. (3) should be modified. However, given that $\delta J/J\ll 1$ (which is usually the case in micropillar structures in standard semiconductor microcavities), it is reasonable to infer that $W_{\text{bif}}$ is only weakly affected by $\delta J$. Thus, the calculated red curve shown in the Fig. 1a serves as a good indicator for the bifurcation threshold of these nontrivial states in the presence of TE-TM splitting. We have performed numerical calculations of Eq. (1) that verify that this is indeed the case.

Band structure and topological states. In the following we discuss the excitation spectra of the stable spin bifurcated condensate configurations. We employ an effective field model, treating the effect of spin-polarized lattice nodes by introduction of local $z$-directed (out of the cavity plane) magnetic fields. The idea is based on expansion of the Bogoliubov dispersion in the two spin components $\sqrt{E_{\pm}(E_{\pm}+U_{\pm})}\approx E_{\pm}+U_{\pm}/2$ with $U_{\pm}\approx\alpha_{1}|\psi_{n\pm}|^{2}$ being the self interaction energy. It yields the effective magnetic field magnitude $\Delta=\alpha_{1}|S_{n}^{z}|$ and the validity region of the approximation as $E_{\pm}\gg U_{\pm}$. The effective field model allows clear analytical investigation of the topological properties of the lowest band gap opening at $E\sim J$ as long as $J\gg\Delta$. Furthermore, as the pseudospin strength is given by $|S_{n}^{z}|=S_{n}=(W-\Gamma)/\eta$ for the fully polarized condensates (see Eq. (S17) in [Sup, ]) one can adjust the strength of $\Delta$ such that effects due to birefringence are negligible. The effective field model then obeys $J\gg\Delta\gg\epsilon$. A more rigorous Bogoliubov treatment is addressed in Ref. [Sup, ].

Firstly, to examine the band structures of the AFM (i) and FM (iv) configurations shown in Fig. 1b, we scrutinize the following $4{\times}4$ tight-binding Hamiltonian:

[TABLE]

The above Hamiltonian is written in the basis of the bispinor $|\text{A}+,\text{A}-,\text{B}+,\text{B}-\rangle$ inner-cell states, where A and B indicate the graphene sublattices, and “$+(-)$” specifies right (left) circular polarization (i.e., the spin of the polaritons). The total Hamiltonian in momentum space is then written $\hat{H}=\sum_{\mathbf{k}}|\mathbf{k}\rangle\langle\mathbf{k}|\otimes H_{\mathbf{k}}$. The first term in Eq. (5) corresponds to the polaritonic graphene with TE-TM splitting. The $2{\times}2$ operator $\hat{\mathcal{J}}_{\mathbf{k}}$, dependent on the quasi-wavevector $\mathbf{k}=(k_{1},k_{1})$, is written,

[TABLE]

with

[TABLE]

and $J$ and $\delta J$ as in Eq. (1). Both $k_{1}$ and $k_{2}$ can be chosen to vary from $-\pi$ to $\pi$ and to cover the whole Brillouin zone. In the second term of Eq. (5), the on-diagonal blocks $\mu_{1(2)}\Delta\,\sigma_{z}$ serve to account for the excitations in the magnetic patterns depicted in Fig. 1b. The magnitude of the Zeeman splitting induced by the polarized condensate reads $\Delta=\alpha_{1}S$ where we have omitted the index $n$ in $S_{n}$ since the condensate density is taken equal at each lattice site. Here $\hat{\sigma}_{z}$ is the $z$-Pauli matrix, and the coefficients $\mu_{1}$ and $\mu_{2}$ define FM ($\mu_{1}=\mu_{2}=1$) and AFM ($\mu_{1}=-\mu_{2}=1$) phases. The translation basis vectors ($\mathbf{a}_{1,2}$) in real space are chosen to be conventional for the graphene lattice.

Since FM phase corresponds to the case of a uniform, out-of-plane, external magnetic field, a gap opens (see Fig. 2d) between the bands, characterized by different Chern numbers, and bridged by chiral edge states Nalitov et al. (2015). Fig. 2d(i) displays the band structure of this phase, which is numerically obtained for $\delta J=2/3$, $\Delta=5/3$, with energy counted in units of $J$. We used the convention for which $K$ and $K^{\prime}$ points in the first Brillouin zone are positioned at ($k_{1},k_{2}$) equal to ($2\pi/3,-2\pi/3$), and vice versa, $\Gamma$ is placed at the origin, and $M$ at ($\pi,\pi$). Fig. 2d(ii) shows a slice of the band structure along $\Gamma\to M\to K\to\Gamma$ pathway shown by the green solid line at the bottom of Fig. 2d(i). Figure 2d(iii) illustrates the band structure of the graphene stripe (not to be confused with the ”stripe” phase in Fig. 1b) with zigzag edges revealing the topologically protected edge states, marked by red and blue colored lines for each edge respectively.

For the AFM phase, it can be shown that a gap opens when $\Delta>3\delta J$, the gap value being $E_{g}=2(\Delta-3\delta J)$ Ohadi et al. (2015). The gapped spectrum in the AFM lattice is characterized by trivial topology (see Fig. 2a) and no edge states connecting the bulk bands. In Fig. 1a the dotted white area in the AFM phase corresponds to the bands touching, whereas the yellow area to the bands being gapped. Here, the boundary separating the two regimes in Fig. 11a is calculated for the case of strong TE-TM splitting $\delta J/J=0.5$. We point out that that gapped-ungapped AFM boundary coincides with the $W_{\text{bif}}$ boundary (red curve in Fig. 1a) when $\delta J\to 0$. The AFM edge states shown in Fig. 2a(iii), separated from the bands and marked by red and blue colors, are doubly degenerate and are not topologically protected. Figures 3(a,d) show phase diagrams for the AFM and FM spin phases as a function of effective Zeeman splitting and TE-TM splitting. The dark orange line marks the boundary between gapless and gapped phases. The points A and D correspond to Fig. 2a and d, respectively. Figure 3e shows the band structures of FM-phase ribbons and serves to illustrate the effect of Chern numbers of the bands (indicated in Fig. 3d) on the dispersion of the chiral edge states.

To examine the band structure of the stripe and dipole phases, the following $8{\times}8$ Hamiltonian should be constructed in reciprocal $(k_{1},k_{2})$ space:

[TABLE]

Note, that a unit cell in these two cases differs from that corresponding to AFM and FM phases and should be constructed as a pair of graphene unit cells taken successively, with the translational vectors being $\mathbf{a}_{1}$ and $2\mathbf{a}_{2}$. The first term in Eq. (8), $\widehat{\mathbb{J}}(\mathbf{k})$, represents polaritonic graphene with TE-TM splitting and is specified in the Supplementary material Sup , the second term is responsible for the magnetic patterns: for the dipole phase one sets $\mu_{1}=\mu_{4}=1$, $\mu_{2}=\mu_{3}=-1$, for the stripe phase $\mu_{1}=\mu_{2}=1$, $\mu_{3}=\mu_{4}=-1$. Both phases are characterized by topologically trivial band structure and edge states localized within the bulk (see Fig. 2b,c). Figure 3(b,c) shows phase diagram for dipole and stripe phases with points B and C corresponding to the Fig. 2b,c. In the trivial case $\delta J=0$ one arrives to the gap opening condition $\Delta>J$ and $\Delta>\sqrt{3}J$ for dipole and stripe spin phases respectively which is plotted in Fig. 1a indicating that gap opening only takes place at higher condensate densities (i.e., higher excitation powers).

To investigate the FM phase in more detail, we performed Bogoliubov linearization of Eq. (1), generalized to all four spin patterns, and applied it to the FM one (see [Sup, ] for more details). Figure 4 demonstrates (A) the topological phase diagram for the effective spin-orbit interaction strength $\delta J/J$ and the effective interaction energy $\Delta$, and (B, C) eigenenergy curves along $\Gamma\to M\to K\to\Gamma$ pathway. In addition, green dotted and red dashed curves show the dispersions of Eq. (5) without and with effective magnetic field respectively, with the corresponding value $\Delta$. Note that the bogolon dispersion overlaps with the one obtained in the effective magnetic field approximation in the region of high energies $E\gg\Delta$ and expectedly deviates for low energies (in the vicinity of $\Gamma$ point). This shows that the effective magnetic field model accurately describes the linear excitations of the condensate in the vicinity of the $K$-point. Moreover, Fig. 4(A) clearly resembles Fig. 3(D).

Conclusions. We have proposed an experimentally friendly geometry for realization of an optical $\mathbb{Z}$-topological insulator based on polaritonic graphene in the spin bifurcation regime. Differently from previous works, our proposal does not require application of an external magnetic field and the topological order appears spontaneously in the mean field picture through many-particle interactions under non-resonant and non-polarized pumping. The proposed effective field model due to interactions produces the topological phase diagram of the system and is consistent with the Bogoliubov model in a wide range of parameters.

Acknowledgements. We thank Dr. O. Kyriienko for valuable discussions. The work was supported by Russian Science Foundation (Project No. 18-72-10110). A.V.N. and I.A.S. acknowledge support from Horizon2020 RISE project COEXAN. H.S. acknowledges support from UK’s Engineering and Physical Sciences Research Council (grant EP/M025330/1 on Hybrid Polaritonics). A.V.N. acknowledges support from Icelandic Research Fund, Grant No. 196301-051.

S1 Bogoliubov theory of spin-bifurcated polariton condensates in a honeycomb lattice

We start with the following generalized Gross-Pitaevskii equation as given by Eq. 1 in main text. It describes the bulk honeycomb lattice of coherent spin-up ($+$) and spin-down ($-$) polaritons belonging to a condensate order parameter $\mathbf{\Psi}=(\psi_{1+},\psi_{1-},\psi_{2+},\psi_{2-},\dots)^{T}$ where $T$ denotes the transpose:

[TABLE]

Here we have omitted a factor $1/2$ appearing on the RHS in Eq. 1 in the main text for more transparent derivation. $\widehat{\mathbb{J}}$ incorporates all the hoppings between the sites of the lattice, $(W-\Gamma)$ is incoherent pumping-to-losses imbalance, $\epsilon+i\gamma$ is the complex splitting in the linear polarizations, $I_{2N}$ is the identity matrix of dimension $2N$ where $N$ is the number of condensates (sites) in the lattice which, for the case of $N=1$, we define $I_{2}\equiv\widehat{\sigma}_{0}$, and $\left\{\widehat{\sigma}_{\lambda}\right\}_{\lambda=x,y,z}$ are the Pauli matrices. The “spin” operator $\widehat{S}_{A}$ is written,

[TABLE]

where $\bar{\alpha}=\alpha_{1}+\alpha_{2}$ and $\alpha=\alpha_{1}-\alpha_{2}$ are spin-anisotropic interaction parameters and

[TABLE]

The block-diagonal $\widehat{S}$- and $S$-operators are further designated as direct sums of the corresponding operator subspaces:

[TABLE]

The Bogoliubov elementary excitation (bogolon) spectrum is obtained by the standard procedure of linearizing the Eq. (S1) around the steady state condensates $\mathbf{\Psi}_{0}=\sqrt{\rho}\mathbf{\mathtt{e}}_{0}e^{-i\omega_{0}t}$. We therefore seek $\mathbf{\Psi}$ in the form,

[TABLE]

where the condensates density $\rho$ is taken to be uniform across the lattice, $\omega_{0}$ is the chemical potential, and $\mathbf{\mathtt{e}}_{0}$ captures their polarization,

[TABLE]

Here, $\mathbf{u}$ and $\mathbf{v}$ are the amplitudes of the excitaions, of the same structure as $\mathbf{\mathtt{e}}_{0}$ (S7). Hereafter, $\mathrm{e}^{-i\omega t}$ and $\mathrm{e}^{i\omega^{*}t}$ are denoted as $\mathrm{e}_{+}$ and $\mathrm{e}_{-}$, respectively.

S2 Spins and $\widehat{S}$-operators

For site $A$ we write,

[TABLE]

We keep only linear order terms of $v_{A}^{\pm},u_{A}^{\pm}$. $S_{A0}$ and $S_{A0}^{z}$ are defined as the normalized total spin and its $z$-projection corresponding to polarization $\mathbf{e}_{A}$: $S_{A0}=\left(\left|e_{A}^{+}\right|^{2}+\left|e_{A}^{-}\right|^{2}\right)/2$ and $S_{A0}=\left(\left|e_{A}^{+}\right|^{2}-\left|e_{A}^{-}\right|^{2}\right)/2$.

Consider $\widehat{S}$-operators.

[TABLE]

S3 Hamiltonian with respect to $\mathrm{e}_{+}$ and $\mathrm{e}_{-}$

We now decompose Eq. (S1), separating the terms with $\mathrm{e}_{+}$ and $\mathrm{e}_{-}$ multipliers.

[TABLE]

Then, we equate the corresponding coefficients in front of $\mathrm{e}_{+}$ and $\mathrm{e}_{-}$:

[TABLE]

Further we reveal the action of operator $\widehat{\mathbb{S}}_{+}$ to $\mathbf{\mathtt{e}}_{0}$, as it contains $\mathbf{u}$ and $\mathbf{v}$ terms.

As we consider totally polarized condensate, $\mathbf{e}_{A}=(\mathrm{e}^{i\phi_{A}},0)^{T}$ or $(0,\mathrm{e}^{i\phi_{A}})^{T}$ with unknown phase $\phi_{A}$ we satisfy the realtion $e_{A}^{+}\,e_{A}^{-}=0$. From Eq. (S7) and (S10d), we get after some algebra

[TABLE]

Inserting (S12) into (S11), we arrive at

[TABLE]

and finally

[TABLE]

Each position in the matrix (S14) is of $(2N{\times}2N)$ dimension, so the final matrix to be diagonalized is of size $(4N{\times}4N)$.

S4 Extracting $\rho$ and $\omega_{0}$ for the stationary condensates

Prior to the diagonalization of Eq. (S14), we should obtain quantities $\omega_{0}$ and $\rho$. To find them, we use the equality

[TABLE]

where $\widehat{H}_{0}$ is defined in Eq. (S10a). Then,

[TABLE]

The density $\rho$ is derived from the requirement of the real-value essence of $\omega_{0}$:

[TABLE]

Thus,

[TABLE]

For $\delta J\approx 0$, the first term can be evaluated due to $\langle\mathbf{\mathtt{e}}_{0}|\widehat{\mathbb{J}}|\mathbf{\mathtt{e}}_{0}\rangle=-3NJ,-2NJ,-NJ$, and [math] for the ferromagnetic, stripe, dipole, and antiferromagnetic phases, respectively.

S5 Band structures of the bogolons

In order to get the bogolon dispersion, we transform the external spatial degrees of freedom into the lattice momentum $\mathbf{k}$ (reciprocal) basis. In Eq. (S14), we substitute $\mathbf{u}$ and $\mathbf{v}$ with separated intra-cell and inter-cell degrees of freedom and Fourier-transform along the inter-cell coordinate. The bulk Hamiltonian is then written,

[TABLE]

The number of sites per unit cell is denoted as $\mathcal{N}$. $\widehat{\mathbb{E}}_{0}$, $\widehat{\mathbb{E}}_{0}^{\prime}$, and $\widehat{\mathbb{E}}_{0}^{*}$ with $\mathcal{N}$ positioned as an argument are taken over a single unit cell, i.e. with the summation $\bigoplus\limits_{A}$ limited to sites $A$ within a single unit cell.

Substituting the density (S17) into $\widehat{H}_{0}$, we arrive at

[TABLE]

To present the explicit form of the bulk Hamiltonian (S19), we denote the hopping $2{\times}2$ blocks by $\widehat{\mathcal{J}}_{i}$:

[TABLE]

with the set of angles $\varphi$ chosen as $(0,2\pi/3,4\pi/3)$. Then,

[TABLE]

where $\widehat{\mathbb{O}}$ is the $2\times 2$ zero matrix. The blocks $\widehat{S}_{A0}$ for each distinct polarization pattern belonging to the spin-bifurcated condensates (see Fig. 1b in main text) are written:

[TABLE]

S6 Connection with the model of effective magnetic field

We will consider the scenario where the in-plane magnetic field, caused by birefringence in the cavity and/or its distributed-Bragg-reflectors, is much smaller than the condensate induced out-of-plane magnetic field, i.e. $\alpha\rho\gg\epsilon$ in frequency units. This can be achieved by considering Eq. (S17) which dictates that $\rho$ can be increased by simply tuning the ratio of pump-decay ($W-\Gamma$) against the saturation rate $(\eta)$. If the saturation rate $\eta$ is taken to be smaller than then same-spin polariton-polariton interaction strength $\alpha_{1}$ then we have,

[TABLE]

The prefactor $1/2$ originally omitted from Eq. (S1) is now returned back. Equation (S23) then describes the spinor bogolon dispersion experiencing an effective out-of-plane magnetic field. Results from diagonalizing the above Hamiltonian are presented in the main text (Fig. 4).

In the linear regime ($\rho\simeq 0$) there is no effective out-of-plane magnetic field and the bogolon dispersions become analogous to the dispersions of $\widehat{\mathbb{J}}(\mathbf{k})/2$ graphene Hamiltonian. The effective-magnetic-field model dispersions presented in Figures 2 and 3 in the main text are then the result of diagonalization of

[TABLE]

Within the notations, the magnitude of the magnetic field is $\rho\alpha/2$, in accordance with the definition of $\Delta$ in the main text.