{'fraction_non_alphanumeric': 0.035753027384826545, 'fraction_numerical': 0.010803391556529632, 'mean_word_length': 4.478714372346877, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Vortex states in a non-Abelian magnetic field

A type-II superconductor survives in an external magnetic field by admitting an Abrikosov lattice of quantized vortices. This is an imprint of the Aharonov-Bohm effect created by the Abelian U(1) gauge field. The simplest non-Abelian analogue of such a gauge field, which belongs to the SU(2) symmetry group, can be found in topological insulators. Here we discover a superconducting ground state with a lattice of SU(2) vortices in a simple two-dimensional model that presents an SU(2)"magnetic"field (invariant under time-reversal) to attractively interacting fermions. The model directly captures the correlated topological insulator quantum well, and approximates one channel for instabilities on the Kondo topological insulator surface. Due to its simplicity, the model might become amenable to cold atom simulations in the foreseeable future. The vitality of low-energy vortex states born out of SU(2)"magnetic"fields is promising for the creation of incompressible vortex liquids with non-Abelian fractional excitations.


I. INTRODUCTION
The ubiquitous Rashba spin-orbit coupling, found in two-dimensional electron systems and topological insulators (TIs), is equivalent to a non-Abelian SU(2) "magnetic" field without dynamics 1 . This effective "magnetic" field can be quite strong 2 and responsible for topologically protected edge states in two-dimensional TIs 3-6just as the ordinary U(1) magnetic field is responsible for edge states in quantum Hall systems. Here we examine another possible analogy between U(1) and SU(2) fields: can a two-dimensional superconductor subjected to the Rashba spin-orbit coupling exhibit a vortex lattice like type-II superconductors in magnetic fields?
Spin-orbit-coupled superfluids have been studied in the context of cold atoms, motivated by the experimental progress in subjecting trapped two-component bosonic atoms to artificial gauge fields [7][8][9][10] . A genuine Rashba spin-orbit coupling has not been created yet, but new technologies based on atom chips could realize it in the foreseeable future. Theoretical studies have mostly scrutinized the Rashba spin-orbit-coupled condensates of bosons that can move freely through the trap. Mean-field calculations favor "plane wave" and "standing wave" condensates formed at finite wavevectors [11][12][13][14] , but quantum fluctuations may give rise to more unconventional states, perhaps even with fractional quasiparticles 15 .
In this paper, we consider a two-dimensional lattice model of attractively interacting fermions that can form a paired spin-triplet superfluid in the presence of a Rashba spin-orbit coupling. The model describes a twodimensional TI placed in contact with a conventional superconductor 2,16 , but it can be reformulated to describe triplet exciton pairing among repulsively interacting electrons. The surface of a Kondo TI SmB 6 can also exhibit instabilities 17 that are addressed by this model. The choice of a lattice rather than a continuum model makes vortex arrays more competitive to the "plane wave" and "standing wave" states at the mean-field level. Indeed, we find condensates with a vortex lattice in a broad parameter regime when vortex cores of small size reside inside lattice plaquettes. The presented results follow from a numerical mean-field calculation at zero temperature. Quantum fluctuations enhanced by the low spatial dimensionality can give rise to vortex liquids.
The significance or our results is two-fold. First, the discovered vortex lattice is a novel state of matter -an Abrikosov-type array of SU(2) topological defects that features spin-supercurrent loops and respects the timereversal (TR) symmetry. The present study demonstrates a realistic microscopic model where such a vortex state can be stable. This is a different yet remarkably similar vortex array to the one identified by the Landau-Ginzburg analysis 18 as a prominent metastable state in the continuum limit. Related but different states have been seen in a few studies of two-component interacting particles with spin-orbit coupling [19][20][21] .
Second, the quantum melting of a vortex lattice can produce an incompressible quantum liquid 22 whose excitations have fractional exchange statistics. This was numerically demonstrated in the context of rotating bosonic atoms, where the quantum melting of a U(1) Abrikosov lattice yields a fractional quantum Hall state [23][24][25][26][27] . Similar physics is expected in the present SU(2) system 2, 28,29 . Both cold atom gases and quantum well heterostructures enable experimental manipulation of the quantum state. Destabilizing an existing vortex lattice at low temperatures by almost any means would allow experimental access to a fractional quantum liquid 2 . The non-Abelian fractional states created by the SU(2) Rashba flux could be useful for quantum computing. This work identifies the precise nature of the SU(2) vortices, which informs about the fractional statistics of quasiparticles in the corresponding vortex liquids 29 .
The paper is organized as follows. Section II A introduces the model, discusses its features and relates it to realistic physical systems. Section II B presents the main results -the mean-field phase diagram of SU(2) vortex states and other competing orders. Section II C qualitatively analyzes fluctuation effects, and explains why incompressible quantum liquids with fractional excitations should also exist in the phase diagram of this model. The main discussion ends with conclusions in Section III. Certain special subjects and various technical details are discussed in Appendices.
The model we study is given by the imaginary-time path integral expressed in terms of the four-component fermion (Grassmann) spinors ψ and three-component bosonic (complex) spinors η. These fields live on a square lattice whose sites are labeled by i. The fermions ψ iστ have spin σ ∈ {↑, ↓}, and another two-state degree of freedom τ = ±. The latter is motivated by the Hamiltonian representation of an ultra-thin TI film 2 . A thick TI slab hosts protected helical Dirac-metal states on its top (+) and bottom (−) surfaces 30 , with opposite helicity of spin-momentum locking. When these surfaces are brought close together in a film, the Dirac spectrum acquires a gap through intersurface tunneling (∆), but the low-energy quasiparticles retain the τ degree of freedom associated with their surface belonging. The non-interacting Hamiltonian (2) describes fermion hopping on the tight-binding square lattice, in the presence of a static SU(2) gauge field defined on lattice links: (3) This gauge field, written in terms of the Pauli matrices σ x,y,z that couple to spin, implements the Rashba spinorbit coupling of strength a. The Pauli matrices τ x,y,z operate on the "surface index" τ . The fermions' SU(2) charge τ z depends on the surface they visit, and the ∆τ x term directly describes inter-surface tunneling that creates a bandgap. The fermion density is controlled by the chemical potential µ in our formalism.
Interestingly, the 100 crystal surface of a Kondo TI SmB 6 has similar features. It hosts three helical Fermi pockets in the vicinity of Γ, X and Y points of the surface's 1 st Brillouin zone [31][32][33][34][35] . The Γ pocket has the opposite helicity of spin-momentum locking than portions of the X and Y pockets 17 . Collective excitations seen by neutron scattering in SmB 6 are precisely poised to couple these pockets in either Cooper or exciton channels 36,37 . Any coupling between the pockets of opposite helicity can lead to triplet pairing of the kind we investigate in this paper.
The SU(2) flux Φ of (3) on a square plaquette, with corners 1, 2, 3, 4 enumerated in a clockwise sense, can be defined as a matrix by: Depending on which of the four corners is labeled by 1, one finds Φ = arccos 1 − 2 sin 4 (a) with four different combinations of ± signs. The leading order term Φ → 2a 2 τ z σ z in the limit of small a is finite and independent of the exact path around a plaquette; this becomes the "magnetic" Yang-Mills flux in the equivalent continuum-limit single-particle Hamiltonian: Here, t = 1/2m, µ ′ = µ + 4t, v = a/m, A → −mv(ẑ × σ) and σ = (σ x , σ y , σ z ), using the units in which and the lattice constant are equal to 1. The eigenvalues of the matrix Φ are SU(2) gauge-invariant and their finite values lead to an Aharonov-Bohm effect just like the analogous U(1) magnetic field, even though tr(Φ) = 0 due to the TR symmetry.
Interactions among the fermions are generated by a Hubbard-Stratonovich field, the S = 1 bosonic spinor η: This is the part of the most general decoupling of short range attraction that leads to Cooper pairing in spintriplet channels 2 . The remaining part, responsible for singlet pairing, is discussed in the Appendix and otherwise neglected -singlets are completely suppressed at the mean-field level in the portion of the phase diagram that we analyze, and, more importantly, they are impaired by spin-orbit SU(2) neutrality. Triplet Cooper pairs are able to gain energy from the Rashba spin-orbit coupling in proportion to the momentum they carry 38 . The model (1) is invariant under TR and has the full point-group symmetry of the square lattice. The only continuous symmetry of the model for generic a is the U(1) symmetry responsible for charge conservation. This symmetry is spontaneously broken in a superfluid state. However, the continuous symmetry group is larger if a = n π 2 , n ∈ Z. At a = 0, the spin "up" and "down" states are decoupled, so one can apply an independent global U(1) transformation on each spin projection. At a = π/2, the nearest-neighbor hopping terms in (2) reduce to ∝ ψ † i σ x,y ψ j , which allows the following sublattice-dependent symmetry transformations on a bipartite lattice: This is a double global U(1) transformation specified by two angles θ A and θ B , where A, B indicate the lattice sites belonging to one or the other sublattice. All properties of the model are periodic under a → a + 2π, while the change a → a + π is equivalent to changing the sign of t. The spectrum is invariant under a → −a, since this corresponds to the lattice inversion through a point (without affecting spin).

B. The mean-field phase diagram
In this paper we analyze the zero-temperature phase diagram of the model (1), beginning with the mean-field approximation. We consider various spatially dependent but static order parameters η = (η ↑ , η 0 , η ↓ ), and search for the one that minimizes the ground state energy at any given values of µ, ∆, a, U t , U t0 , etc 39 . This is equivalent to the usual self-consistent solution of the Bogoliubov-de Gennes equations for a superconductor. We assume that the order parameter η belongs to one of these categories: 1. spatially periodic with a unit-cell of l x × l y sites: freely vary η i at every site of the unit-cell; 2. commensurate pair-density wave (PDW) at a wavevector Q = (q x , q y ): vary Q and the Fourier amplitudes of η; 3. incommensurate plane wave state (PWS) at a wavevector Q = (q x , q y ): vary Q and the amplitude of η.
The first category is the most general one among all periodic states, but computationally demanding if the unitcell is large. For example, we need 12 complex numbers to specify a periodic state with a l x = l y = 2 unit-cell. The second category is a subset of the first category, but faster to compute. The PDW order parameter is restricted to the form: where r i = (x i , y i ) are the integer-valued coordinates of the lattice site i, and Q is a commensurate wavevector in the first Brillouin zone (i.e. there is an integer n such that nQ is a reciprocal lattice vector). It takes 6 complex numbers to specify the two triplet amplitudes η ± , and another two rational numbers to specify Q. Finally, the third category captures the plane wave subset of all category-2 states: It takes 3 complex numbers to specify the amplitude η pw , and two real numbers to specify Q. A plane-wave state can be handled even when it is incommensurate, since it respects the lattice translation symmetries up to a global gauge transformation. Within each of these categories, we construct the Bogoliubov-de Gennes Hamiltonian and diagonalize it to determine the superfluid meanfield ground state. Then, we vary all the quantities that specify the order parameter in its category until we find a local minimum of the ground state energy at any point in the phase diagram. At the end, the absolute ground state energies in all categories are compared to identify the best order parameter. A more technical discussion of this procedure and its accuracy can be found in the Appendix.
The main results are illustrated in Fig.1. This is the phase diagram of non-trivial superfluid states that take advantage of the spin-orbit coupling a, plotted as a function of the two interaction couplings U t , U t0 . The chemical potential µ = 0 is kept in the middle of the bandgap created by ∆ > 0, but interactions are strong enough to cause pairing. The first striking feature is the existence of several competing orders that reflect the frustrated dynamics. Frustration is created by the competition between the lattice constant (equal to 1 in our units), the SU(2) "magnetic" length (∼ a −1 ), and the pairing coherence length. The discovered orders have very similar energies, and differ by no more than 1 − 2% of the absolute energy value in many scrutinized parts of the phase diagram. The achieved absolute error of any calculated energy is much smaller than the energy differences between the found states.
The line U t ≈ U t0 separates two different regimes. For U t > ∼ U t0 , the symmetric triplet component η 0 of the order parameter costs less energy than the spinful components η ↑ and η ↓ . The ensuing η 0 condensates are found to be uniform and with completely suppressed η ↑ = η ↓ = 0. Conversely, the region U t < ∼ U t0 hosts condensates dominated by η ↑ and η ↓ , typically with η 0 = 0. The phase transition between the two types of condensates is first order within numerical resolution. The vanishing of η 0 = 0 when U t < ∼ U t0 explains the tendency The order parameter is dominated by the spinful components η ↑ , η ↓ for Ut < ∼ Ut0 and the "symmetric triplet" component η0 for Ut > ∼ Ut0. The interaction couplings Ut, Ut0 are expressed in the units of t −1 and µ = 0, ∆ = t, where t is the hopping constant of the Hamiltonian (1). The order parameter magnitude in the same units of energy is |η| ∼ t, implying that the results correspond to a strong coupling regime.
of the other phase boundaries to be vertical (the quasiparticles end up decoupled from U t0 in the mean-field approximation, so the phase boundaries do not depend on U t0 ).
The region U t < ∼ U t0 with η ↑↓ = 0 hosts TR-invariant vortex lattice condensates at a ≥ 0.8. The appearance of dense vortex lattices is related to the large SU(2) flux (5), maximized at a ≈ 0.96. A TR-invariant vortex lattice needs both η ↑ and η ↓ order parameter components to establish a pattern of spin currents that gain energy through the Rashba spin-orbit coupling. It is, then, natural that vortex states live only in the U t < ∼ U t0 regions. The vortex lattice states were found by the category-1 search using the minimal unit-cell of 2 × 2 sites. Their order parameter, shown in Fig.2, is TR-invariant and features a periodic structure of spin currents on lattice links: where S k (k ∈ {x, y, z}) are spin projection matrices in the S = 1 representation. Partial understanding of this state follows from the Landau-Ginzburg theory of twodimensional S = 1 superfluids with Rashba spin-orbit coupling 18 , which predicts that similar checkerboard vortices of J z spin-currents are metastable. In the continuum Landau-Ginzburg picture, the spin-orbit energy is actually gained through "helical" spin currents J x,ythe pattern of J x,y becomes singular and then drives the appearance of J z spin-current vortices. However, in the present calculation, η 0 = 0 renders J x,y = 0 in the vortex state, so the spin-orbit energy is gained by an intrinsically lattice mechanism. By gauge symmetry 2,38 , the triplets are effectively coupled to the SU(2) gauge field: where A ij = A k ij S k is of the same kind as (3) but constructed with the S = 1 spin operators S k instead of the S = 1 2 SU(2) generators 1 2 σ k , and A ij = A k ij A k ij (the summation over repeated upper indices is implied). The compact form of this coupling on the lattice has an additional term in the S = 1 representation involving bond scalars: where k, l ∈ {x, y}. Note that K kl ij = K kl ji = K lk ij are symmetric in both lower and upper indices. Only K xx and K yy couple to the gauge field (3), but this is enough to generate a lattice pattern shown in Fig.2, which breaks lattice translation and rotation symmetries. This is ultimately responsible for the spin-orbit energy gain and the stability of the J z vortex lattice. Quantum fluctuations inevitably lower the average density of singularities by vortex-antivortex annihilations and make room for vortex motion on the lattice, but cannot remove all singularities because that would produce a higher-energy pair-density-wave condensate.
The main competitors to the TR-invariant vortex lattices are various pair-density waves (PDWs). All found PDW states have ordering wavevectors along a lattice direction, e.g. Q = (Q, 0). The magnitude Q tends to grow with the spin-orbit coupling a as expected. However, the lattice commensuration effects introduce large  (14), with the same magnitude on all bonds but different sign on horizontal and vertical bonds; red and blue denote K > 0 and K < 0 on a given bond respectively. Note that K yy is the same as K xx rotated by 90 degrees.  continuum limit is approached when a is small and the ordering wavevector of the lowest-energy PDW is sufficiently small to be in the smooth (non-oscillatory) region of E a (Q). Among all these PDW states, the only one that ever wins as a mean-field ground state is the TRinvariant order at Q = π/2, seen for 0.7 ≤ a ≤ 0.875. The next strongest PDW is also TR-invariant at large a and Q = π, but it never wins. All other prominent PDWs found by category-2 calculations break the TR symmetry. Among the category-3 results, only the states at Q → 0 (small a) and Q ≈ π (large a) can minimize energy, but not enough to win as ground states; still, only such PWS can be TR-invariant, so the states with TR symmetry consistently have lower energy than those that break it. A separate category-1 calculation discovered additional metastable TR-invariant vortex lattices with a 4×4 unit-cell, in the range 0.725 ≤ a ≤ 0.775 where the Q = (π/2, 0) PDW states are stable. The most prominent vortex pattern is the same checkerboard as before, but with a twice larger spatial separation between vortices shown in Fig.3. It is interesting to note that the category-1 numerics fairly consistently found this particular type of local energy minimum instead of the global PDW minimum. Note that the explored phase-space was very large, since 16 × 3 × 2 = 96 real quantities were needed to represent the order parameter on a 4 × 4 unitcell. In contrast, only 12 real numbers are needed to specify the PDW order parameter of the true ground state discovered by the category-2 calculation. This probably indicates that vortex lattices have larger entropy (shallower and broader energy minimum) than the true PDW ground state. Quantum fluctuations beyond the meanfield approximation are expected to energetically enhance the states with larger entropy.
In the remaining scrutinized region 0.6 ≤ a ≤ 0.675, the lowest energy condensates are found to be uniform at all U t , U t0 . Even category-2,3 calculations settled for uniform states when a ≤ 0.65, despite having freedom to access smaller and incommensurate wavevectors. This is a consequence of being in the strong-coupling limit, with interaction energy scales U −1 t , U −1 t0 exceeding the bandgap ∆. It was necessary to work in this limit in order to numerically resolve the tiny energy differences between various competing states (see Appendix), but the price to pay was a strong suppression of either η 0 or η ↑↓ order parameter components in any given phase (interactions compete for electrons to form pairs so the dominant interaction decides the pairing channel, in agreement with the first order transition at U t ≈ U t0 ). Consequently, energy gain from the spin-orbit coupling requires finite bond scalars K kl (14) when η 0 = 0. Orders of K kl are hard to arrange at larger wavevectors when a becomes smaller, and their spin-orbit energy scales as a 2 . Seeing finite-momentum ordering at smaller a is probably possible only in the weaker coupling limit, which is not accessible with the current numerics. If interactions are not too strong, then the spin-orbit energy gain proportional to a in (13) can be sufficient to overcome the tendency of interactions to suppress η 0 . The resulting states can have finite helical spin-currents J x,y , and perhaps establish the vortex arrays envisioned in the Landau-Ginzburg analysis.

C. Fluctuation effects
The system studied here is two-dimensional, so one can anticipate significant corrections to the mean-field phase diagram. The continuous U(1) symmetry associated with superfluidity can be broken only at T = 0, while it gives way to algebraic correlations at low T > 0. All other symmetries of the model are discrete, and hence remain broken at low T > 0.
At special values a = n π 2 of the spin-orbit coupling, the model (1) has another continuous symmetry. The spinorbit coupling is absent if n is even, so the corresponding ground-states are trivially uniform. However, the condensates at odd n break the lattice translation symmetry, generally with a 2×2 unit-cell. The 2×2 vortex state that wins at a < π 2 becomes degenerate with the Q = (π, 0) PWS and infinitely many other states when a = π 2 . This follows from the symmetry (9), which allows free global rotations of the staggered θ pattern from Fig.2(b) in the opposite directions θ A = −θ B on the sublattices A and B. Since (9) is a continuous symmetry at a = π 2 , the degenerate condensates (e.g. vortex lattice and PWS) can be long-range ordered only at T = 0 and become algebraically correlated at low finite temperatures. Note that this special symmetry is unrealistic because it requires the absence of next-nearest and further-neighbor hopping in the model.
The mean-field energy advantage of the vortex lattice over its related PWS gradually disappears as a → π 2 . Since the related PWS state is obtained by annihilating all singularities of the vortex lattice, one may expect prominent local quantum fluctuations that annihilate and re-create vortex-antivortex pairs when a is close to π 2 . The annihilations reduce the average density of vortices, and hence create room for vortex motion on the lattice. One possible outcome is the condensation of vortices, which by duality leads to a Mott insulator 40,41 . However, if vortices fail to condense due to their large density, still comparable to the density of particles, then a melted vortex lattice is an incompressible quantum liquid with topological order -analogous to fractional quantum Hall states, but without TR symmetry breaking and quantized spin-Hall response 29 .
The crystal lattice helps to pin the vortices, so a quantum phase transition between the vortex array condensate and a band-insulator can be of 2nd order. The mean-field phase diagram features precisely such a transition when the intrinsic fermion bandgap ∆ becomes sufficiently large (comparable to the pairing energy scale). If vortices must rely only on mutual interactions to establish a long-range positional order, then a quantum Lindemann criterion 22 suggests that the mentioned 2nd order transition must be preempted by a 1st order quantum melting of the vortex array. Such melting occurs because the kinetic energy cost of localizing vortices exceeds the potential energy cost of deforming the vortex array when the order parameter amplitude becomes sufficiently small. For this reason, we expect that the actual phase diagram contains an intermediate vortex liquid phase between a vortex lattice condensate and a band-insulator.
Obviously, more analysis going beyond the mean-field approximation is needed to establish the accurate phase diagram of this two-dimensional system. However, the mean-field approximation is sufficient to reveal the existence of vortex ground states in the phase diagram of certain microscopic models like (1), which ultimately could be either crystalline or liquid. It is also adequate for establishing or confirming the predicted structure of vortices on the lattice, i.e. the patterns of charge and spin currents. The latter is a very useful information for anticipating the possible kinds of fractional quasiparticle statistics in putative incompressible vortex liquids. Field-theoretical arguments generally point to a particleflux attachment in incompressible quantum liquids 29 , so the non-trivial exchange statistics of quasiparticle excitations follows from the Aharonov-Bohm phase that one particle acquires on its path around the flux quanta attached to another quasiparticle. This simple principle may need to be combined with some deeper understanding of vortex dynamics in order to classify the possible topological orders in this system.

III. CONCLUSIONS
We studied the mean-field phase diagram of a simple two-dimensional model of attractively interacting fermions on the square lattice. The fermions experience a strong Rashba spin-orbit coupling in the manner characteristic for two-dimensional topological insulators. We discovered that a TR-invariant dense array of vortices has the lowest energy at the mean-field level among many competing orders in a large section of the phase diagram. This state is an analogue of the Abrikosov vortex lattice for the particular type of the SU(2) Yang-Mills flux that realizes the Rashba spin-orbit coupling. Other competing orders include pair-density waves and plane-wave states, some of which break the TR-symmetry. The energy competition between these states is fierce, so we expect that quantum and thermal fluctuations can easily destabilize some orders in favor of incompressible quantum liquids.
It should be pointed out that much of the observed phase diagram features are not universal. The phase transitions between the found competing orders are 1st order, and hence model-dependent. A somewhat different model may suppress or enhance the vortex states. Sorting out which microscopic details are favorable to vortex states is not a goal of the present work, and is probably not worth pursuing at present when correlated two-dimensional TIs are not experimentally available. However, our model is the simplest description of the physics that can be realistically expected in interacting TI quantum wells, on the boundaries of Kondo TIs, or engineered with cold atoms. In that sense, TI quantum wells are a promising system that may be able to realize the vortex arrays or liquids that we discussed in this paper.
Acknowledgements. This work was supported by the National Science Foundation under Grant No. PHY-1205571, and also in part by the National Science Foundation under Grant No. PHY-1066293 with hospitality of the Aspen Center for Physics.
The order parameters φ i+ and φ i− are singlet pairs of two electrons on the same Dirac cone, τ = +1 or τ = −1, while φ i0 is the inter-cone singlet. In general, no special symmetry governs the values of the interaction couplings U s+ , U s− , U s0 . Category-1 calculations were carried out to examine the competition between the singlet φ and triplet η order parameters in the present model. A general conclusion is as expected: the spin channel with strongest coupling (smallest order parameter suppression U ) wins. The transition between singlet and triplet orders is first order. The only interesting new feature of this model is that the spin-orbit coupling tends to strengthen and enhance the triplet orders. This is illustrated in the typical singlet-triplet phase diagram, Fig.5.
Note that even the η 0 -dominated phase in the region U t > ∼ U t0 is enhanced by the spin-orbit coupling, despite being uniform. This is a consequence of the quadratic spin-orbit effect in (13): both K xx and K yy are finite. One consequence is that the interesting vortex lattice and PDW phases are not helped by the spin-orbit coupling in their competition with the uniform η 0 phase. However, all triplet orders have a certain advantage over the singlet orders when the spin-orbit coupling is large.
Appendix B: The competition between pair density waves Fig.6 shows the phase diagram of competing PDW and PWS states. Arrows indicate (by length and orientation) the ordering wavevector Q of PWS and PDW states that were found to have the lowest energy by a category-2 or -3 calculation. Some PDW/PWS states were discovered by a category-1 calculation, in which case they are not accompanied by arrows; being defined in a 2 × 2 unitcell, their wavevector is then Q ∈ {(0, π), (π, 0), (π, π)}. The arrow color (not clearly visible everywhere even in the online color plots) indicates the relative amplitudes of the triplet spinor components. The amounts of R(red), G(green), B(blue) pigments in the RGB arrow color is set by the normalized amplitudes of η ↑ , η 0 , and η ↓ components respectively. Therefore, purple=red+blue arrows (hard to see) denote the PDW/PWS states dominated by η ↑ and η ↓ in roughly equal amounts when U t < ∼ U t0 . The sampled points in the phase diagram are represented by painted blocks with the following color scheme. Green colors are vortex arrays with four singularities in the unit-cell of 2 × 2 lattice sites. Blue colors are PDW states with non-zero order parameter amplitudes at two physically different commensurate wavevectors. Hence, a PDW state breaks translation symmetry. If a blue block is accompanied by an arrow, then the two wavevectors are Q and −Q (and not high-symmetry ones in the 1st Brillouin zone). Red colors are PWS states with a non-zero amplitude at a single wavevector Q (which can be a non-zero high-symmetry wavevector). A PWS respects translation symmetry up to a global U(1) gauge transformation. Light purple represents a conventional uniform superfluid condensed at Q = 0, and white is a band-insulator with a vanishing order parameter. Regions with spontaneous magnetization are painted with a lighter color, but follow the same color scheme. All states with spontaneous TR symmetry breaking are separated by a pair of black lines, where the thinner line faces the TR-invariant state. The size of blocks reveals the resolution at which the phase diagram was scanned.
The dense vortex array (Fig.2) discovered by the category-1 calculation can be regarded as a special kind of PDW states. Its lattice Fourier transform has finite amplitudes only at the wavevectors (π, 0) and (0, π). The metastable vortex array in a 4 × 4 unit-cell (Fig.3) has equal finite Fourier amplitudes at the wavevectors (± π 2 , 0) and (0, ± π 2 ). The qualitative distinction from the conventional PDW states is made primarily by the presence of circulating currents. The vorticity of such currents is counted by first extrapolating the order parameter to real-valued coordinates, and then measuring FIG. 6: A sample of the mean-field phase diagram of metastable pair density waves (PDW) and plain wave states (PWS), obtained by comparing the best results of category-2 and category-3 calculations. Stars indicate TR symmetry breaking, and thick black lines indicate a transition between TR-invariant and TR-broken states. Some regions of the phase diagram are noisy due to numerical minimization errors in a frustrated energy landscape, especially when the same PWS order at Q = (π, 0) is detected by both category 2 and 3 runs. However, almost all noise issues are resolved by the winning vortex lattice phase, not shown here (see Fig.1). the winding number of its resulting smoothly varying phase on the continuous paths around lattice plaquettes. This particular interpretation of vorticity in the lattice order parameter (whose phase may jump from one site to another) is somewhat arbitrary, but the pattern of circulating currents is physically measurable and imprints an Aharonov-Bohm phase on any degree of freedom that couples to it.

Appendix C: Numerical procedure and error estimates
The three category calculations carry out the same numerical procedure, and share the same code for common tasks. They mainly differ by the manner they define and handle the order parameter. Each type of the order parameter requires its own representation of the Bogoliubov de-Gennes (BdG) Hamiltonian for quasiparticle excitations.
The category-1 search specifies the order parameter on a real-space unit-cell of l x × l y sites. Both l x and l y are fixed before the search. The BdG Hamiltonian is represented in the Bloch-state basis corresponding to this unit-cell, while the unit-cell internals are handled in real-space. The category-2 search specifies its order parameter by Fourier triplet amplitudes at the commensurate wavevectors Q and −Q (which may be physically the same if Q is a high-symmetry wavevector). Its BdG Hamiltonian is also represented in momentum space, using the largest reduced 1st Brillouin zone allowed by Q. The category-3 BdG Hamiltonian is the simplest and structurally similar to that constructed for uniform condensates in the momentum space (pairing occurs between two fermions at momenta displaced by Q). The technical details of these BdG Hamiltonians are presented in the next section.
Once a BdG Hamiltonian is constructed, all three category searches diagonalize it and then obtain the energy of the many-body system (formally non-interacting after the mean-field approximation) by integrating out the single-particle BdG spectrum up to the Fermi energy µ in the appropriately reduced 1st Brillouin zone. Since the model has well-defined cut-offs, this numerical inte-FIG. 7: The difference ∆Ei = Ei − E between energy Ei calculated by the category i = 1, 2, 3 method and the reference accurate energy E of 36 randomly chosen states, plotted as a function of the number N of sampled points per spatial dimension for the 1st Brillouin zone integration (see text for the explanation). Each connected sequence of circles in a single color corresponds to a different order parameter (with varied ordering wavevectors and unit-cell sizes), but the energy of every chosen state was calculated with all category methods that can handle such a state. The category-1 energy at N = 100 is used as a reference E. The plot is meant to illustrate the pace and manner in which the calculated energies converge to the accurate values.
gration is the only source of errors in the calculations of energy. The integration employs Gauss-Legendre integration formulas, based on sampling the integrated functions at a set of discrete points in the 1st Brillouin zone; the sampled function is interpolated by a polynomial that can be integrated out exactly. One controls the absolute error of the calculated energy by the number N of sampled points per spatial dimension in the 1st Brillouin zone (1BZ), as shown in Fig.7.
The ground state search employs the simplex method to minimize the state energy as a function of state parameters. Each category represents the superfluid state by its own limited set of n variable parameters, as discussed earlier. The simplex method keeps track of the function values at n + 1 points (vertices of the simplex) in the abstract n-dimensional variable space. One vertex of the simplex is moved to a new position in each iteration of the method, according to a set of rules that gradually make the simplex drift toward and shrink around a local minimum in the parameter space. A search for the local minimum begins by specifying an initial "seed" state (the center of the simplex) and the overall size of the simplex (the distance between vertices). The search ends after a number of iterations when the simplex shrinks below a requested minimum size. The typical accuracy of the obtained order parameters is ∼ 0.1% of the absolute order parameter magnitude. In order to speed up the phase diagram mapping, all discovered ground states are recorded for future use as "seeds" on the nearby points of the phase diagram. A few independent phase diagram scans, and many independent calculations on arbitrarily selected points, were conducted with random "seeds" in FIG. 8: A typical energy scale E of the lowest energy condensates in all three categories on a path through the phase diagram. The units for E are given by the hopping t. The best states in the category 1,2,3 are plotted by the green, blue and red lines respectively, but these lines can hardly be distinguished because they coincide when plotted on the overall energy scales.
an attempt to reduce the danger of falling into metastable states.
The search for the minimum-energy states was carried out with a reasonable level of accuracy that did not cost an extreme computation time. This was N 1 = 8, N 2 = 8 and N 3 = 100 sampled 1BZ points per dimension for the categories 1, 2, 3 respectively. The category 1 is found to be very accurate even at N ∼ 10, while the categories 2 and 3 reach a similar accuracy only at N > 50 (see Fig.7). For the great majority of these order parameters, the three category calculations of energy agree to a fraction of a percent of the absolute energy value with N = 100, making absolute errors of about (1 − 2) × 10 −4 t. Energy minimization at high accuracy is prohibitively slow (except in category-3), so the optimal order parameters were obtained by using a manageable accuracy N = 8. However, the energies of discovered order parameters were recalculated at the end with high accuracy N = 100. All presented conclusions about the phase diagram follow from the comparisons of such accurate order-parameter energies in different categories. It is very unlikely that a highly accurate minimization would shift the order parameters enough to qualitatively change the phase diagram.
The typical energy scale of the paired ground states is shown in Fig.8. This plot compares the best condensates found in all three categories, but their energies cannot be distinguished on the absolute energy scale. This illustrates a significant frustration of the model and a fierce competition between different orders. The absolute energy difference between competing orders is often just ∼ 10 −3 t, but this is sufficiently larger than the absolute error to reliably determine the ground state.
The category-2 search is particularly complicated and tedious because one variable property of the order parameter is its commensurate ordering wavevector Q. The discrete components of Q must be varied outside of the simplex energy minimization. We systematically scan for the metastable states at different Q in the order that generally favors smaller unit-cells. The category-2 algorithm searches several times for an energy minimum at every checked Q. These searches are at first done with a few restricted forms of the category-2 order parameter that were observed to often minimize energy. Then the best found state is used as a "seed" for the final leastrestricted category-2 search. Each search starts with the same finite simplex size, to look for the minimum within a significant neighborhood of the "seed". Only the best overall result is ultimately reported for each Q as a metastable state energy. New wavevectors Q are checked down the line until the lowest discovered metastable energy has not changed in a while. This often reaches the unit-cells with up to 4 × 5 or 3 × 7 sites.
The category-1 search also tests several candidate states at each point in the phase diagram. It first obtains the lowest possible energy of all PDW/PWS states that can fit in the specified unit-cell of l x × l y sites (this step was skipped in the calculations on the 4 × 4 unitcell). Then it constructs restricted order parameters that loosely correspond to type-1 and type-2 states discussed in Ref. 18 , and minimizes their energy. Each new attempt uses the previous best result as a "seed". The final minimization is restricted only by the specified unit-cell dimensions, and only the best overall state is ultimately reported.
Given this aggressive procedure, one can be reasonably confident that the ordered states with small unitcells have been captured correctly. For example, we can be reasonably confident that a PDW or PWS with a 2×2 unit-cell is not the global energy minimum in the regions of the phase diagram where the TR-invariant 2 × 2 vortex lattice was reported as a consistent winner. We can be still fairly confident that no PDW/PWS state with a slightly larger unit-cell is the global minimum in these regions. Searching as long as it takes for the most optimal ordering wavevector Q of a PDW is futile, so we can never rule out the existence of global PDW minimums with larger unit-cells. However, considering also the local quantum or thermal fluctuations, the ordered states with larger unit-cells are almost invariably destabilized more than the states with smaller unit-cells. Fluctuations are certain to be prominent in our very frustrated system. The category-1 periodic order parameter is specified by complex triplet amplitudes η = (η ↑ , η 0 , η ↓ ) at every site of a fixed unit-cell. The fermion fields ψ are organized into four-component spinors, but representing their triplet pairing requires an eight-component Nambu spinor at each site i: Expressing all terms of (2) and (8) that contain fermion fields as a product of the Nambu spinors Ψ † and Ψ yields the BdG Hamiltonian: where the matrices h ij are defined only locally (i = j): and for the pairs of nearest-neighbor sites h i,i+x = −th x , h i,i+ŷ = −th y : with C a = cos(a), S a = sin(a). Including interactions in the uniform singlet pairing channels of equation A1 amounts to adding to h ii . Next, we carry out the Fourier transform of (D1) by introducing a crystal momentum k that lives in the reduced 1st Brillouin zone k x ∈ (−π/l x , π/l x ), k y ∈ (−π/l y , π/l y ) corresponding to the order parameter's unit-cell. We must formally join l x l y Nambu spinors (D1) into a single super-spinor, whose Fourier trans-formΨ k corresponds to a Bloch state. The Fouriertransformed Hamiltonian is: h is constructed as a block-matrix with blocks h ii (diagonal) and h ij (off-diagonal) for every two sites ij in the same unit-cell. c µ , µ ∈ {x, y} are the k-dependent blockmatrices made from h ij that connect the sites belonging to different unit-cells.
The full spectrum of the BdG Hamiltonian can now be obtained by diagonalizing the 8l x l y × 8l x l y matrix for any k in (D2). Since the Nambu representation had to be doubled for triplet pairing, the spectrum has an unphysical double degeneracy for every state. Only one of the two formally degenerate eigenstates is physical.

Category 2
A category-2 order parameter is generally a pairdensity wave given by: where η ± = (η ↑± , η 0± , η ↓± ) are three-component complex spinors and Q is a commensurate wavevector in the 1st Brillouin zone (1BZ) of the crystal lattice. In order to construct a BdG representation of the Hamiltonian, we must first identify the smallest positive integer q for which qQ is a reciprocal lattice vector; q specifies the size of a unit-cell in lattice sites for this periodic order parameter. Triplet pairing will mix fermion states at momenta separated by ±Q. In order to capture this, we define the Nambu spinor block in momentum space: and divide the 1BZ zone into q equal segments so that p = nQ + k with n = 0, . . . , q − 1. The residual wavevector k should live inside a reduced Brillouin zone whose area is q times smaller than that of 1BZ. However, the shape of the reduced zone is restricted by Q in somewhat complicated ways. It is more practical to allow k to vary inside the full 1BZ and then divide any 1BZ integrals (e.g. in the calculation of the ground state energy) by q to undo the overcounting of the BdG Hamiltonian eigenvalues. The BdG Hamiltonian is again the Fourier transform of all fermion terms in (2) and (8): +h + δ n,n ′ +1 + h − δ n,n ′ −1 Ψ n ′ Q+k , where: with E p = −2t cos(p x ) + cos(p y ) cos(a) − µ (D5) V p = 2t sin(p y ) + i sin(p x ) sin(a) .
We used the properties E −p = E p and V −p = −V p to construct the Hamiltonian. Once (D4) is diagonalized also with respect to the index n, one obtains the full set of Bloch states for the PDW order parameter.

Category 3
The category-3 order parameter is a plane wave: which is a subset of the category-2 order parameters (D3) with η − = 0 for commensurate Q. However, since the triplet now lives at a single wavevector, we can diagonalize the BdG Hamiltonian even when Q is incommensurate. To that end, we define the Nambu spinor in momentum space as: