Multi-orbital tight binding model for cavity-polariton lattices

TL;DR

This paper introduces a minimal parameter tight-binding model for exciton-polariton lattices that accurately reproduces experimental band structures, highlighting effects of non-orthogonality, interactions, and spin-orbit coupling.

Contribution

The paper develops a novel tight-binding model based on non-orthogonal orbitals to describe polariton lattice band structures with high fidelity.

Findings

Model accurately reproduces experimental band structures.

Non-orthogonality and interactions significantly influence polarization and dispersion.

Spin-orbit coupling affects bulk and edge state properties.

Abstract

In this work we present a tight-binding model that allows to describe with a minimal amount of parameters the band structure of exciton-polariton lattices. This model based on $s$ and $p$ non-orthogonal photonic orbitals faithfully reproduces experimental results reported for polariton graphene ribbons. We analyze in particular the influence of the non-orthogonality, the inter-orbitals interaction and the photonic spin-orbit coupling on the polarization and dispersion of bulk bands and edge states.

Figures

Equations

$$|\Psi ⟩=\sum _{i=1}^N{c}_i|{\psi }_i⟩.$$

$$\mathbf{H}\mathbf{c}=\epsilon \mathbf{S}\mathbf{c},$$

$$H_{ij}=⟨{\psi }_i|\hat{H}|{\psi }_j⟩,S_{ij}=⟨{\psi }_i⟩{\psi }_j.$$

$$H_{ij}=\frac{\left({\epsilon }_i+{\epsilon }_j\right)}{2}{S}_{ij}+t_{ij},$$

$$\tilde{\mathbf{H}}\tilde{\mathbf{c}}=\epsilon \tilde{\mathbf{c}},$$

$$\tilde{\mathbf{H}}={\mathbf{S}}^{-\frac{1}{2}}\mathbf{H}{\mathbf{S}}^{-\frac{1}{2}}.$$

$$\hat{H}=\sum _{ij}\left[\frac{\left({\epsilon }_i+{\epsilon }_j\right)}{2}{S}_{ij}+t_{ij}\right]{\hat{b}}_i^{†}{\hat{b}}_j.$$

$${\hat{b}}_i^{†}=\sum _j{S}_{ij}^{-1}{\hat{a}}_j^{†},$$

$$\hat{H}={\hat{H}}_s+{\hat{H}}_p+{\hat{H}}_{sp},$$

$${\hat{H}}_s$$

$${\hat{H}}_p$$

$$+$$

$${\hat{H}}_{sp}$$

$${\hat{\mathbf{b}}}_{ip\sigma }={\hat{b}}_{i{p}_x\sigma }\hat{\mathbf{x}}+{\hat{b}}_{i{p}_y\sigma }\hat{\mathbf{y}}.$$

$${\mathbf{e}}_{ij}^L$$

$${\mathbf{e}}_{ij}^T$$

$$t_{\beta }=t{S}_{\beta }.$$

$$n(r)=\left\{\begin{array}[]{ll}n_{1}\textup{ \, if \, }r\leq R,\\ n_{2}\textup{ \, if \, }r>R,\end{array}\right.$$

$$E_{lm}(\bm{r})=\left\{\begin{array}[]{ll}\frac{1}{\sqrt{N}}J_{l}\left(q_{lm}r\right)\,e^{il\phi}\textup{ \, if \, }r\leq R,\\ \\ \frac{A}{\sqrt{N}}\,K_{l}\left(\kappa_{lm}r\right)\,e^{il\phi}\textup{ \, if \, }r>R.\end{array}\right.$$

$$q_{lm}R\frac{J_{|l|+1}(q_{lm}R)}{J_{|l|}(q_{lm}R)}={\kappa }_{lm}R\frac{K_{|l|+1}({\kappa }_{lm}R)}{K_{|l|}({\kappa }_{lm}R)},$$

$$n_1^2{k}_0^2$$

$$n_2^2{k}_0^2$$

$$q_{lm}R=\sqrt{V^2-{\left(\frac{n_2}{n_1}\right)}^2{\left({\kappa }_{lm}R\right)}^2},$$

$$V=\frac{k_{z}R}{n_1}\sqrt{n_1^2-n_2^2}=\frac{2\pi R}{{\lambda }_{\text{cav}}}\sqrt{n_1^2-n_2^2}.$$

$${\omega }_{lm}=\hslash c\frac{2\pi }{{\lambda }_{\text{cav}}}+\hslash c\frac{{\lambda }_{\text{cav}}}{2\pi }\frac{(q_{lm}R{)}^2}{{\left(n_{1}R\right)}^2}.$$

$${\psi }_s(\mathbf{r})$$

$${\psi }_{p_x}(\mathbf{r})$$

$${\psi }_{p_y}(\mathbf{r})$$

$$S_s$$

$$S_L$$

$$S_T$$

$$S_{sp}$$

$$\mathbf{H}\mathbf{c}=\epsilon \mathbf{c}+\epsilon \tilde{\mathbf{S}}\mathbf{c},$$

$$\mathbf{H}\mathbf{c}=\frac{\epsilon }{1-\lambda \epsilon }\mathbf{c},$$

$${\epsilon }_{\mathbf{k}}=\frac{{\bar{\epsilon }}_{\mathbf{k}}}{1+\lambda {\bar{\epsilon }}_{\mathbf{k}}},$$

Full text

I Introduction

Coupled photonic resonators have appeared in the past few years as an excellent platform to engineer lattice Hamiltonians Nolte and Stefan (2010); Houck et al. (2012); Bellec et al. (2013); Jacqmin et al. (2014). The possibility of controlling the geometry, on-site energy and hopping along with the spectroscopic access to the momentum- and real-space distributions of the wavefunctions are opening new perspectives in the study of elaborate solid-state Hamiltonians in the photonics realm. In addition, the engineering of gain and losses and the presence of Kerr nonlinearities are unveiling genuinely photonic phenomena in lattices, which include lasing in topological states St-Jean et al. (2017); Bahari et al. (2017); Bandres et al. (2018); Zhao et al. (2018); Parto et al. (2018); Klembt et al. (2018), PT-symmetric Weimann et al. (2017) and charge conjugated phases Poli et al. (2015) and the observation of dissipative phase transitions Fitzpatrick et al. (2017); Rodriguez et al. (2017).

Lattices of polariton resonators in semiconductor microcavities provide one of the most versatile platforms to implement this kind of Hamiltonians Schneider et al. (2017). Polaritons are hybrid light-matter quasiparticles that arise from the strong coupling of quantum-well excitons and photons confined in a micron-scale Fabry-Perot cavity. Their excitonic component results in significant polariton interactions and in sensitivity to external magnetic fields. The first feature has allowed the observation of bi-stability Baas et al. (2004), polariton superfluidity Amo et al. (2009) and solitons Amo et al. (2011); Sich et al. (2011) in planar structures, while the second has been used to demonstrate lasing in circularly polarised states Sturm et al. (2015) and in chiral edge states Klembt et al. (2018).

A very convenient way to implement lattices of polariton resonators is by confining their photonic component in fully or partially etched structures. The building block of these lattices is typically a resonator of cylindrical symmetry, in which photons are confined in the three spatial directions. The polariton resonators show confined modes separated by a gap, each of them with a particular geometry: the ground state is formed by cylindrically symmetric $s$-modes, the first excited state is doubly degenerate with $p$-type modes, the next states have $d$-symmetry and so on. Such confined modes have been realized by fully etching the semiconductor structure Bayer et al. (1998); Bajoni et al. (2008), by partial etching of the upper cavity mirror Klembt et al. (2017); Whittaker et al. (2018), by growth interruption and etching of the cavity spacer Kaitouni et al. (2006); Winkler et al. (2015), and in half cavities closed by an external mirror Zhang et al. (2014); Dufferwiel et al. (2014); Besga et al. (2015). By laterally coupling the photonic modes of the resonators, lattices of different geometries have been implemented, including one-dimensional regular Bayer et al. (1999); Tanese et al. (2013); Winkler et al. (2016), Stub Baboux et al. (2016), Su-Schrieffer-Heeger (SSH) lattices St-Jean et al. (2017) and aperiodic lattices Tanese et al. (2014), and two-dimensional honeycomb Jacqmin et al. (2014); Klembt et al. (2018) and Lieb lattices Klembt et al. (2017); Whittaker et al. (2018), showing a wide variety of dispersions and topological features. One of the great assets of this system is the possibility of designing lattices with synthetic strain, which have been recently employed to engineer new types of Dirac cones Milićević et al. (2018) and are promising to engineer artificial gauge fields Rechtsman et al. (2013); Salerno et al. (2015).

The design of polariton lattices and the interpretation of the polariton bands measured in photoluminescence studies have so far largely relied on the mapping to a tight-binding model. In this model, each orbital mode of each cylindrical microresonator is independent from the other orbitals and plays the role of a point-like tight-binding site, all of them with identical on-site energy, coupled to their nearest neighbors. In SSH, Lieb and honeycomb geometries, this kind of tight-binding Hamiltonian presents chiral symmetry and, therefore, the upper and lower bands of eigenvalues should be mirror symmetric with respect to the value of the on-site energy. However, this simple model shows significant deviations from the experimentally observed dispersions, both in 1D and in 2D lattices Jacqmin et al. (2014); Klembt et al. (2018); Winkler et al. (2016); Baboux et al. (2016); Klembt et al. (2017); Whittaker et al. (2018). In particular, in experimental observations, a significant asymmetry between upper and lower bands is systematically observed. An efficient way to fit this band asymmetry is to add a next-nearest neighbor coupling to the tight-binding model. This technique was used, for instance, in the works of Jacqmin, Baboux and coworkers Jacqmin et al. (2014); Baboux et al. (2016).

Despite the apparent success of the fits, the question of the physical relevance of the actual next-nearest neighbor coupling remains, particularly in structures based on complete etching of the semiconductor microcavities, for which the photonic confinement is expected to be very strong within the physical dimensions of the micropillar. Therefore, the observed band asymmetries call for other corrections to the tight binding description. One of them is the coupling between modes of different symmetry belonging to nearest neighbor sites: $s$ and $p$-modes or $p$ and $d$-modes in adjacent micropillars. Simultaneously, the significant spatial overlap between adjacent micropillars in real structures raises questions about the accuracy of the tight binding model, which assumes the limit of weak overlaps. When the overlaps are significant, the original basis made of the individual uncoupled resonators is far from an orthogonal basis, and non-orthogonal corrections need to be added to the original tight-binding Hamiltonian. In models like the honeycomb lattice, these corrections have been shown to result in band asymmetries quite similar to those induced by next-nearest neighbors McKinnon and Choy (1995). Understanding the effects of these corrections is of crucial importance to interpret a number of physical observations within this model.

In this article, we show that, indeed, the experimental dispersion of lattices of polariton micro-pillars can be described with very high accuracy using a realistic tight-binding model that takes into account both the non-orthogonality of the micro-pillar basis and the coupling between $s$- and $p$-bands. We show in this way that direct next-nearest neighbors coupling is not necessary to fully reproduce all experimentally observed phenomenology. To complete our description we take into account the TE-TM splitting characteristic of dielectric microcavities. We compare our model to experimental dispersions obtained in a honeycomb lattice of coupled micro-pillars. Our results should improve significantly polariton tight-binding models.

The rest of the paper is organized as follows: in Section II we introduce the basics of our non-orthogonal tight binding model and a simple variational approach based on low contrast refraction indices that implements an effective modeling of the single pillar photonic modes inside the lattice. We apply our model to the case of a honeycomb lattice in Section III and compare our results with experimental data showing that they agree quite well, even in the case of distorted lattices. Finally, we conclude in Section IV.

II A minimal tight-binding description for cavity-polariton lattices

In this Section we present a simple tight-binding (TB) approach to describe cavity-polariton lattices made out of single cavity micropillars with several polaritonic modes of different symmetries. At the core of the method lies the fact that we will consider the case of weakly coupled cavities where the photonic modes of a single cavity are a good starting point of the calculation. We will explicitly take into account the overlap between photonic modes at nearest neighbors cavities, including those with different symmetries, since this turns out to be very important to describe the experimental data. Our approach is similar in spirit to the one developed in Ref. Kamalakis et al. (2006) for photonic crystals, where the global photonic field was written as a linear combination of the modes corresponding to isolated pillars located at each lattice site. Here, however, we will show that to effectively capture the behavior of the real photonic modes of a pillar due to the spatial overlap with its neighbors it is important not to consider the modes of an isolated pillar in vacuum but those of a pillar surrounded by an effective media.

II.1 Non-orthogonal tight-binding approach

We first summarize the basics of the usual tight-binding (TB) approach involving a non-orthogonal set of localized orbitals in a lattice (see for instance Refs. McKinnon and Choy (1995) and Paxton (2009)). For simplicity, we start by considering the case of a single orbital per site. Generalization to multi-orbital sites is done at the end the section.

For a system with $N$ sites positioned at $\bm{R}_{i}$ with $i=1,...,N$, the single-particle wave-function in the TB approximation is given by the linear combination,

[TABLE]

Here $\ket{\psi_{i}}$ is the orbital state of a single pillar at site $\bm{R}_{i}$ which is assumed to be normalized. The Schrödinger equation $\hat{H}\ket{\Psi}=\varepsilon\ket{\Psi}$ can then be reduced to the following matrix equation

[TABLE]

where we have introduced the notation $\bm{c}=(c_{1},c_{2},\dots,c_{N})^{\mathrm{T}}$. The elements of the Hamiltonian ($\bm{H}$) and overlap ($\bm{S}$) matrices are given by

[TABLE]

In the case of an orthogonal basis, $H_{ij}$ is simply either the on-site energy $\varepsilon_{i}$ of the orbital $\ket{\psi_{i}}$ (for $i=j$) or the so-called hopping matrix element $t_{ij}$ (for $i\neq j$) that comes from the inter-site potential term. Both are usually taken as independent parameters. In the non-orthogonal case however, these terms are mixed and $H_{ij}$ can be parameterized in different (equivalent) ways. We choose the following one

[TABLE]

with $t_{ii}=0$. This symmetric way to represent $H_{ij}$ will allow us later on, to make a simple approximation to the inter-orbital coupling ($t_{ij}$) and to consider the cases where the on-site energy changes from site to site. Notice also that this parametrization explicitly takes into account that a global shift of the site energies translates in a global shift of the bands. On the other hand, since in our case the orbital states are taken to be real functions we have that $S_{ij}=S_{ji}$. Equation (2) can be solved by making the substitution $\bm{c}=\bm{S}^{-\frac{1}{2}}\tilde{\bm{c}}$ to get the more familiar orthogonal eigenvalue problem

[TABLE]

with

[TABLE]

Once solved, $\bm{c}$ can be recovered from $\tilde{\bm{c}}$ by back-substitution. A very convenient way of writing the Hamiltonian is in terms of second quantization operators. Hence, we introduce a set of bosonic creation and annihilation operators, $\hat{b}_{i}^{\dagger}$ and $\hat{b}_{i}$, respectively, which allow us to rewrite the Hamiltonian as

[TABLE]

These operators, in the non-orthogonal case, are not those that create or annihilate a particle in the state $\ket{\psi_{i}}$. In fact, they are a linear combination of the latter. If we denote such operators by $\hat{a}_{i}^{\dagger}$, so that $\ket{\psi_{i}}=\hat{a}_{i}^{\dagger}\ket{0}$ where $\ket{0}$ is the vacuum state, then we have that

[TABLE]

with $\left[\hat{a}_{i}^{\dagger},\hat{a}_{j}\right]=S_{ij}$ and \left[\hat{a}_{i}^{\dagger},\hat{a}_{j}^{\dagger}\right]=\big{[}\hat{a}_{i},\hat{a}_{j}\big{]}=0.

II.2 Multi-orbital model

Let us now consider the more realistic case of a lattice in which each micropillar supports several polariton modes of different symmetry related to its polar angle distribution. We restrict ourselves to the case of $s$ and $p$ orbitals because it will be enough to explain the experimental data described below. Generalization to more orbitals is straightforward by applying the same procedure. Using the non-orthogonal tight-binding approximation described above, and considering only nearest-neighbors (NNs) overlap and hopping terms, the Hamiltonian (7) can be written as

[TABLE]

where the first two terms describe the coupling between the same type of orbitals ($s$ and $p$), and the last one the coupling between different types of orbitals ($s-p$ coupling). In its most general form, and accounting for the two polarization modes of each orbital, these three terms are given by

[TABLE]

This Hamiltonian is a generalization of the models studied in Refs. Jacqmin et al. (2014); Milićević et al. (2015, 2017); Gulevich et al. (2016, 2017); Nalitov et al. (2015); Baboux et al. (2016); Whittaker et al. (2018); Li et al. (2018) to account for the $s$-$p$ inter-orbital coupling and overlap. Here, the operator $\hat{b}_{i\ell\sigma}^{\dagger}$ ($\hat{b}_{i\ell\sigma}$) creates (annihilates) a polariton at site $i$, and in its NNs, according to Eq. (8). The index $\ell$ labels the considered orbital in the single pillar eigenmodes ($s$, $p_{x}$ and $p_{y}$), while $\sigma=\pm$ indicates the polarization of the photon component in the circular polarization basis and $\sigma/\bar{\sigma}$ indicate opposite polarization. The coupling term for the $s$-bands, $t_{s}$ is spatially isotropic, while for the $p_{x}$ and $p_{y}$ orbitals (Eq. 11) it is different depending on the orientation of the orbital with respect to the direction of the link between adjacent micropillars Wu et al. (2007); Jacqmin et al. (2014): $t_{L}$ when orbitals are oriented parallel to the link, and $t_{T}$ for orbitals oriented perpendicular to the link (usually $|t_{L}|\gg|t_{T}|$). To describe this feature we have used a compact vector notation, similar to the one employed in Sala et al. (2015) to represent the operators that act over the $p$ orbitals, namely,

[TABLE]

In this way the expression $(\bm{e}_{ij}^{L/T}\cdot\,\hat{\bm{b}}_{jp\sigma})$ selects the component of $\hat{\bm{b}}_{jp\sigma}$ in the direction specified by the unit vectors

[TABLE]

where $L$ ($T$) indicates whether the unit vector points in the longitudinal (transverse) direction to the link $ij$, whose orientation is given by the angle $\varphi_{ij}$ (see Fig. 1(a)). Therefore, these vectors select the projection of the $p_{x/y}$ orbitals parallel (perpendicular) to the lattice bond. In addition to the overlap and hopping terms that conserve the polarization we also include spin-orbit coupling (SOC) terms that flip it Nalitov et al. (2015). For the sake of simplicity we have assumed that the SOC strength is proportional to the corresponding direct coupling, and we model it by the adimensional parameter $\Delta$. This SOC term arises from the fact that the coupling between pillars depends on whether the polariton polarization is parallel or perpendicular to the bond, owing to the fact that the two polarization modes experience different tunnel barriers in the presence of TE-TM splitting Sala et al. (2015).

Finally, Eq. 12 describes the coupling between $s$ and $p$ orbitals in adjacent sites, given by the coupling strength $t_{sp}$ (see Fig. 1(a)). The summation $\braket{i,j}$ in Eqs. 10-12 runs over each site of the lattice and its NNs which is very appropriate for most lattices (see for instance Fig. 2(a)). An estimate of the magnitude of the second nearest neighbors hopping and overlap terms using the effective model for the photonic modes presented in the next section shows that they are roughly a factor $10^{-3}$ smaller than the ones corresponding to NNs (see also Appendix A).

At this point, to describe the experimental data, all the parameters of the Hamiltonian could be taken as free fitting parameters. This has been the most common approach used in the literature so far. While this gives reasonable results in many cases, it would be desirable to reduce the number of parameters, even though the approach might result less flexible, to gain a better physical insight and gain some predictability power on the design of different polaritonic lattices. In order to do so, we will assume in what follows that all hopping elements $t_{\beta}$, where $\beta\in\left\{s,L,T,sp\right\}$, can be written as

[TABLE]

This is a reasonable assumption if one notices that for a system of exciton-polariton microcavity pillars spatially overlapped (see Fig. 2) the inter-site potential entering the definition of $t_{ij}$ (see Eq. (4)) might be considered as being constant within the microstructure where the two (photonic) modes mainly overlap. It is important for this to be valid to have used the symmetrized parametrization shown before so that $t_{ij}=t_{ji}$. Equation (16) immediately eliminates the need to distinguish between $t_{s}$, $t_{L}$, $t_{T}$ or $t_{sp}$ as they are all determined by the same parameter $t$ and the corresponding overlap matrix element (which itself is fixed by the choice of a single parameter $V$ that we will introduce in the following section).

II.3 Single pillar mode: a simplifying approximation

Following the spirit of reducing the number of free parameters in the model, we present here a simple way to approximate the photonic eigenmodes of a single micropillar in order to calculate the overlap integrals that appear in the definition of Hamiltonian (9). To this end, and for reasons that will become clear below, we consider a cylindrical microcavity defined as an infinite long ($z$-axis) circular dielectric waveguide with a step refraction index profile Gérard et al. (1996),

[TABLE]

where $n_{1}$ is the refractive index of a core of radius $R$ and $n_{2}$ is the refractive index of the surrounding material. In doing this, we have ignored the $3$D nature of the problem and treat it as effectively $2$D. This is a valid assumption as far as the wavelength of the confined modes are much larger than the micropillar resonant wavelength.

The calculation of the electromagnetic modes of this type of waveguide is a well-known problem in the literature (see for instance Yariv and Yeh (2007); Kapany and Burke (1972); Gérard et al. (1996)). The exact solutions for the propagating modes are in general a mixture of transverse electric (TE) and transverse magnetic (TM) waves. Finding these hybrid modes —usually refereed to as EH and HE modes depending on which of the electric E or magnetic H of the photon field is non-zero along the propagation direction ($z$)— is somehow involved as it requires solving the wave equations in cylindrical coordinates with the different components of the fields being coupled. However, a good approximation for both the fields and the mode equation that determines the mode frequencies can be obtained if we assume that the core refractive index ($n_{1}$) is only slightly higher than that of the surrounding medium ($n_{2}$). This is an approximation that is often used for describing optical fibers. At first glance, this approach may seem very far from the real scenario of an isolated pillar surrounded by vacuum since $n_{1}\approx 3.5$ and $n_{2}=n_{0}\approx 1$. However, when dealing with lattices like the one depicted in Fig. 2(a), where each micropillar is spatially overlapped with their NN’s, it is reasonable to expect that those neighboring pillars will provide a substantially different environment to the central one, and modify the photon confinement as compared to an isolated micropillar. We propose that the presence of the neighboring pillars, which in general tends to delocalize the single pillar mode, can be effectively described by considering an effective refractive index for the surrounding medium. In this sense, our approach is variational. Because all the micropillars are made of the same material we expect $n_{2}\lesssim n_{1}$. Its precise value, of course, may depend on the lattice geometry: the presence of a higher/lower number of nearest neighbors will result in greater/smaller delocalization of the modes in a considered micropillar.

It is important to note at this point that even in this approach the photonic modes are well confined within the micropillar and continue to be a good starting point to build a TB model (see Appendix A). Therefore, by assuming that $n_{1}-n_{2}\ll 1$ we simplify the equations for matching the field components at the $r=R$ interface (see Ref. Yariv and Yeh (2007) for details).

In this limit, the modes become linearly polarized (say, along the $\hat{x}$ and $\hat{y}$ directions), the two polarizations being degenerated, and the electric field amplitude is given by

[TABLE]

where $r$ and $\phi$ are the corresponding polar coordinates of $\bm{r}$, and $J_{l}(x)$ and $K_{l}(x)$ are the Bessel functions of the first kind and modified of the second kind, respectively. Here we have ignored the $z$-dependence of the fields (plane wave) since it is irrelevant for our purpose, $A=\frac{J_{l}\left(q_{lm}R\right)}{K_{l}\left(\kappa_{lm}R\right)}$, whereas $N$ is determined by the normalization condition. The transverse wavenumbers inside and outside of the waveguide, $\left(q_{lm},\kappa_{lm}\right)$, are determined by the following equation

[TABLE]

where the orbital index $l=0,\pm 1,\pm 2,...$ and the subscript $m$ indicates the $m$-th root of this transcendent equation. Equation (18) can be solved numerically if one takes into account that $q_{lm}$ and $\kappa_{lm}$ are related by

[TABLE]

with $k_{0}=\omega_{lm}/c$ and where $\omega_{lm}$ is the frequency of the mode $lm$. It is clear then that Eqs (19) and (20) can be rewritten as

[TABLE]

where

[TABLE]

Here we have used the usual experimental condition for isolated micropillars, $k_{z}=\frac{2\pi}{\lambda_{\text{cav}}}n_{1}$, where $\lambda_{\text{cav}}$ is the micropillar resonant wavelength (see for instance Refs. Bayer et al. (1998); Bajoni et al. (2008)).

At this point, for simplicity, we can approximate $\left(\frac{n_{2}}{n_{1}}\right)^{2}\approx 1$ in Eq. (21), and hence the solutions of Eq. (18) can be parameterized in a very convenient way with a single parameter $V$ that effectively captures how well the electromagnetic field is confined within the micropillar. Note that in this limit, the classical electromagnetic problem is analogous to the quantum problem of a particle of mass $m$ in a finite circular potential well of radius $R$ and magnitude $U=\hbar^{2}V^{2}/2mR^{2}$, if we interpret the electric field as the wave function amplitude. Note also that the energies of the modes in the electromagnetic problem can be calculated as $\omega_{lm}=\frac{\hbar c}{n_{1}}\sqrt{k_{z}^{2}+q_{lm}^{2}}$, which in the limit $k_{z}\gg q_{lm}$ can be rewritten as

[TABLE]

Finally, we can define the $s$ and $p$ modes (for each polarization, linear or circular) as

[TABLE]

The ($2$D) overlap integrals involved in our model can be calculated as follows

[TABLE]

We emphasize that in this approximation all the overlaps are determined by $V$, which then plays the role of a variational parameter that effectively describes the delocalization of the photonic modes due to the penetration in the adjacent overlapping micropillars as compared with those of an isolated pillar. Note that in the limit $n_{1}-n_{2}\ll 1$ we are considering, the confined modes are polarization degenerate, and Eqs. (27) are polarization independent. Polarization effects related to the spin-orbit coupling (TE-TM splitting) is phenomenologically incorporated in our model via the $\Delta$ terms in Eqs. (10), (11), (12).

III The honeycomb lattice

The extraordinary transport and topological properties of graphene have stimulated a number of experimental and theoretical studies of the polariton honeycomb lattice Kusudo et al. (2013); Jacqmin et al. (2014); Milićević et al. (2015, 2017, 2018, 2018); Bleu et al. (2016, 2017); Nalitov et al. (2015, 2015); Ozawa et al. (2017); Solnyshkov et al. (2016); Kartashov and Skryabin (2017); Solnyshkov et al. (2018); Klembt et al. (2018). Here we analyze the bulk band structure and the edge states spectrum based on the complete tight-binding model presented in the previous section, highlighting the role of its different physical ingredients. We compare our numerical results with experimental data and show that they provide a very good description of the band structure. We also point out some specific signatures of the spectrum related to the photon polarization that might be relevant for future experiments. Finally, we reproduce recently published experimental results on the emergence of tilted Dirac cones in polariton graphene lattices under strain Milićević et al. (2018), showing that we capture correctly the dependence of parameters with distance. This gives our model certain predictive capability that could be useful to engineer different effects on artificial microcavity-polariton lattices.

III.1 Bulk bands

To achieve a better understanding of the influence of the different terms in Hamiltonian (9) let us consider first the bulk bands along a specific direction of high symmetry of the underlying lattice. For that, we define the lattice vectors $\bm{a}_{1}=\sqrt{3}\,a\,\hat{\bm{x}}$, $\bm{a}_{2}=\frac{\sqrt{3}}{2}\,a\,(\hat{\bm{x}}+\sqrt{3}\,\hat{\bm{y}})$ and the relative position of the basis sites $A$ and $B$, $\bm{\delta}_{AB}=\frac{a}{2}\,(\sqrt{3}\,\hat{\bm{x}}+\hat{\bm{y}})$ (see Fig. 1(b)). Here $a$ is the distance between two NNs pillars. The calculated spectrum for a lattice of micropillars of a diameter $D=2R=3~{}\mu$m and a center-to-center distance $a=2.4~{}\mu$m, along $k_{y}=0$, is show in the Fig. 3. In the different panels of the figure, we analyze the contribution of the different terms of the model separately, that is, considering only one of them at a time. In each panel, the black dashed line represents the bands in the absence of $s$-$p$ coupling ($t_{sp}=0$), non-orthogonality ($S_{\beta}=0$), and SOC ($\Delta=0$). In this case each band is clearly particle-hole symmetric and the upper and lowermost $p$-bands present a very small dispersion –this bands would be completely flat for $t_{T}=0$. The red solid lines in each panel of Fig. 3(a)-(c) include each contribution separately: (a) only $s$-$p$ coupling ($\Delta=0$, $S_{\beta}=0$), (b) only non-orthogonality ($\Delta=0$, $S_{sp}=0$ and $t_{sp}=0$), and (c) only SOC ($S_{\beta}=0$, $t_{sp}=0$) –see the figure caption for the value of all the parameters.

As clearly seen in the figure, each term leads to a different effect on the bands. The inter-orbital coupling (Fig. 3(a)) plays a very important role on the deformation of the bands as it tends to join them, stretching the top of the $s$-band and the bottom of the $p$-band, and making them to acquire a V-like shape in the neighborhood of the $\Gamma$ point. Notice also that the uppermost and lowermost $p$-bands are not affected by this coupling. This is to be expected as those bands involved the $p$-orbitals that are perpendicular to the bonds and hence they do not couple to the $s$-bands. On the contrary, one of the main effects of the non-orthogonality between orbitals in different sites (Fig. 3(b)) is to produce a clear asymmetry between those quasi-flat bands, making the uppermost wider and the lowermost narrower. This point is quite relevant as it will allow us to reproduce the experimental data quite well without the need to include an energy-dependent hopping (as previously done in  Jacqmin et al. (2014)). We can highlight that the non-orthogonality induces the opposite asymmetry of the bands as compared to the $s$-$p$ coupling. This effect is significant for the $p$-bands, while it remains negligible for the $s$-bands. Indeed, from the overlap integrals Eq. (27) we estimate the overlap between $p$ orbitals in adjacent micropillars to be of the order of $15$% for typical lattices, while it is only $4$% for the $s$-bands. The effect of the SOC (Fig. 3(c)) is simply to split the bands, as expected, leading to the appearance of a polarization (spin) texture. Finally, Fig. 3(d) shows the bands including all terms.

III.2 A comment on the effect of non-orthogonality

We would like to emphasize some important aspect of our model for cavity polariton lattices. On the one hand, we neglect the second nearest neighbors ($2$NN) hopping terms. This is so because in most of the lattices there is no overlap between $2$NNs micropillars and hence the coupling of the photonic modes goes through the extremely weak evanescent field in vacuum, out of the micropillars. On the other hand, we do include the non-orthogonality between orbitals located in adjacent sites. Its effect on the asymmetry of the bulk bands (Fig. 3(b)) could be qualitatively reproduced in a phenomenological way by including an effective $2$NN hopping term (see, for instance, Ref. Jacqmin et al. (2014)). However, $2$NN hopping would have a very different effect on the flat band states localized in zigzag and bearded edges: the $2$NNs hopping would destroy the flatness of the edge states band, while the non-orthogonality preserves it. This can be easily understood as follows. We re-write Eq. (2) as

[TABLE]

where $\tilde{\bm{S}}$ is the non-diagonal part of $\bm{S}$. Now, we will assume, for the sake of simplicity and to make the argument clear, that we are only considering a set of equivalent orbitals so that all energy sites can be taken to be equal to zero without any loss of generality. In addition, we continue to assume that the hopping terms are proportional to the overlap between NNs sites. Under these assumptions $\tilde{\bm{S}}=\lambda\bm{H}$, with $\lambda<0$ and $|\lambda|\ll 1$. Therefore

[TABLE]

and so the band structure is given by

[TABLE]

with $\bar{\varepsilon}_{\bm{k}}$ the energy dispersion of the orthogonal case ($\lambda=0$). It is then clear that: i) flat bands remain flat when the non-orthogonality is included. In particular, the ones at $\bar{\varepsilon}_{\bm{k}}=0$ do not move when non-orthogonality is included. Therefore, flat band edge states characteristics of zigzag and armchair edges are not affected by non orthogonal effects; ii) upper (lower) bands that correspond to $\bar{\varepsilon}_{\bm{k}}>0$ ($\bar{\varepsilon}_{\bm{k}}<0$) gets broader (narrower) as the factor $(1+\lambda\,\bar{\varepsilon}_{\bm{k}})^{-1}$ is bigger (smaller) than $1$. Of course, in addition to this there is also some deformation of the original band. We emphasize once again that this effect is opposite to the one induced by the $s$-$p$ coupling (see Fig. 3(b)). Therefore, the presence in the experimental data of a clear asymmetry between the lowest and uppermost $p$-bands is an evidence of the importance of the non-orthogonality while the asymmetry of the inner middle bands is an indication of the relevance of the $s$-$p$ coupling.

When more orbitals are involved and coupled, as in our numerics, the above picture gives only a qualitative description of the non-orthogonality effect. Furthermore, we expect this picture to hold even if we extend our model and include (slightly) different proportionality constants between the different hopping terms and the corresponding overlap matrix elements.

III.3 Polariton graphene ribbon and comparison with the experiment

We now calculate the band structure for a polariton graphene ribbon (PGR) with zigzag edges (see Fig. 1(b)), taken to be infinite along the $x$ direction (hence we can use Bloch theorem) and containing $30$ unit cells along the transverse direction (defined by $\bm{a}_{2}$). We use the same TB parameters as above.

To analyze the polarization properties of the bulk (center of the ribbon) and edge states, and compare with the data obtained from photoluminescence experiments, we define a quantity that describes the probability of detecting a state (via the emission of a photon) located at the site $\mathbf{R}_{n\alpha}=n\mathbf{a}_{2}+\bm{\delta}_{\alpha}$ (where $\bm{\delta}_{\alpha}$ gives the position of the non-equivalent micropillar $\alpha\in\left\{A,B\right\}$ in the unit cell) with polarization $\sigma$ and quasi-momentum $\bm{k}=k_{x}\,\hat{\bm{x}}+k_{y}\,\hat{\bm{y}}$. Namely,

[TABLE]

where $G_{n^{\prime}\alpha^{\prime}n\alpha}=\exp\left[\frac{-\left[\left(\mathbf{R}_{n^{\prime}\alpha^{\prime}}-\mathbf{R}_{n\alpha}\right).\hat{\bm{y}}\right]^{2}}{2\Sigma^{2}}\right]$ is a Gaussian weight function focused on the site $\mathbf{R}_{n\alpha}$ with standard deviation $\Sigma$, which allows to replicate the spatial dependence of the emitted light in experiments. Indeed, in photoluminescence experiments of photonic graphene micropillars, the escape of photons out of the microcavity results in a Gaussian-like distribution of the emission centered at the position of the pump spot Jacqmin et al. (2014). $\bm{S}_{k_{x}}$ and $\bm{c}_{k_{x}}$ are the ribbon overlap matrix and eigenvectors, respectively (see Appendix C), and $C$ is determined for each $k_{x}$ by normalization condition.

Figures 4(a) and 4(c) show the calculated emission spectrum ($p_{n\alpha\bm{k}}^{\sigma}$) for the case of excitation of a bulk site located at the center of the ribbon, $\left(n=15,\alpha=A\right)$— along the path $\bm{k}=k_{x}\,\hat{\bm{x}}$, and setting $\Sigma=3a$. We included in the simulations an artificial broadening parameter ($\gamma=0.07$ meV) with the sole purpose of reproducing the experimental linewidth. Figures 4(b) and (d), show the corresponding experimental results for the photon emission polarized parallel or perpendicular to the ribbon’s edge, respectively. The experimental conditions are those of Ref. Milićević et al., 2017: an AlGaAs-based microcavity, with 28(40) Bragg pairs in the upper(lower) mirror, with 12 GaAs quantum wells etched into a lattice of micropillars of a diameter of $3~{}\mu$m and a center-to-center distance $a$ of $2.4~{}\mu$m; a non-resonant laser at 740 nm excites the bulk of the lattice in a spot of $4~{}\mu$m in diameter; the light emitted from the polariton bands is collected as a function of the linear polarization direction, emitted angle (in-plane momentum) and wavelength. In Fig. 4 the bands are measured along the $k_{x}$ direction for an angle of emission in the $y$ direction corresponding to $k_{y}=4\pi/3a$, passing through the centre of the second Brillouin zone. The reason for selecting the emission through the second Brillouin zone is to avoid destructive interference effects characteristics of bipartite lattices that prevent clear observation of the bands at the center of the first Brillouin zone Shirley et al. (1995); Jacqmin et al. (2014). In order to make a detailed comparison between experiment and theory, in Fig. 4(b) and (d), on top of the experimental data we plot the numerical results (dotted lines) corresponding to the bulk system polarized in the $x$ or $y$ direction as appropriate (same data as in Fig. 3(d)). Many details of the experimental spectrum are clearly captured by our simple model. We emphasize that our approach for this (bulk) case involves only a few fitting parameters: $\Delta E$, $t$, $V$ and $\Delta$.

The corresponding figures for the case of photoluminescence from a site at the edge of the ribbon ($n=1,\alpha=A$) are presented in Fig. 5. Here the spectrum is slightly more complex as several new features appear so a more detailed analysis is needed. Figure 5(a) shows the complete ribbon band spectrum calculated with our TB model. Red and blue lines correspond to the parallel ($\parallel$) and perpendicular ($\perp$) polarization with respect to the ribbon edge, respectively. Localized edge states (indicated by the arrows) appear both at the $s$ and the $p$ bands. In the latter case, there are two types of edge states, as discussed in Ref.  Milićević et al. (2017): (i) the usual edge states (open arrows), similar to the ones on the $s$ bands, that appear near the Dirac cones and are usually flat—here the dispersion observed for one of the polarizations is mainly due to the difference on the onsite energy of the edge site, see discussion below—, and (ii) the ones with a non-trivial dispersion (solid arrows). In the latter case the splitting is caused mainly by the inclusion of the SOC coupling, being the bands rather polarized.

Figures 5(b) and 5(d) show the corresponding calculated spectrum ($p_{1Ak_{x}}^{\sigma}$) while Figs. 5(c) and 5(e) shows the measured bands along the $k_{x}$ direction for an angle of emission in the $y$ direction corresponding to $k_{y}=4\pi/3a$. A careful analysis of the latter shows that : (i) there is a clear difference between both polarizations in the case of the ‘flat’ edge states (the Dirac cone edge states, highlighted with open arrows in Fig. 5(a)): polarization parallel to the edge (Fig. 5(c)) shows a flat edge state (as naively expected) while for the perpendicular polarization (Fig. 5(e)) it is dispersive —we emphasize here that this effect cannot be accounted for by the inclusion of a second NNs hopping as the later is negligible; (ii) the SOC induced splitting of the ‘dispersive’ edge states (those highlighted with solid arrows in Fig. 5(a)) is a bit stronger for the lower edge bands as compared with the upper edge bands.

We have found that these features can be accounted for in our model by modifying the site energy of the surface (edge) pillars as compared with those of the bulk orbitals. In particular, only the energy of the $p_{y}$ orbital needs to be modified, being different for each polarization. Hence, all the results shown in Fig. 5 include such a change, which is given by $\delta\varepsilon^{\prime}_{p_{y}\parallel}=0.4$ meV and $\delta\varepsilon^{\prime}_{p_{y}\perp}=0.6$ meV. A possible origin for this energy shift might be the combination of the presence of excitonic stress at the edge pillars of the lattice and different confinement of photonic modes when the number of nearest neighbors is reduced with respect to the bulk pillars (micropillars in the zigzag edge have two NNs while in the bulk all pillars have three NNs).

III.4 Strain induced merging of Dirac cones

One of our goals in developing this tight binding model is not only to account correctly for all the different couplings and non-orthogonality effects, and establish their relative importance, but also to be able to predict the band structure up to some fine details as this would be very useful in the design of future experiments. To show this potentiality, we have analyze the case of a distorted honeycomb lattice as the one used in Ref. Milićević et al. (2018), where the length of the bond perpendicular to the zigzag edge was changed to be $a^{\prime}\geq a$. Figure 6 shows experimental luminescence measured at the center of the lattices under similar conditions as Fig. 5, at an exciton photon detuning at the bottom of the $p$ bands of $-10$ meV. Here the micropillars are $2.75~{}\mu$m in diameter and a center-to-center distance of the undistorted bonds of $a=2.4~{}\mu$m, with the strained bond being $a^{\prime}=2.4~{}\mu$m (undistorted lattice, (a)-(b)), $2.6~{}\mu$m (c)-(d), $2.7~{}\mu$m (e)-(f) and $2.72~{}\mu$m (g)-(h). The left/right column corresponds to the emission linearly polarized parallel/perpendicular to the edge. Solid white lines show the calculated dispersions. All the parameters of our model (except for $\Delta E$, which was slightly corrected for each distance) where changed only through their dependence with the bond distance. Quite notably, in agreement with the experimental data, we find that the merging of the Dirac cones occurs for $a^{\prime}=2.7$ $\mu$m for the parallel polarization—for the perpendicular polarization we observe a similar behavior but for a larger distortion due to the presence of the SOC. The latter is a signal that the magnitudes of the overlap and SOC are correct.

IV Final remarks

We have presented a relatively simple tight-binding model to describe generic cavity polariton lattices including the most relevant physical ingredients. Namely, the coupling between single pillar modes of different symmetry ($s$-$p$ coupling) and the non-orthogonality between different sites. A careful analysis and comparison with the experimental data allowed us to identify the most prominent features each contribution introduces and, although they change the band structure with similar magnitudes, it turns out that the $s$-$p$ coupling leads to the most distinguishable effects. This coupling substantially reshapes the bands, particularly the $s$-band, a feature so far neglected in experimental and theoretical polariton studies.

The non-orthogonality plays an important role in the $p$-bands, resulting in an asymmetry in the dispersion of the uppermost and lowermost $p$-bands. Our estimates show that the $p$ orbitals have an overlap between adjacent micropillars of the order of $15$% –for typical lattices–, while for the $s$ ones it is only $4$% and, hence, non-orthogonal effects can be safely ignored for the $s$-bands.

In concordance with this, it is important to emphasize that second NNs hopping is negligible as there is essentially no overlap between second nearest micropillars, and the evanescent field out of the etched micropillars decreases extremely fast. Note that this might not be the case in polariton lattices fabricated with other techniques. In particular, lattices fabricated by partial etching of the structure (upper mirror) Klembt et al. (2017); Whittaker et al. (2018), metallic deposition on the surface, or intracavity mesa techniques  Kaitouni et al. (2006); Winkler et al. (2015), might present deeper evanescent fields and may result in significant second NNs couplings. We stress, however, that second NN hopping and non-orthogonality act very differently on the flat band states localized in zigzag and bearded edges: while the former destroys the flatness of the edge states band, the latter preserves it.

While here we restricted ourselves to consider only the $s$ and $p$ modes, it is rather natural to ask whether higher energy modes should also be included. Calculations show that by adding the $d$ modes similar results are obtained. However, comparisons with experiments reveal that, although qualitatively the structure of the $d$ bands is well captured by the model, it overestimates its bandwidth and its coupling with the lower modes. That is, the measurements show a greater confinement for the $d$ bands than the expected for the model. We argue that this may be due to the proximity of $d$ bands to the exciton-energy, which in the experiments shown here amounts to $-10$ meV for the $p$ bands and about to 0 for the $d$ bands, much closer to the exciton resonance. Therefore, the excitonic contribution to the $d$ polariton states is greater than for the lower modes and the photonic component smaller, then reducing the hopping between pillars and correspondingly the bandwidth. Preliminary calculations using a different parameter $V^{\prime}$ for $d$ bands show slightly better description of the experiment. Yet, that comes at the price of increasing the number of parameters of our simplified parameterization with very few fitting parameters and it does not result in a relevant improvement for the $s$ and $p$ bands.

The theoretical results here presented provide accurate guidelines to describe the band structure of lattices of polariton micropillars, and explain the break up of the particle-hole symmetry observed experimentally and assigned, up to now, to second nearest neighbors effects.

Acknowledgements.

We acknowledge financial support from ANPCyT (grants PICTs 2013-1045 and 2016-0791), from CONICET (grant PIP 11220150100506), from SeCyT-UNCuyo (grant 06/C526), the ERC grant Honeypol, the Quantera grant Inerpol, the FETFLAG grant PhoQus, the French National Research Agency (ANR) project Quantum Fluids of Light (ANR-16-CE30-0021), the Labex CEMPI (ANR-11-LABX-0007) and NanoSaclay (ICQOQS, Grant No. ANR-10-LABX-0035), the French RENATECH network, the CPER Photonics for Society P4S, the I-Site ULNE via the project NONTOP and the Métropole Européenne de Lille via the project TFlight. GU acknowledges support from the ICTP associateship program and thanks the Simons Foundation.

Appendix A Wave functions profiles and penetration length

The approach presented in Sec. II.3 might be, at first glance, somehow anti-intuitive about how well confined is the wave function within the micropillar. To try to clarify this point we show in Fig. 7(a) and 7(b) the profiles of the wave functions $\psi_{s}(x,0)$ and $\psi_{px}(x,0)$, respectively. In each case we plot the wave function centered at $x=0$, $x=a$ (NNs distance for the honeycomb lattice) and $x=\sqrt{3}a$ (NNNs distance for the honeycomb lattice). We use the similar parameters as in the main text, $a=2.40$ $\mu$m, $D=2R=3.00~{}\mu$m and $V=11.0$. Using these values, and approximate experimental values for the micropillar refraction index ($n_{1}\approx 3.5$) and the resonance wavelength of the cavity ($\lambda_{\text{cav}}\approx 740$ nm) we calculate from Eq. (22) the value of the effective refractive index of the external medium $n_{2}\approx 3.4733$. As we see in the figures, although $n_{2}$ has a value very close to $n_{1}$, the wave functions are still well confined, so the values of the NNs overlap integrals are small $S_{s}\approx 0.04$, $S_{L}\approx-0.15$, $S_{T}\approx 0.02$ and $S_{sp}\approx-0.08$ while for NNNs these integrals are negligible $S_{s}\approx 3\times 10^{-5}$, $S_{L}\approx-2\times 10^{-4}$, $S_{T}\approx 1\times 10^{-5}$ and $S_{sp}\approx-9\times 10^{-5}$. In addition, we show the penetration length for the $s$ and $p$ modes (Fig. 7(c)) defined as $1/\kappa_{s}$ and $1/\kappa_{p}$, respectively, for $n_{2}\in[3.400,3.498]$—here $\kappa_{s}=\kappa_{01}$ and $\kappa_{p}=\kappa_{11}$, see Eqs. (24) and (25). Notably, the penetration lengths are very small even for values of $n_{2}$ very close to $n_{1}$. In the case of the $s$ mode $\kappa_{s}^{-1}\propto e^{\frac{2}{V^{2}}}$ when $V\rightarrow 0$ ($n_{2}\rightarrow n_{1}$). This non-analytic behavior is rather particular for a confinement potential in $2$D. For the $p$ modes the behavior is slightly more complex and $\kappa_{p}^{-1}\rightarrow\infty$ when $V\rightarrow z_{01}$ ($n_{2}\rightarrow n_{c}=3.49831$)—here $z_{01}$ is the first zero of the $J_{0}$ Bessel function. The fact that the penetration length diverges for a finite value of the potential means that for $n_{2}>n_{c}$ the $p$ modes are not confined in the micropillar.

Appendix B Two couple pillars: effective media approximation

To test the range of validity of our approach we consider here the case of two coupled pillars, the so-called polaritonic molecule (PM) Michaelis de Vasconcellos et al. (2011), as those shown in Fig. 8(a). We notice that this is not the best scenario for our variational approximation to the single micropillar photonic mode as the effective index $n_{2}$ used to represent the surrounding pillars is somehow harder to justify in this configuration. Yet, we will show that even in this case it provides a very satisfying phenomenological description of the energy separation between the first two modes of the PM, the bonding (B) and anti-bonding (AB) modes.

We use the Hamiltonian (9), without considering the effect of SOC ($\Delta=0$), to fit the experimental results Michaelis de Vasconcellos et al. (2011) for the splitting ($E_{AB}-E_{B}$) as a function of diameter ($D=2R$) of the micropillars. This is done in Fig. 8(a) keeping the parameter $G=a/D-1$ constant. Here $a$ is the separation between the centers of two micropillars. This corresponds to maintain the normalized overlap between the pillars unchanged. Note that for $G<0$, $G=0$, $G>0$, the micropillars overlap, are tangent, and are separated, respectively. Note that to fit correctly the experiment we have to consider the $R$ dependence of the parameters $V$ and $\Delta E$. Namely, from Eq. (22) we have

[TABLE]

and from Eq. (23)

[TABLE]

where $v$ and $\delta E$ where taken as adjusting parameters. The results for the model using $t=-5.2$ meV, $v=2.1\frac{1}{\mu\text{m}}$ and $\Delta E=40.5\frac{\text{meV}}{\mu\text{m}^{2}}$ and the experiment value $G=-0.20$ are shown in Fig. 8(a). The agreement with the experimental data is very good. Using the same value of the parameters we show in Fig. 8(b) the comparison with the experimental splitting ($E_{AB}-E_{B}$) as a function of the distance between the centers of the micropillars , keeping $D$ constant (that is, modifying the value of $G$), for $D=2.5\mu$m and $D=3.0\mu$m. In this case we can see a good agreement between model and experiment for values of $G$ between [math] and $-0.25$. The discrepancies observed for $G>0$ are expected since for this condition the pillars are separated and the model loses validity. In the other case, for $G<-0.25$, the discrepancies can be understood by noting that the approximation given by Eq. (16) overestimates the value of the real hopping integral.

Appendix C Ribbon Hamiltonian

When considering an infinite long ribbon along the $x$ direction, with a finite width ($y$ direction), it is better decompose the eigenstates in plane waves along the ribbon’s direction and define a crystal momentum $k_{x}$. In this case the eigenstates of the system can be written as

[TABLE]

where the Bloch wave functions are given by

[TABLE]

Here $m$ is the index that lists the transverse layers that make up the ribbon and $n\alpha$ is a composite index that labels the intra layer elements, so that the position of each micropillar is given by $\left(m\,\bm{a}_{1}+n\,\bm{a}_{2}\right)+\bm{\delta}_{\alpha}$ where $\bm{a}_{1}=a_{1x}\hat{\bm{x}}$ and $\bm{a}_{2}$ are the primitive vectors, while $\bm{\delta}_{\alpha}$ gives the position of the non-equivalent micropillar $\alpha\in\left\{A,B\right\}$ in the unit cell. The other two indices, $\ell$ and $\sigma$, refer to the orbital and polarization degrees of freedom, respectively.

As mentioned in Sec. II.1, the problem is then reduced to solving the matrix equation

[TABLE]

where in this case the ribbon Hamiltonian and overlap matrix can be written as

[TABLE]

The matrix elements of the layer matrices defined above are given by

[TABLE]