Half Quantum Vortices in Nematic Superconductor

Motivated by the superconductivity of $M_x$Bi$_2$Se$_3$, we study topological excitations in a nematic superconductor using Ginzburg-Landau theory. An isolated excitation at low field is shown to be either a distorted phase vortex or a tightly-bounded pair of half quantum vortices. Close to upper critical field $H_{c2}$, the vortex lattice is shown to be always hexagonal in the extreme type-II limit. Due to the different symmetries of the vortex lattice states, at least two phase transitions must take place when the external field is lowered from $H_{c2}$.

Superconductivity of M x Bi 2 Se 3 has attracted much attention [1]. Observations of in-plane anisotropy [2][3][4][5][6][7][8] lead to suggestion that the superconducting order parameter is multi-component and spontaneously breaks lattice rotational symmetry [9]. However, the leading theoretical candidate (odd-parity E u representation) is necessarily topological [10], yet experimental search for Majorana bound states turned out mixed results [11][12][13]. Besides, the current theoretical model calls for a preexisiting pinning field, and simple anisotropy measurements cannot rule out an anisotropic s-wave model as the alternative [14,15]. By theoretically investigating the topological excitations of a nematic superconductor, we seek to propose yet another way to prove (or reject) this "nematic hypothesis".
Half quantum vortex (HQV) is predicted to exist in a superfluid with a multi-component order parameter [16][17][18]. In contrast with the conventional phase vortex (PV), the phase and orientation of the order parameter each winds by π around an HQV, resulting in a singlevalued wavefunction topologically distinct from a PV. Evidence of HQV has been reported in the polar phase of 3 He [19], ultracold BEC [20][21][22][23], and superconducting Sr 2 RuO 4 [24]. Numerical solution of Ginzburg-Landau equation shows that tightly-bounded HQV pairs may ex-ist in M x Bi 2 Se 3 at low applied field [25]. Due to gauge screening and spin-orbit coupling unique to unconventional superconductors, the HQV physics is expected to be quite different from the counterpart in neutral superfluid.
In this letter, we take the phenomenological Ginzburg-Landau (GL) free energy of a two-component superconducting order parameter as our starting point, and analyze its solutions at both low field and near H c2 . We show that, at low field, an isolated topological excitation can be a PV or a bounded HQV pair. Near H c2 , a complete extension of the Abrikosov vortex lattice [26][27][28] to the multi-component case is presented, and we find a solution consisting of an array of defects with topological charges identical to low-field HQVs. Based on the difference in symmetry, we argue that at least two phase transitions are required at intermediate field to connect the two limits.
The Model. We assume the system is uniform along the c-direction of the lattice, effectively two-dimensional. The superconducting order parameter η i (indices i, j = x, y) forms a two-dimensional representation of the symmetry group (D 3d for M x Bi 2 Se 3 ). We write down the general quartic GL free energy density coupled to the gauge field A i : where p i = −i∂ i + A i . We work in the unith = (−e * /c) = 1. It is often useful to consider the alternative basis η ± = η x ± iη y . Let p ± = p x ± ip y likewise, and one has arXiv:2005.03366v1 [cond-mat.supr-con] 7 May 2020 where K ij stands for K i + K j , and δβ = 1/2 + β.
One can identify ± η up to a global phase shift of π: the orientation acts as a nematic director. In two dimensions, a nematic director is equivalently labeled by the direction normal to it. Consequently this theory enjoys a duality: F remains invariant under the replacement of In the absence of external gauge field, this theory has two distinct phases. These uniform solutions are of the form: We will focus on the nematic phase in this letter. Upon integration by part, it can be seen that the difference (K 2 − K 3 ) only couples to (∇ × A)( η) * × η. We expect this term to be generally small since it signifies particle-hole asymmetry [18,29]. We will set K 2 = K 3 throughout.
Zhitomirskii [27] investigated the behavior of this model. It will be helpful to define C = K 23 /2K 1 [41] Stability conditions are HQV pair at low field. Consider the nematic solution in (4). An HQV centered at the origin resembles: where φ is the azimuthal angle, and (n p , n o ) = (±1, ±1) marks its topological charge. In the ±-basis introduced earlier, (7) becomes A HQV can thus be viewed as a single vortex of one of the η ± components, in a non-zero background of the other. In a superconductor, only the phase fluctuation is screened by the gauge field. A configuration with nonzero total n o charge will have an energy divergent with system size. HQVs must therefore come in pairs with opposite n o charges. Does the pair reach equilibrium at a non-zero separation, or does it collapse to form one PV? To proceed, we consider the extreme London limit, where the magnitude | η| is everywhere uniform except at the point-like core of a topological defect.
When the HQV pair is separated by d λ the penetration depth, the gauge interaction is screened out, while the orientation gradient ∇θ still significantly deviates from zero in the region between the two cores. Reducing a evidently reduces the gradient energy. The HQV pair is attractive at this very large separation.
Next, we consider the intermediate case where the cores are still well separated, but the separation d < λ.
As A varies with length scale λ, it can be treated as approximately uniform. We will set A = 0 for this part. The effective free energy density to be considered is: where gives overall dimensions, g = C/(1 + C), andû = (cos θ, sin θ) is the orientation unit vector. The duality for this effective theory is The total free energy of the pair F pair can be split into three parts: the energy for two isolated HQV cores 2F core , the logarithmic interaction energy F log , and an dipole potential F dip that depends on the relative alignment of the pair to the background η at infinity. One expects F log to be purely logarithmic in d because (9) is scaleinvariant. As d → 0, F pair must match F PV , the core energy of an isolated PV.
The core energy F core is calculated[43] by seeking solution to the theory (9) of the form (7). Noting that δθ = δχ = 0 is an exact solution at g = 0, one may solve the GL equations order-by-order in g, and compute the perturbation series for the free energy: Short-and large-distance cutoffs ξ and Λ are required in the integration over space. The series contains only even powers of g because the free energy is self-dual.
To calculate F dip , consider a pair of HQVs centered at (±d/2, 0) respectively. See Fig 1. One can write: The angle between η at infinity and the pair is θ ∞ . The dipole energy turns out to be The self-dual requirement implies the absence of even powers of g. Similar calculation is performed for F PV , taking χ = φ + δχ and θ = 0 + δθ. The result is compared with 2F HQV , and the mismatch in logarithmic part is attributed to F log : This is an attraction: the pair continues to collapse at intermediate separation. The dipole potential (13) aligns the pair to be parallel (perpendicular) to the background η when C is positive (negative). Finally, if the cores repel at very short distance, the HQVs can avoid total collapse and stay as a tightlybounded pair. An HQV here is interpreted as a single vortex of η + or η − . The geometry remains identical to where the amplitude . The function f (r) represents the amplitude variation around the core starting out at zero and saturated when r > ξ. We assume d ξ. One must employ the GL free energy (2) directly. Only the terms that couple η + and η − will be dependent on d. Let us define two integrals: The d-dependent part of the free energy F int is then Both I 23 and I β are decreasing function of d. The I 23 term always favors separation, and the pair is aligned parallel (perpendicular) to the η at infinity if C is positive (negative). The I β term favors separation only if δβ > 0. If δβ < 0, the competition between the two terms determines whether the pair collapses into one single PV. The qualitative behavior is summarized in Fig 2. Vortex Lattice at high field. Near H c2 , the external magnetic field dominates, and there is no distinction between the nematic and chiral phases. We return to the free energy (2), as it is natural to consider η ± here.
[44]. The first step is to solve the linearized GL equation at applied field H = H c2 . We look for a periodic solution with lattice parameter a, b, and angle α between them. These parameters satisfy the constraint ab sin α = 2π/H c2 , and it suffices to consider distinct ρ = b/a and α. We adopt the gauge (A x , A y ) = (−H c2 y, 0), define the Landau level wavefunctions: and finally form the periodic ansatz: where the magnitude √ Ω is a minimization parameter. We consider the Gibb's free energy density Using the linear solution (19) as variational ansatz, the Gibbs free energy to be minimized is where f 4 is the part of (2) quartic in η ± , and h s = (∇× A − H). Angular brackets denote spatial average. The relation (21) is general, and yields exactly the Abrikosov formula when applied to an s-wave superconductor. At leading order in ∆H, the Maxwell equation ∇×h s = −4π (δF/δ A) can be integrated analytically to give an exact expression for h s . We obtain[45] To our knowledge, this crucial point has been missed in existing literature [30,31]. Lacking a key ingredient, previous authors took additional simplifying assumptions to derive their minimization targets. Our treatment, in contrast, exactly mirrors the original Abrikosov method for s-wave superconductor.
With analytic expressions of both f 4 and h s , R in (21) is numerically minimized with respect to ρ and α to determine the most favorable lattice structure. Apart from C and β, the quantity κ introduced in (6) is the only other dimensionless parameter that controls the physics. We identify three regimes of distinct behaviors across a range of C, β and κ.
Given the duality (3), only C ≥ 0 needs to be considered. At C = 0, ansatz (19) is effectively singlecomponent, and hexagonal lattice is always favored, independent of β and κ. For C > 0, we first consider the limit κ → 1/ √ 2. We find that near the minimal value β → −1, the hexagonal phase is stable for all C.
Consider increasing β. At β = β s ≈ −0.20, the square phase first appears at C = 1. Further increasing β, and the square phase spread to smaller C. The transition between the square and hexagonal phases is a sharp jump.
Further increasing β destabilizes the new square phase, too. For β > β i ≈ 0.05, again starting from C = 1 and spreading inward, the free energy landscape becomes extremely flat and develops multiple near-degenerate local minima (without apparently any symmetry) on the (ρ, α) plane. We are unable to resolve them with the available numerical precision. [46]. The minimum that is connected to the square lattice also drifts into ρ = 1, though α remains locked at π/2. We term this regime the irregular phase. The behavior is summarize in Fig 3. The model exhibits the same qualitative behavior at larger κ, though β s and β i are increasing functions of κ. At κ = 100, which is in the range relevant to M x Bi 2 Se 3 [32,33], we found β s ≈ 0.32 and β i ≈ 0.36. This means Fig 3 is largely irrelevant to the physics of M x Bi 2 Se 3 , which has β < 0 and must be in the hexagonal phase. However, we speculate that the irregular phase found here may be related to the reported complicated and disordered vortex arrangement in Sr 2 RuO 4 [34], which is described by a very similar GL theory with β > 0 and a moderate κ of 2.6 [35].
Assuming a hexagonal lattice, let us examine the structure of (19). The η − component is identical to the hexagonal solution of an s-wave superconductor, and contains exactly one zero per unit cell, each with winding number −1 (clockwise). In contrast, the η + component has three zeroes per unit cell. One of them coincides with the zero of η − , and carries winding number +1 (counterclockwise). The other two zeroes are located at the centers of the two equilateral triangles, respectively, and each has winding number −1. The total winding number per unit cell for the η + component still sums to −1. The relative phase between the η ± components is periodic. Considering both components, the total topological charge per unit cell is equal to that of one PV, or one HQV pair at low field. See Fig (4). Bridging the two limits. As the external field H is lowered, the high-field vortex lattice must evolve into an array of well-isolated topological excitations. These can be single PVs or tightly-bounded HQV pairs, as previously discussed.
Despite having the correct total topological charge, however, the η + component of the vortex lattice solution has too many zeroes (three) per unit cell. Annihilation must occur, leaving only one negative-winding zero. At least two phase transitions take place during the process [47].
Near H c2 , the vortex lattice enjoys full hexagonal symmetry. Near H c1 , however, the background superconducting order parameter defines a special nematic orientation, and the well-isolated topological excitations form a hexagonal array distorted along the nematic direction. If PVs are favored, the resultant array has C 2v symmetry. If HQV pairs are preferred, the dipolar nature of the pair further reduces the symmetry down to only one mirror reflection. We conclude that the hexagonal symmetry near H c2 must be broken in a phase transition as the field is lowered.
The locations of zeroes of the η + component are protected by the hexagonal symmetry near H c2 , and the symmetry-breaking transition allows the zeroes to drift. Eventually the positive-winding zero will meet and annihilate with a negative-winding one, and this annihilation is another phase transition.
Discussion and conclusion. In this letter we study the topological excitations of a nematic superconductor. At low applied field, an isolated excitation is either a PV or a tightly-bounded HQV pair. This is in agreement with existing numerical simulation [25] that observed tightlybounded HQV pairs at low field. Available images of vortices in Cu x Bi 2 Se 3 [36] shows only blurry, elongated defects. It would be interesting if individual HQVs can be visualized in experiments. For instance, local density of states probed by STM may show a twin-peak structure due to the HQV pair.
We construct the phase diagram for the vortex lattice state near H c2 . For the range of parameters relevant to M x Bi 2 Se 3 , the vortex lattice is always hexagonal. This agrees with the finding in [37]. The vortex lattice solution at H c2 developed in this letter is the full solution for an unconventional superconductor without extra simplifying assumptions. Our method can be adopted for other unconventional superconductors, and may prompt the reexaminations of prior results.
Irrespective of whether the M x Bi 2 Se 3 family of materials supports HQV pairs at low field, we predict at least two phase transitions between high-field and lowfield states. The existence of these transitions is a direct validation of an unconventional order parameter. [41] The value of C can be estimated from the anisotropic ratio of in-plane Hc2; see [38] and supplemental material of [7]. For Cu doping, we find C ≈ 0.6 using [5]. For Sr however, the reported ratios in [3,[6][7][8]

I. GL PARAMETERS AND LENGTH SCALES
In the s-wave case, the GL parameter κ s-wave plays a double role. On one hand, it is the ratio of magnetic penetration depth to the coherence length, i.e. κ s-wave = λ/ξ; on the other, it is the sole parameter that distinguishes between type-I and type-II behavior, and vortices are only allowed if κ s-wave > 1/ √ 2. For the present model, the two roles are no longer played by the same quantity.
The criterion for type-II superconductivity is the existence of an upper critical field H c2 higher than the thermodynamic critical field H c . For the present theory, H c2 was calculated in [1]. In our notation it reads: where we have introduceκ H c can be easily read off of the free energy (1) in the main text. It is different for chiral and nematic phases.
The requirement of H c2 > H c is captured by This κ parameter, however, is not directly accessible experimentally. The ratio of length scales λ/ξ is the quantity reported in experiments. The experimental estimations put this ratio typically at the order of 100 for M x Bi 2 Se 3 . We adopt the s-wave-inspired for estimating ξ, as is usually done in the literature. Penetration depth λ is extracted from the GL equation for A expanded around the uniform state solution. It is anisotropic in the nematic phase, since the uniform background of order parameter itself breaks the rotational symmetry.
We define the ratio of length scales For the nematic phase, the geometric mean of the two anisotropic penetration depth is used in the ratio. Assuming 1 + β is not extremely small, as far as order of magnitude is concerned, one can see that κ r ≈ κ for the entire range of −1/3 < C < 1 in the nematic phase. The above discussion informs our choice of range of κ in the vortex lattice calculation based on the experimentally reported κ r . arXiv:2005.03366v1 [cond-mat.supr-con] 7 May 2020

II. FREE ENERGY OF HQV PAIR AT INTERMEDIATE DISTANCE
We elaborate on the steps taken to obtain the equations (11), (13) and (14) in the main text. Given the effective free energy density we want to evaluate and compare the free energy cost various field configurations. Before going on to the analysis, we would like to remark on the difference between (S8) and the effective model for spin fluctuation commonly employed in the context of liquid 3 He. Spin-orbit coupling in liquid 3 He is negligible. Consequently the terms proportional to g in (S8), where the spin orientation of the Cooper pairû couples with the partial derivative, do not appear in the model for 3 He. On the other hand, the strong coupling Fermi liquid correction in 3 He leads to |∇θ| 2 and |∇θ| 2 having unequal coefficients.
The GL equations for θ and χ are These equations can be solved order-by-order in g, and the results plugged into (S8) and integrated over the entire space to obtain a naive perturbation series in g of the corresponding free energy. For one isolated HQV centered at the origin, the solution turns out to be where φ is the azimuthal angle and (n o , n p ) = (±1, ±1). This solution is consistent with the duality, which rotates θ HQV by π/2 on the left hand side, and sends g → −g and rotates θ ∞ by π/2 on the right hand side. All four types of HQVs turns out to have the same free energy, at least up to O(g 3 ). Similar to a PV in a conventional neutral superfluid, the energy is logarithmically diverging at both large and small distance, and we introduce the respective cutoffs Λ and ξ. At the end, the free energy for an isolated (neutral) HQV is with all the odd powers of g vanish due to the self-dual requirement as discussed in the main text. The short-distance cutoff ξ masks all the physics of the core region, where the magnitude of the order parameter | η| has significant variation and the effective free energy (S8) fails to be a good approximation. Along the same vein, Λ can be viewed as a finite size cutoff that is smaller than the electromagnetic screening length scale λ (and any other large distance scale naturally present in the full theory). All the large-distance physics gets masked by this one single cutoff parameter.
For an isolated PV, the field configuration is and the free energy is (S13) The logarithmic part of the free energy of an HQV pair should interpolate between F PV and 2F core as the pair separation d is varied. When d approaches core size ξ, the pair really has the appearance of one single P V . On the other hand, two well-separated HQVs should appear isolated, and here we take the criterion for isolation to be d ∼ Λ. These considerations lead to the identification of the logarithmic interaction (S14) The dipole energy F dip must be calculated using the field configuration of an HQV pair. The zeroth order solution is as in the main text. The O(C) correction to free energy comes solely from substituting θ 0 and χ 0 into (S8). Let u 0 = (cos θ 0 , sin θ 0 ), it is easy to see that (S16) We wish to draw the reader's attention to the form of this integral. It is proportional to the integral I 23 defined in (16) in the main text, if one sets the core function f (r) = 1, corresponds to the extreme London limit where the core shrinks to one point. The dipole energy turns out to be: (S17) Since cos(2θ ∞ ) changes sign as well under duality, the self-dual requirement forces all even powers of g to vanish.
III. VORTEX LATTICE NEAR Hc2 FOR UNCONVENTIONAL ORDER PARAMETER

A. Minimization Target
As we have mentioned in the main text, the original Abrikosov method for finding the most favorable vortex lattice solution at high field can be executed in the present case without further approximation. We briefly recount the method as applied to the present case here.
Let us denote the Gibbs free energy density as where F 2 is quadratic in η, and F 4 is quartic. The GL equations reads We will work in the Landau gauge and set A 0 = (−H c2 y, 0). Then one identifies H c2 by demanding that the linearized equation has solutions. Indeed, it has many degenerate solutions at this point. Let η 0 be any one of these degenerate solutions of (S21). When the external field H is lowered from H c2 , the massive degeneracy is broken by the quartic term, and each degenerate configuration evolves into a saddle point of the free energy. Such a saddle point solution can be perturbatively calculated: Here both Ω and A 1 are directly proportional to ∆H ≡ H − H c2 . These are required to calculate the Gibbs free energy G up to O(∆H 2 ), and can be solved by plugging (S22) into (S19) and (S20).
For the moment assume that (S20) can be integrated (we will show below explicitly how this can be done) to give Using this, the GL equation (S19) and (S20), and the fact that η 0 solves the linearized equation (S21), one may obtain the corresponding Gibbs free energy G: Angular bracket denotes spatial average. The true minimum of the free energy is found by minimizing (S24) with respect to the initial ansatz η 0 . Note that one can skip the actual evaluation of Ω when using (S24): it is canceled in the ratio. From this point on, it is more natural to switch to the η ± = η x ±iη y basis. The readers are directed to Zhitomirskii's treatment for the solution of linearized equation (S21), bearing in mind that our assumption of K 2 = K 3 means Zhitomirskii's D coefficient vanishes. The periodic solution was quoted in the main text and we shall repeat it here: The ratio ρ = b/a. This solution has the periodicity of a lattice with primitive lattice vectors that satisfy ab sin α = 2π/H c2 . It is worth noting that ψ j is a jth Landau level wavefunction, and the lowest Landau level ψ 0 is annihilated by p − . The spatial average F 4 can be numerically computed by first going to the Fourier representation of F 4 . It is the local variation of the magnetic flux h s that requires further attention. For the s-wave case, it was originally argue that h s ∝ |ψ 0 | 2 using special identities available only because the linearized solution is the lowest Landau level. For the present case, the ansatz (S25) contains both ψ 0 and ψ 2 , but the Landau level structure still allows us to integrate (S20) analytically.
From (S20) one may write down Integrating with respect to y, in principle, undoes the partial derivative on h s . To determine the constant of integration, we replace ψ j with the regularizedψ j , where the infinite sum over n in the original definition of ψ j is instead cut off at ±N for some N . The newψ j vanishes exponentially as y → ±∞, but remains a jth Landau level wavefunction. With the regularized ansatz, the system is completely normal at y → −∞, and h s must vanish there. We therefore choose to begin the integration from y = −∞. After the integration is done, the cutoff N is sent to infinity to restore the spatially extended periodic ansatz. Let x denotes the regularized version of j x , and the above discussion amounts to The exact form of the current j x = δF 2 /δA x is: Utilizing the commutator [p + , p − ] = 2H c2 , and the lowest Landau level property one may trade all occurrences of the ladder operators p ± in (S29) for some combination of partial derivative ∂ y . The integration with respect to y is now trivial. The y → −∞ part vanishes due to our regularization. And then the same trick is employed again to trade all ∂ y in the result for p + . Finally the cutoff N is allowed to go to infinity again. The end result is It is now straightforward to numerically evaluate h s 2 and h 2 s for any given ρ and α. When C = 0, ω → ∞, and the ansatz (19) in the main text becomes effectively single-component, containing η − only. In this limit the problem is formally identical to the s-wave case, and the vortex lattice is always hexagonal.

B. Helmholtz Free Energies
So far we have framed the discussion exclusively in terms of the Gibbs free energy, as we deem it conceptually more appropriate for the present case with an externally applied magnetic field. Historically, however, the procedure has been carried out in terms of the Helmholtz free energy. We will apply Legendre transformation and establish the equivalence between our result and the conventional lore.
Let B = ∇× A be the average magnetic flux, and ∆B = B − H c2 will be the new expansion parameter. Recall that our Gibbs free energy density was defined as with the addition of the constant H 2 /8π. To get the usual Helmholtz free energy F , one needs Definition (S32) also implies that ∂ G /∂H = −(B − H)/4π. Using (S24) for G , one may write down: This implies ∆H and ∆B are of the same order. One can now identify from (S24) and (S34) which turns (S33) into Finally, by using (S34) again, ∆H can be eliminated in favor of ∆B. The end result is in agreement with the existing literature. For the s-wave case with a scalar order parameter ψ, integrating the current yields simply h s ∝ |ψ| 2 , while f n ∝ |ψ| n for n = 2, 4. One instantly recovers the standard s-wave result

IV. VORTEX LATTICE NEAR Hc1
In this section, we explicitly argue that the topological excitations form a distorted hexagonal array near H c1 . This is regardless of whether these excitations are single PVs or HQV pairs. Let us first revisit the ansatz (S12) and (S15) for isolated PV and HQV pair, respectively. It is clear that in either case the orientation angle θ is roughly a constant when viewed from a distance R much greater than core size ξ: where the ξ/R term is the dipolar distortion that simply vanishes for PV. We have taken the liberty to re-align the coordinate axes so that θ ∞ = 0. Consequently, the unit vectorû 0 ≈x, and ∇θ ≈ 0. Now we re-derive the effective free energy in the extreme London limit, with the gauge sector included, using the above approximation. The free energy reads: where M is the appropriate constant coefficient. By rescaling one brings the free energy into an effectively isotropic form It is a well known fact that, for this isotropic theory, single PVs at finite density arrange themselves into a regular hexagonal array. Assuming g > 0 (and consequently C > 0), to account for the anisotropic scaling (S41), the hexagonal lattice and the vortices themselves are elongated in the direction parallel to the background η, and compressed in the orthogonal direction. When g < 0, the distortion is the exact opposite, and this is consistent with the duality where g is mapped to −g if η is rotated by 90 degrees. We expect hexagonal vortex lattice and the background η to both be aligned to some high-symmetry direction of the underlying trigonal lattice of M x Bi 2 Se 3 .
We shall assume g > 0 in all the subsequent discussion. The solution (S12) indicates that the shape of an isolated PV receives a quadrupolar distortion due to the background. More generally, the background (even allowing for periodic modulation due to the vortex lattice) is always invariant under π spatial rotation followed by π phase shift. Therefore a distorted PV retains the symmetry of an ellipse, with its semi-major axis aligned with the background nematic direction. The symmetry point group of this low-field PV state is C 2v , containing a two-fold axis going through the ellipse's center, and two mirror planes perpendicular to the ellipse's axes. We will denote this as the P V phase. See Fig S1 for an illustration.
In the discussion of lattice symmetry, any mirror reflection is implicitly followed by time-reversal to compensate for the reversal of winding. Furthermore, the two-fold rotation and one of the reflection are defined with a global phase shift of π to compensate for the reversal of the nematic orientation angle θ.
Suppose the HQV pair is preferred instead. The pair is dipolar in nature, but at a distance R from its center, the dipole behavior is suppressed by the small ratio ξ/R. At low density, these pairs should likewise arrange themselves into a distorted hexagonal array. But the dipole nonetheless breaks the two-fold rotation and the mirror plane not parallel to it. At the end, there is only one mirror plane left in the symmetry point group. This phase is denoted as HQV .

V. POSSIBLE β-H PHASE DIAGRAM
We extend the discussion on how the vortex lattice near H c2 can be connected to the low-field limit in the main text. When going from high field to low field, we only know for certain the starting point (H c2 vortex lattice) and two different end points (P V or HQV phase), and every detail in the middle is only speculation. Thus we refrain from including the discussion in the main text. Our best guess is guided by two principles: that the sequence of events should be as simple as possible, and that at all time the lattice symmetry should reflect the actual symmetry of the array.
Let T denote the high-field triangular/hexagonal phase near H c2 . As described in the main text, in the T phase the zeros of the η − component all have negative winding and form a regular hexagonal array. The eta + component has both positive-and negative-winding zeros. The positive-winding zeros coincides with the zeros of η − , and the negative-winding zeros sit at the centers of the triangles. We introduce two more intermediate phases D 1 and D 2 . In D 1 , pairs of negative-winding zeroes of η + gravitate toward the hexagonal lattice point between them, breaking the lattice symmetry down to C 2v . In D 2 , in addition to the converging of the negative-winding zero pairs, the  The D 1 phase. The unfilled circles drift away, breaking the hexagonal symmetry, and the distortion of the background hexagon is also indicated.
The unbroken symmetry group is C 2v . (c) The D 2 phase. The filled circles also move, and the only remaining symmetry is the x → −x reflection.
positive-winding zero originally resting on the lattice site also drifts away toward one of the two negative-winding zeros, breaking all symmetry except the mirror plane parallel to the trio of zeros. The phases are illustrated in Fig  S2. It is clear that D 1 evolves naturally into P V when the three zeroes merge at the lattice site, and that D 2 turns into HQV when the pair of approaching zeros annihilate, leaving one negative-winding zero for η − and η + each to form a dipole. Finally, even at low field the HQV phase should cease to exist when β changes sign, because it only makes sense in a nematic background.
When β is positive, the uniform chiral superconducting state breaks the time-reversal symmetry. Near H c1 , it has been shown [2] that the isolated "crescent vortex" solution is dipole-like and has only a reflection symmetry. In our proposed D 2 phase, the positive-winding zero of η + and the zero of η − can pair up to form a dipole-like structure, while the remaining zeros of η + can broaden out so that the interior of the unit cell between lattice points sees a background strongly dominated by the chiral η − . In other words, our D 2 phase can evolve smoothly into a crescent vortex lattice in a chiral background on the β > 0 side. A plausible phase diagram on the β − H plane is given in S3.
The phase boundary T -D 1 and T -D 2 sees the breaking of three-fold rotational symmetry, among other things. This is most likely a first order transition [3]. There is no reduction in symmetry across D 1 -P V and D 2 -HQV phase boundary, so the transition too can only be first order. The breaking of two-fold rotation across the D 1 -D 2 and P V -HQV boundary means the transition may be second order.