Robust quantum valley Hall effect for vortices in an interacting bosonic quantum fluid

Topologically protected pseudospin transport, analogous to the quantum spin Hall effect, cannot be strictly implemented for photons and in general bosons because of the lack of symmetry-protected pseudospins. Here we show that the required protection can be provided by the real-space topological excitation of an interacting quantum fluid: a quantum vortex. We consider a Bose-Einstein condensate at the Γ point of the Brillouin zone of a quantum valley Hall system based on two staggered honeycomb lattices. We demonstrate the existence of a coupling between the vortex winding and the valley of the bulk Bloch band. This leads to chiral vortex propagation on each side of the zigzag interface between two regions of inverted staggering. The topological protection provided by the vortex winding prevents valley pseudospin mixing and resonant backscattering, allowing a truly topologically protected valley pseudospin transport.

Since the eighties, the concept of topology has been applied to reciprocal space. The topology of Landau levels [3][4][5] and more generally of Bloch bands [6] has been shown to determine the spectacular properties of topological insulators. In this case, the single-particle energy bands of the system are described by topological invariants [5] (such as the Chern number). The field expanded even further with the discovery of the quantum spin Hall effect and of the associated class of Z 2 topological insulators [7,8]. Indeed, if one considers spinor particles in a lattice (electrons for instance), the Chern number computed using only one spin component is not a topological invariant. On the other hand, the difference between two spin Chern numbers is a Z 2 topological invariant for a Hamiltonian verifying Time-Reversal Symmetry (TRS) [9]. In that case the bulk-boundary correspondence applies and guarantees on the interface with a trivial insulator the presence of a pair of counter-propagating spinpolarized states, which because of TRS do not couple the one to the other. This triumph of topology was followed by the attempts to extend the concept of Z 2 topological insulators to other types of two-level systems which can be mapped to a pseudospin representing either an internal degree of freedom (the polarization of a photon) or an external one (angular momentum states, valleys of a honeycomb lattice [10], etc). However, for photons, TRS acts differently from fermions [11] and rigorously, there is no symmetryprotected Z 2 photonic topological insulator. This can be clearly visualized by explicitly considering the pho-tonic spin-orbit coupling due to the energy splitting between TE and TM modes [12][13][14]. It respects TRS, but it has a double winding number which couples counterpropagating spin-polarized photonic modes. The realization of a Z 2 topological insulator analog for light therefore requires to fabricate a structure where the TE-TM splitting is weak, which is possible but very demanding [15,16]. Other degrees of freedom, like the angular momentum of photons in lattices of ring cavities have been considered [17] with the formal problem that no specific symmetry protects this pseudospin which is affected by disorder. Finally, the quantum valley Hall (QVH) effect in staggered honeycomb lattices uses the valley pseudospin [10,18]. It has been evidenced experimentally in electronic systems [19] and recently considered in a large series of works in topological photonics [20][21][22][23][24][25]. Here, the mechanism of dissipation is inter-valley scattering [26]. Even if it is argued to be weak in electronic systems and to be zero for certain types of defects respecting the lattice symmetry in photonics, it formally leads to the Anderson localization of the 1D edge states.
The topology of the quantum fluid in real space and of the band in the reciprocal space have already been fruitfully combined in topological superconductors and superfluids [27][28][29]. The collective excitations of the fluid, described by the Bogoliubov-de Gennes equation, are split off by the superconducting gap, which can become topologically non-trivial for specific shapes of the pairing, creating topological edge states. A vortex, whose core remains in the normal phase, necessarily contains such edge states, which can be Majorana fermions [30] protected by the particle-hole symmetry. Many other solitonic [31][32][33][34][35][36][37] and vortex [38] solutions were found in non-trivial topologies, but the chiral behavior has been mostly discussed for weak Bogoliubov excitations [39][40][41][42][43][44][45].
In this work, we propose an original combination of real and reciprocal space topologies, creating a truly protected pseudospin current in a bosonic system. Here, the topological phase and the edge spin currents are not protected by a symmetry of the Hamiltonian, but by the real-space topology of the quantum vortices. We consider a BEC at the Γ point of the Brillouin zone of a QVH system based on two staggered honeycomb lattices. We demonstrate the existence of a coupling between the winding number of a vortex and the valley of the bulk Bloch band. This coupling leads to chiral vortex propagation at an interface between two regions with inverted staggering, where the winding-valley coupling provides true topological protection against backscattering, contrary to the interface states of the non-interacting Hamiltonian. This configuration can be seen as a Z 2 topological insulator, similar to the quantum spin Hall effect [9], but where the role of spin is played by the winding of the vortices. Our results apply to polariton condensates in recently fabricated polariton honeycomb lattices [46] and to atomic BECs in optical lattices [47].
Non-interacting QVH. We consider an interface between two honeycomb lattices with opposite staggering, each of them being well described by a tight-binding (TB) Hamiltonian: where 2∆ = E B − E A is the energy difference between the ground states of A and B sites and J is the nearest neighbour tunnelling coefficient. A non-zero ∆ leads to the opening of a bandgap and implies the presence of opposite Berry curvatures in K and K valleys. If the gap is sufficiently small, the Berry curvature is strongly localized in each valley which allows to compute the valley Chern numbers: C K,K = ±0.5. The number of chiral states in each valley at the zigzag interface between the opposite lattices is defined by the domain wall topological invariant [48]: N K,K = C K,K (l) − C K,K (r) = ±1 (where l and r stand for the left and right domains). This results in the presence of one chiral state in each valley with opposite group velocities (QVH effect). However, these valley states, degenerate in energy, are not protected by some specific symmetry, which means that the backscattering due to diffusion from one valley to the other is not forbidden for single particles. Quantum vortices. The BEC can be described by a single-particle wavefunction (WF) ψ (the order parameter). In the mean-field approximation, ψ is the solution of the Gross-Pitaevskii equation (GPE), including interparticle interactions: where m is the particle mass, α is the interaction constant, U is the external potential, and µ is the chemical potential of the condensate. The existence of ψ imposes the irrotationality of this bosonic quantum fluid: ∇ × v = 0 everywhere, except zero-density points. The condensate velocity is given by v = ∇ϕ/m (ϕ = arg ψ). The phase winding around the zero-density points where ψ = 0 is fixed by the single-valuedness of ψ: ∇ϕdl = 2πp, where p is the winding number. The solutions with non-zero p are called vortices, and their characteristic size is determined by the healing length ξ = / √ 2αnm. We will consider only single-winding vortices (p = ±1), because they are energetically stable. We are going to study such vortex solutions in a staggered honeycomb lattice.
Winding-valley coupling. First, we shall demonstrate that the core of a vortex with a given winding corresponds to a certain valley (K or K') of the single-particle dispersion of staggered graphene, that is, the existence of winding-valley coupling for vortices.
Let us consider the core of a sufficiently large vortex (ξ a, where a is the distance between nearest neighbors), where the density is necessarily small and the interactions can be neglected. To minimize the on-site energy given by E = E A |ψ A | 2 +E B |ψ B | 2 , the WF is mostly localized on the atoms of the A type, which have lower energy (assuming E A < E B ). In the limit of a large gap, ∆ J, only the A-atoms are populated. The WFs in a periodic lattice can be written as a product of a Bloch function and a plane wave. For the hexagonal lattice, the Bloch part of the WF determines the densities and phases on the A and B atoms. Therefore the Bloch function in the vicinity of the vortex center is (1, 0) T . We can obtain the corresponding plane wave by Fourier transform of the WF ψ (k) analytically in the TB approximation (see [49] for details). We find that the maximum value of the WF is achieved for k = K and k = K , depending on the vortex winding p. Thus, both the Bloch wave and the plane wave part of the WF in the core of a vortex of a given winding define a state corresponding to a certain well-defined valley of the single-particle dispersion, and there is a winding-valley coupling for sufficiently large vortices which reads: where τ = ±1 is the valley number and s is the lattice staggering (s = +1 for E A < E B and s = −1 for . This result is linked with the well-known optical selection rules in Transitional Metal Dichalcogenides [50] where the phase pattern at the K point exhibits an angular momentum for each unit cell, which determines the angular momentum of photons absorbed for a given valley. To confirm our analytical TB solution, we have performed numerical simulations by solving the GPE beyond the TB approximation, with an explicit honeycomb lattice potential U (r). To find the WF of the vortex, we have introduced the relaxation term [51], preserving zeroes of the WF. The results of these calculations are  Figure 1: Numerical density profile of the vortex stationary solution in real(a,b,c) and reciprocal (d,e,f) space for different filtering scales (w = 7, 3, 1 µm, respectively). All parameters as in Ref. [26].
shown in Fig. 1. We filter the vortex core using a Gaussian function of size w (panels (a)-(c)). For large w, the image in the reciprocal space ( Fig. 1(d)) is dominated by the condensate centered at the ground state (Γ point).
The ground state itself is empty, because the vortex imposes v = 0 everywhere. For smaller w ( Fig. 1(e,f)), the core of the vortex is centered at the K points of the reciprocal space, while the K valleys are empty. Opposite results are obtained for opposite winding, confirming the valley-winding coupling for vortices. Vortex at the interface. We have shown that the vortex WF in the reciprocal space is composed of 2 important contributions. Most of the particles of the condensate, far from the vortex core, are concentrated around the Γ point (small k). These particles are practically unaffected neither by the presence of the lattice, nor by any possible interfaces. On the other hand, the core of the vortex is at the K point, and we can expect interesting effects linked with the interfaces, where in the linear regime the states from the bulk K points give rise to chiral propagative interface states (QVH states). We shall therefore calculate analytically the energy of the vortex as a function of its position and wavevector of the core, using the TB approximation.
In a general case, the energy of the vortex can be calculated using the grand canonical expression [2]: Qualitatively, this expression is the difference between the energy of a system with a vortex and the energy of a system without a vortex (but with a condensate in the ground state with the unperturbed density n). The first step is to split the integral into 2 regions: the core (|R| ≤ ξ) and the outside zone (|R| > ξ). In the second region, |ψ| 2 ≈ n, and the only contribution to the vortex energy comes from the kinetic energy term, which gives the well-known logarithmic expression [53] E r>ξ v = πn 2 ln (1.46R 0 /ξ) /m (R 0 is the system size). In the vortex core, the presence of the lattice has to be taken into account. As we have shown above both analytically and numerically, the core of the vortex is a wavepacket centered at a wavevector k 0 close to either K or K (we take a Gaussian wavepacket ψ G ). We calculate its energy versus k 0 using the TB dispersion E(k). However, the X spatial direction, perpendicular to the interface, has to be treated in the real space (x 0 is the vortex center). The contribution to the kinetic energy is calculated as: E kin,r<ξ v (x 0 , k 0 ) = x0+ξ x0−ξ dx dkψ * G ψ * 0Ĥ ψ 0 ψ G , where ψ 0 (x, k) are the single-particle eigenstates of the lattice. These eigenstates are quantized in the X direction. Their spatial overlap with the vortex core plays an important role. For the delocalized bulk states the overlap tends to zero with increase of the stripe width. On the other hand, the state localized at the interface (width κ) has a non-vanishing overlap and the contribution of this state dominates the dispersion of the vortex core. An example of such dispersion in the vicinity of the K and K points is shown in Fig. 2(a,b): the dispersion of the core (blue line) inherits the dispersion of the linear eigenstates at the interface (red dots), and therefore their valley-dependent propagation direction (chirality), as compared with the non-propagating bulk states with zero group velocity exactly at K or K (black points).
The kinetic energy of the core also depends on the position of its center x 0 : if the core is perfectly superposed with the interface state (centered at the interface), the energy at k 0 = K is exactly the same as that of the interface state. On the other hand, if the core is located in the bulk, its energy is that of the top of the valence band, determined by the energy splitting E kin (x 0 , k 0 ) = −∆. The interface therefore represents a barrier, if only the kinetic energy is taken into account.
The contribution of the interactions to the vortex core comes from the sensitivity of the vortex to the local changes of the density in the condensate. In the vortex core, the density |ψ| 2 is small as compared with the background density n(r), and the integral reads: E int,r<ξ v = ξ 0 αn 2 πrdr. Thus, the vortices are attracted to lower-density regions minimizing the total energy of the system. The density of the condensate without a vortex depends on the local potential, which affects the density of the condensate at the scale given by the healing length ξ. Considering the interface as a Delta barrier V 0 δ(x), the density of the condensate in its presence can be found as [52]: n(x) = n 0 (1 − cosh 2 ((x c + |x|)/ξ )), where x c and ξ depend on V 0 . The interaction energy of the vortex core as a function of x 0 therefore exhibits a minimum of the width ξ ≈ ξ.
The sum of the kinetic and interaction energy depends on the parameters of the system. An example of such dependence as a function of x 0 is shown in Fig. 2(c) for ξ > κ. In this case, the vortex can be localized on either side of the interface, the latter acting as a barrier preventing the vortex to go to the other side of the interface and change valley. We see that the properties of the single-particle dispersion of the interface states are inherited by the vortex solution of the non-linear equation via the core.
Our analytical results are again fully confirmed by numerical simulations of vortex propagation along the interface using Eq. (2). The snapshots of one of such simulations are shown in Fig. 3 (see [49] for movies). We see that the vortex remains attached to the interface and propagates along, without being scattered backwards on the corners. An additional defect of 1 meV has been added on an interface pillar for comparison with the lin- ear case, where it led to strong backscattering [26], which allows us to check that the vortex is indeed immune to backscattering thanks to the additional topological protection provided by its winding via the winding-valley coupling. For direct comparison, all parameters used in numerical simulations were exactly the same as in [26] (except interactions [49]). It allows to obtain the group velocity of the interface states v g = ∂E/∂k = 0.7 × 10 6 m/s or 0.7 µm/ps. This is the velocity with which the WPs at the interface can be expected to propagate in this particular lattice. Interestingly, the vortex velocity is different from v g . We stress that it is also different from what can be calculated for the vortex rolling effect (see [49]). Indeed, in our calculation we were assuming that only one type of the atoms is occupied for a given staggering. However, as shown in a scheme in Fig. 4(a), the interface represents a violation of a perfect staggering, and thus the higher-energy sublattice acquires a density estimated as n = 2n/ 1 + ∆ + √ ∆ 2 + 4J 2 2 /4J 2 (see [49]). The resulting velocity, reduced with respect to that of the linear interface states, is given by: We plot the dependence of v on the pillar size ratio ∆R/R (determining the gap size ∆) in Fig. 4(b). Red dots show the results of numerical simulations. Black line is the analytical solution given by Eq. (5), where v g and ∆ are taken from numerical simulations in linear regime. We see that it corresponds almost perfectly to the points (exact numerical solution) while there are no fitting parameters. This confirms the validity of our interpretation.
Conclusions. We demonstrate a vortex-valley coupling for a BEC in a staggered honeycomb potential. The main consequence of this property is the robust chiral propagation of vortices at the zigzag interface between two lattices with opposite staggering. The vortices, contrary to the linear WPs, are immune from backscattering thanks to their real-space topological protection. Hence, this work highlights a new combination of real and momentum space topology. These results are promising for the development of a new field of vortextronics, where the information will be carried by vortices. The possibility to create chiral pathways for vortices and to automatically sort them according to their winding is crucial for information treatment.
We acknowledge the support of the project "Quantum Fluids of Light" (ANR-16-CE30-0021), of the ANR Labex Ganex (ANR-11-LABX-0014), and of the ANR Labex IMobS3 (ANR-10-LABX-16-01). D.D.S. acknowledges the support of IUF (Institut Universitaire de France). tice, because it corresponds to large distances and small wavevectors (long-wavelength approximation). It does not depend on the vortex position, because the local potential does not affect the overall rotation of the particles dominating this energy. Neither does it depend on the propagation of the vortex, because the condensate far from the core remains globally unperturbed by this propagation by definition, otherwise the situation would correspond to the propagation of the vortex with a flow.

SUPPLEMENTAL MATERIAL
In this supplemental material, we present additional details on the derivation of results of the main text. We discuss the winding-valley coupling and the velocity of a vortex at an interface. Finally, we comment the supplemental video files.

Vortex-valley coupling
The calculation of the Fourier transform ψ(k) from the main text is carried out as follows. In the TB approxi-mation, ψ(r) is defined only in discrete points in space, and the integration is replaced by summation. Studying the core only, we take into account only the 3 atoms of the A type of the central hexagon.
This gives the following sum: where p = ±1 is the vortex winding. This expression can be rewritten as To simplify the expressions, let us define the arguments of the two exponents as separate variables: We can then find the position of the maximal probability density in the reciprocal space |ψ(k)| 2 , which writes (by separating the real and imaginary parts): ψ p (k) 2 = 1 + cos 2 η p + cos 2 ζ p + 2 cos η p + 2 cos ζ p (10) + 2 cos η p cos ζ p + sin 2 η p + sin 2 ζ p + 2 sin η p sin ζ p which can be simplified to ψ p (k) 2 = 3 + 2 (cos η p + cos ζ p + 2 cos η p cos ζ p ) (11) The maximal value of this expression is achieved when both η p = 2πν and ζ p = 2πµ, where ν and µ are integer numbers. From the latter, taking for example ν = 0, it is easy to obtain, for p = 1, k y = K (where K = 4π/3 √ 3a), and k x = 0, and for p = −1, k y = −K and k x = 0.

Vortex velocity
We have studied how the vortex velocity depends on the parameters of the system in order to check that the propagation along the interface is not linked with the well-known vortex rolling effect. First, let us see that the vortex really follows the interface, and its core is located exactly within the unit cell, which separates the two inverted materials. Figure S5 shows a snapshot of the phase of the wavefunction with a vortex. A 2π phase jump line is clearly visible, and the core of the vortex is located at the end of this line. The rotation direction of the vortex is shown with a red arrow, and the green arrow indicates the propagation direction of the vortex along the interface (white dashed line). We see that the edge of the phase jump line is within the unit cell located at the interface.
One might think that the vortex is simply rolling along the interface, like a wheel, converting rotation into propagation. The characteristic distance at which the density can vary in the condensate is given by the healing length ξ and therefore the center of the vortex in this "rolling wheel" image has to be located at a distance ξ from the wall, which allows to find the speed of rotation of the particles where they meet with the wall (and therefore the vortex propagation speed) using the expression v = m where one takes r = ξ, which gives simply that the vortex propagates with a velocity roughly equal to the speed of sound in the condensate v = αn/m. In this model, one could therefore expect a pronounced dependence of the vortex propagation velocity on the particle density.
Another alternative could be that the vortex simply propagates with the group velocity of linear states at the interface, which can be calculated from the dispersion, as discussed in the main text. Figure S6 compares the predictions of these models as a function of interaction energy αn with numerical results (black squares). Clearly, the simple predictions of the two naive models (red circles for rolling effect and black dashed line for the linear group velocity) strongly deviate from numerics. The  model of the rolling wheel (red dots) predicts a dependence on the density which is not observed at all (the interaction energy changes by a factor 5, and there is no significant change of the vortex velocity). The group velocity of the interface states strongly overestimates the real vortex propagation speed (also by a factor 5).
To calculate the vortex velocity, we analyze the currents that take place within its core (concentrated in a given valley because of winding-valley coupling). In the bulk, the valley states are not propagating, but rotating, because the 3 quantum-mechanical current terms between the 3 pillars of the same type which have different phases (0, 2π/3, 4π/3) exactly compensate each other, as these are three identical vectors rotated at 120 degrees. Indeed, j = n m ∇ϕ where n is the particle density, and therefore, to calculate current in the tight-binding approach we need to consider only pillars with nonzero density and take into account the phase difference between each pair. At the interface the situation changes, as can be seen in Fig. 4(a) of the main text. The A pillars on the left of the interface are not large pillars (with lower energy) but small pillars (with higher energy), and therefore, the 3 current terms (blue arrows) do not have the same prefactor. The phase differences are the same, but the density on the pillars on the left of the interface is smaller (it is not zero as it would be in the bulk, because the presence of the interface mixes the Bloch states), and therefore the current term marked as a dashed line has a smaller magnitude than the other two. This results in a net current pointing upwards, and this is what leads to the propaga-tive nature of the interface states.
The total current reads Assuming that the density on the A pillars on the right of the interface is n and the density on the A pillars on the left of the interface is n , we can write the magnitude of the current terms as: and The orientation of the vectors makes that the X projection of j 3 is 0, while the X projections of j 1 and j 2 are opposite, and so they compensate each other. The Y projections give: which finally gives Without the interface, n = n and j = 0, as expected. The presence of the interface makes n < n. If we consider an isolated problem of two pillars with coupling J and energy splitting ∆ (which determines the gap in the bulk TMD analog), we can estimate n as n = 2n 1 + ∆ + √ ∆ 2 + 4J 2 2 /4J 2 (13) which finally gives the expression for the group velocity of the main text, because 2π /m/3 √ 3a is simply an estimate of the group velocity v g in terms of the tight-binding parameters.
We can also calculate an approximated expression, assuming that ∆ J, n = n 1 − ∆ 2J (14) which gives for the net velocity along the interface The corresponding calculated velocity shown in Fig. S2 by a solid black line corresponds well to the numerical results, contrary to the predictions of the simple models.
In the opposite limit of very large ∆, This expression also increases with the increase of ∆. It is interesting to see that this expression is bounded from above by a limiting value, which cannot be exceeded by changing ∆ (but only by changing J, which affects m).

SUPPLEMENTAL VIDEO
In the supplemental video file vortexdefect.avi (also available at https://www.youtube.com/watch?v= PNsDF5xUvH4), we show the temporal evolution of the spatial density distribution of the condensate |ψ(r, t)| 2 , obtained by direct solution of the Gross-Pitaevskii equation with U (r) being the lattice potential, without the tight-binding approximation. The snapshots from this movie are shown in Fig. 3 of the main text. The vortex is attached to one side of the interface and propagates along it, passing around two corners and a defect.
A second movie linwp.avi (also available at https://www.youtube.com/watch?v=M7nbL5i9l44) demonstrates that a linear Gauss-Laguerre wavepacket with a non-zero angular momentum does not at all exhibit the same behavior as the vortex in an interacting condensate: the wavepacket is unstable and expands rapidly, preventing the observer to keep trace of the propagation of its center. The features of the interacting BEC maintaining the vortex are therefore crucial for the results obtained in the main text.