Vortices and composite order in $\mathrm{SU}(N)$ theories coupled to Abelian gauge field

We consider $\mathrm{SU}(N)$-symmetric Ginzburg-Landau models coupled to non-compact Abelian gauge field focusing on the case $N>2$ at finite temperature. We show that, at least for sufficiently large gauge-field coupling constants, these models have two phase transitions. The intermediate phase between the symmetric and low-temperature phases is a state with composite neutral order and no Meissner effect. In this neutral phase the system spontaneously breaks only the symmetry associated with phase differences and density differences between components. For $N>2$, in contrast to the $SU(2)$ case, the neutral state cannot be mapped onto an $\mathrm{O}(M)$ model. We term this state ${\mathbb{C}{P}}^{N-1}$-neutral phase. We also show that while $\mathrm{SU}(N)$-symmetric Ginzburg-Landau models are not superconductors or superfluids in the usual sense, their state in external field at sufficiently low temperature is a vortex lattice.


I. INTRODUCTION
Phase transitions in U(1) Abelian gauge theories, such as superconductors, is an old problem which is fairly well understood in the single-component case. 1,2 The phase diagrams of multicomponent superconductors beyond mean-field approximation are more complex and less studied. Multicomponent gauge theories with Abelian gauge fields appear in the context of multicomponent electronic and nuclear superconductivity (see e.g. Refs. [3][4][5][6], and as effective field theories for other quantum systems. [7][8][9][10][11] The principal difference is that multicomponent models feature additional phases in which there is order only in various products of the complex fields, whilst individual phases are disordered. The most investigated cases are two-and threedimensional U (1) × U (1) superconductors described by two complex fields ψ 1,2 . They have, under certain conditions, a phase transition to a "paired" state where order remains only in the products of the fields such as ψ 1 ψ 2 = 0 or ψ 1 ψ * 2 = 0 while individual phases are disordered: ψ i = 0. 3,8,[12][13][14][15] The new states represent (i) neutral superfluids, where only counter-flow of charged components are allowed (termed metallic superfluid) and (ii) charge-4e superconductors. Analogous states of composite order occur for different microscopic reasons in strongly interacting superfluid mixtures on a lattice [16][17][18] and in pair-density-wave superconductors. [19][20][21] The origin of these transitions and additional phases is the fact that in these systems the composite topological defects, i.e. bound states of vortices in different components, have lower energy than the simplest vortices in the individual fields ψ 1,2 . The entropically driven proliferation of composite vortices disorders the phases of the individual fields but does not restore the symmetry completely: the system retains order in products of the fields up to a higher temperature where the elementary toplogical defects proliferate (for a discussion of the principle see Ref. 22).
Besides representing new types of states, paired phases have attracted interest in models connected with the discussion of the order of the phase transitions in the so-called deconfined quantum criticality problem. [8][9][10][11]23 Phases with composite order have also been found in SU(2) models. 9,10,24,25 A theory for an SU(2) double of complex fields coupled to Abelian gauge field can be thought of as a model with U(1) local × O(3) global symmetry. 26,27 The situation is more subtle in the SU(2) case than in the U(1) × U(1) case since it is believed that there are no energetically stable composite vortex-like topological defects in type-2 SU(2) gauge theories. 28 Nonetheless, for sufficiently strong coupling constant it has been shown that such SU(2) models have two transitions: at lower temperature the Meissner effect disappears but the system retains broken global O(3) symmetry which is restored only at a higher temperature. Currently there is no known duality argument relating phase transitions to statistical mechanics of topological defects in SU(N ) cases of the kind that exist for U(1)-systems. 1,2 However, it should be noted that even the lack of energetic stability of the composite vortices in the ground state does not necessarily exclude entropically excited such vortices driving the phase transitions. In that respect an instructive example is the magnetic response of SU(2) gauge theory: 29 while strictly speaking SU(2) gauge theory is not a superconductor [i.e. does not have a conserved U(1) topological invariant 30 ] and individual vortices carrying one flux quantum are not stable in the type 2 regime, 28 nonetheless in an external filed the system forms a vortex lattice 29 stabilized by intervortex forces. Likewise vortices can exist under certain conditions in neutral SU(2)systems. 31 The above raises the question of what the properties are of topological excitations and phase diagrams in SU(N ) gauge theories with N > 2. Such higher-N theories have attracted recent interest in the context of effective field theory for the quantium J-Q Heisenberg model, [32][33][34] which catalyzed renewed interest in SU(N ) gauge theories [35][36][37][38] . Another motivation to study these states is the prospect of the creation of artificial gauage fields in ultracold atoms, which can open up new possibilities of studying these generalizations of superconducting states.
The paper is structured as follows: First we present the models that we consider and the Monte Carlo simulation methods that we use. Thereafter we present our work on phase diagrams, and finally we present our work on response to external magnetic field.

II. MODELS WITH FIXED TOTAL DENSITY
We consider SU(N )-symmetric Ginzburg-Landau models in three spatial dimensions, given by the Hamiltonian density Here A is the magnetic vector potential and the ψ i = |ψ i |e iφi are matter fields corresponding to the superconducting components. The amplitudes are subjected to the constraint that the total superconducting density i |ψ i | 2 = 1, but are otherwise allowed to fluctuate. In order to elucidate the properties of the model (1) we note that it can be rewritten as where φ ij = φ j − φ i and j is the density of charged supercurrent: Here the first line gives the the terms associated with electrically charged currents and magnetic field energy.
On the second line we list those terms that give the energy from electrically neutral currents consisting of counterflows of charged components and gradients of relative density. Note that, quite generically, when all components have non-zero density, only the vortices that have winding in the phase of each component (composite vortices) will have finite energy per unit length 39 , as do vortices in ordinary single-component superconductors. These are therefore the energetically cheapest vortices that can be thermally excited. However, if such objects proliferate, they cannot restore symmetries associated with the phase and density differences between components. Below we investigate numerically the possibility of this scenario.

III. MONTE CARLO SIMULATION METHODS
We discretize the model (1) on a three-dimensional simple cubic lattice with L 3 sites and lattice constant a = 1. The discretized model is given by the Hamiltonian density is a lattice curl, is a gauge-invariant phase difference, k and l signify coordinate directions, and k is a vector pointing from a lattice site to the next site in the k-direction. Periodic boundary conditions are applied in all three spatial directions. The thermal probability distribution for states of the system at inverse temperature β is given by the Boltzmann weight We generate representative samples from these thermal distributions using Monte Carlo simulation. The simulations are performed using the Metropolis-Hastings algorithm with local updates of each of the degrees of freedom. The sizes of the local updates of the vector potential are adjusted during the initial part of the equilibration in order to make the acceptance probability 0.5. Each phase is updated separately by simply selecting a new value for the phase with uniform distribution on the range from −π to π. The amplitude degrees of freedom on a given site are updated together by randomly choosing a new density distribution. Care must be taken to use the correct measure when updating the amplitudes subject to the constraint of fixed total density.
For the determination of phase diagrams, we also use parallel-tempering swaps between systems with neighboring temperatures; typically one set of swaps is proposed every 24 sweeps. The simulated temperatures are adjusted within a fixed interval in order to make the acceptance ratios for parallel-tempering swaps equal for all pairs of neighboring temperatures. We use reweighting between simulated temperatures in order to improve our data (more on this below). For the determination of lowtemperature vortex configurations we use simulated annealing.
Equilibration is checked by comparing results obtained using the first and second halves of the data gathered after equilibration, and by comparing inverse-temperature derivatives obtained using finite differences and statistical estimators. Errors are determined by bootstrapping, and error bars correspond to one standard error.

IV. PHASE DIAGRAMS
In this section we consider the phase diagrams, in terms of electric charge q and inverse temperature β, for the cases of two, three and four components. These phase diagrams involve two types of order: superconducting order and CP N −1 order. We first describe the methods we use to locate phase transitions of each type, and then present and discuss our results on phase diagrams.

A. Locating superconducting transitions
While systems that do not have local U(1) symmetry are not superconductors in the usual sense [i.e. do not have a conserved U(1) topological invariant 22 ], they can exhibit magnetic field screening (the Meissner effect). Correspondingly we define a superconducting phase transition as a transition where the Meissner effect disappears. In order to locate superconducting transitions, we use the dual stiffness 9,25,40 where µνλ is the Levi-Civita symbol, ∆ ν is a difference operator and · is a thermal expectation value. Specifically, we measure the dual stiffness in the z-direction evaluated at the smallest relevant wave vector in the xdirection q x min = (2π/L, 0, 0), i.e. ρ z (q x min ); we denote this quantity simply as ρ. The quantity ρ will (in the thermodynamic limit) be zero in the Meissner state in which fluctuations of the magnetic field are suppressed, and non-zero in the normal phase. In the sense that it is zero in the low-temperature phase and non-zero in the high-temperature phase, it is thus a dual order parameter.
The quantity ρ is expected to scale as the inverse system size 1/L at the critical point of a continuous superconducting transition. Consequently, Lρ is a universal quantity, the finite-size crossings of which (extrapolated to the thermodynamic limit) give the critical temperature of a superconducting transition.

B. Locating neutral transitions
The model has a neutral sector that we call CP N −1 neutral sector. Again, because there is no conserved U(1) topological invariant the order in that sector is not strictly speaking superfluid, in contrast to the composite order in U(1) N models. 12,22 In order to locate neutral transitions we use the helicity modulus, which measures the free-energy cost of an infinitesimal phase twist.
Consider imposing a twist in a certain linear combination i a i φ i of the phases, i.e. replacing the the phase φ i by The helicity modulus for this phase combination is then where F is the free energy. Using the fundamental relations one can derive an expression for the second derivative of the free energy F in terms of the first and second derivatives of the Hamiltonian H with respect to the phasetwist parameter δ: In deriving the above, we have used that ∂H/∂δ = 0. Looking at the Hamiltonian (4) it is clear that any first derivative ∂H/∂δ will be a linear combination of sine functions, and any second derivative ∂H 2 /∂ 2 δ will be a linear combination of cosine functions. In both cases, the coefficients will depend on the a i . Because of the square in the second term in (12), the helicity modulus for a given linear combination of phases will depend not only on the corresponding single-phase helicity moduli, but also on certain cross terms that involve pairs of components. 15,41 In our case the components are all equivalent, and consequently each of the single-phase helicity moduli are equal, as are each of the cross-term helicity moduli. Furthermore, because the phase sum couples to the vector potential, its helicity modulus is necessarily zero. Form this is follows that the ratio between the value of a single-phase helicity modulus and the value of a cross-term helicity modulus is given by simple combinatorics, and that all helicity moduli are proportional to each other (the phase-sum helicity modulus being the special case where the constant of proportionality is zero).
Because of the aforementioned proportionality of all helicity moduli, we could in principle equivalently use any helicity modulus (except that for the phase sum) to locate neutral transitions. We choose to use the average of all single-phase moduli and the negatives of all cross-term moduli, as this would apear to minimize statistical error. We denote this quantity, when rescaled to equal a singlephase modulus, simple as Υ. In order to improve statistical precision, we determine a helicity modulus for a given temperature as as weighted average of the modulus determined from the data for that temperature and moduli obtained by reweighting from neighboring temperatures. We do this as follows: When reweighting from the nearest n temperatures on either side, the unreweighted modulus has relative weight n + 1 and the relative weights of the other moduli decrease by 1 for each step away from the temperature to which reweigting is performed. Since the true value of the phase-sum helicity modulus is zero, the average over all temperatures of the squared phase-sum helicity modulus is a measure of the statistical error in the helicity moduli. The number n is chosen to be the first natural number for which n + 1 gives a larger error. However, we limit n to a maximal value of n = 8 in order to limit the computational power required. Also, in Figure 1. Phase diagrams for N = 2 (blue, lower diagram) N = 3 (red, middle diagram) and N = 4 (green, upper diagram). The phase diagrams show that for high enough electric charge q there are new CP N −1 -neutral phases in which order is retained in the phase differences and relative-density degrees of freedom whilst superconducting order is absent. Errors are smaller than symbol sizes, and lines are a guide to the eye. some cases we measure helicity moduli in each coordinate direction and average them in order to further improve statistical precision. All this is motivated by the fact that for the systems studied here helicity moduli tend to have large statistical errors compared to other quantities that we study.
At the critical point of a continuous neutral transition, the helicity modulus is expected to scale as 1/L, so that LΥ is a universal quantity. We use finite-size crossings of LΥ, extrapolated to the thermodynamic limit, in order to locate neutral transitions.

C. Results
Phase diagrams for the two-, three-and fourcomponent cases are shown in Fig. 1. Our main conclusion is that the systems have split transitions and phases with composite order for N > 2. This occurs for high values of the coupling constant; if the temperature is increased the system first undergoes a phase transition where the Meissner effect disappears.
The remaining broken symmetry in the system is described by the following neutral terms in the original model: The new physics associated with these phases, in contrast to the U(1) N and SU(2) cases, is that these phases cannot be mapped onto O(N ) models. Instead these are CP N −1 -neutral phases where there is order associated with the spontaneous breaking of symmetry only in the phase differences and relative densities between components. For the highest charge simulated (q = 7) the shift in temperature for the superconducting phase transition to the completely ordered low-temperature state is much smaller than the corresponding shifts for the direct transitions at lower charges. The reason for this is that it is integer-flux excitations that govern the charged transitions, and increased charge makes these objects more tightly bound composites of fractional flux excitations [cf. the U (1) N case 13,39,42 ].
We now present more detailed results for some representative examples of the phase transitions in the aforementioned phase diagrams. We choose to consider the case N = 3, and the charges q = 3 for which there is a single transition and q = 6 for which there are two separate transitions. In Fig. 2 we show the helicity moduli, dual stiffnesses and heat capacities for q = 3 and system sizes in the range L = 12 − 48. Within our numerical uncertainty, the two transitions occur at the same temperature. Also, there are clear peaks in heat capacity that increase with system size. In Fig. 3 we show the helicity moduli and heat capacities for the neutral transition for q = 6 with system sizes in the range L = 12 − 40, and in Fig. 4 we show the heat capacities and dual stifnesses for the charged transition in systems with the same charge and sizes. The two transitions are clearly separated in temperature. Again, the transitions are associated with clear peaks in heat capacity that become more pronounced with increasing system size.

V. MAGNETIC RESPONSE OF SU (N ) SYSTEMS
The CP N −1 -composite neutral order demonstrated in the previous section is only possible because of the existence of integer-flux excitations. In a U(1) N case these are vortices where each complex component has 2π winding around a shared core. 39 As mentioned above, in contrast in the SU(N ) case the integer flux vortices are not energetically stable. 28 In certain cases it is possible to stabilize these objects by breaking symmetry explicitly to U (1)×Z 2 ; 43,44 then vortices are characterized not only by an overall phase winding but also by a CP N −1 skyrmionic topological charge when the cores of constituent fractional vortices do not coincide in space.
In this section we study the behavior of the models we consider in external magnetic field and at finite temperature. In particular we determine low-temperature vortex configurations. We find that vortices can be stabilized by the presence of external magnetic field.
This section consists of two parts: In the first part we describe how we implement external magnetic field and measure magnetic response numerically. In the second part we present and discuss our results involving external magnetic field.  16,20,24,32,40,48. Note that the scale for the inverse-temperature axis is different for the heat-capacity plot than for the other two plots; the crossings show much smaller finite-size corrections than the heat-capacity peaks. Errors are indicated by (narrow) shaded error regions.

A. Numerical implementation
For the purpose of implementing external magnetic field, we decompose the vector potential into a sum of two terms: A(r) = A 0 (r) + A 1 (r). Of these two terms, the first corresponds to a uniform magnetic field in the z direction. This term, which is written in the Landau gauge, is held constant: A 0 (r) = (0, 2πxf, 0). The second term, A 1 (r), is updated by Metropolis-Hastings updates and thus gives thermal fluctuations. Periodic boundary conditions are imposed on A 1 (r), whence the contribution to total magnetic flux from this term is zero. Thus the total  flux is equal to that given by the constant term A 0 (r). Given that periodic boundary conditions are used, the value of f must be such that hqLf is an integer.
For the purpose of describing the response of the system to external magnetic field, we consider two quantities: vorticities and magnetic flux density. For both quantities, we average over the z direction, which is the direction of applied field; this gives 2D images. In addition to these real-space images, we measure the absolute values of the Fourier transforms of these, i.e. structure factors. For all four types of quantity (direct images and structure factors of vorticities and magnetic field) we measure thermal averages. We remove the zerowavevector component of the structure factors for clarity, and normalize the remaining components to the zerowavevector component.
We now describe how vorticity is defined on a lattice. For a phase field defined on continuous space, it is clear how to determine if there is a vortex at a certain point: if so, the phase winds around this point. For a phase field defined only on the discrete points of a lattice, it is less clear what is means for there to be a vortex. The conventional way of counting the number of vortices on a given plaquette is this: The net vorticity of the plaquette is obtained by adding contributions from each of the links of the plaquette. For each link, consider the gaugeinvariant phase difference χ i,k (r). The contribution from a given link is the number of multiples of 2π that must be added to χ i,k (r) in order to bring it into the primary interval (−π, π]; this can be positive or negative.
Note that on a continuum the circulation integral of the gradient of phase is quantized, and that this does not depend on the vector potential. However, for small enough integrals around the vortex core the inclusion of the vector potential makes negligible difference, and thus the concept of vorticity used here is consistent with the usual concept on a continuum.
In the limit of zero charge q → 0 the magnetic field becomes completely uniform. We implement this numerically by choosing an initial condition for which the fluctuating part of the vector potential is zero, and then simply not updating the vector potential.

B. Vortex structure in external field
Vortex structures in neutral SU(2) systems under rotation have been considered. 31 This regime is equivalent to the limit of infinite magnetic field penetration length or q → 0 in the system under consideration. Vortex patterns in external magnetic field in the two-component zero-charge case are shown in Figs. 5-6. These vortex patterns are consistent with those found previously 31 in a model without a hard total-density constraint. Fig. 7 shows how these types of vortex state generalize to larger numbers of components. More precisely, it displays a q → 0 SU(3) superconductor in external field. In this three-component case the pattern is such that the Figure 5. Vortex patterns in external magnetic field in the two-component q = 0 case. Both real-space thermal averages (left) and thermally averaged structure factors (right) are shown. The vortices form stripe patterns in each component, and the stripes of the two components are interlaced. The inverse temperature is β = 1.5. The external field is the weakest that can be applied in the Landau gauge for L = 60; thus there are a total of 60 vortices in each component. Figure 6. Same disaplayed quantities and parameters as in Fig. 5, except that here the inverse temperature β = 2.5. The vortices form patterns that are of different types than in Fig. 5, and which are not symmetric between the two components.
vortices can be divided into groups of three. Such a pattern can be fully present only if the number of vortices is divisible by 3. This is why we consider the system size L = 60 (rather than, say, L = 64). Clearly, despite the instability of an individual vortex in the SU(3) case, under rotation, each component forms a regular vortex lattice, even in the presence of thermal fluctuations. Next we move to the effect of finite screening lengths. Vortex patterns in external magnetic field in the case of two components and finite screening length (q = 1) are shown in Figs. [8][9][10]. Since the charge is nonzero the � Figure 10. Same disaplayed quantities and parameters as in Fig. 8, except that here the inverse temperature β = 3.5.

VI. CONCLUSION
SU(N ) theories with non-compact Abelian gauge fields are of interest as effective theories of various condensed-matter systems; however, the most investigated case is that of SU(2). We have considered SU(N )-symmetric Ginzburg-Landau models at finite temperature with a focus on the case N > 2. We found that non-topological vortices play important roles in these systems. First we pointed out that at sufficiently large coupling constants the system has states with composite order which cannot be mapped onto S 1 or S 2 neutral systems like the previously discussed paired phases in U(1) × U(1) and SU(2) theories. The new phases display no Meissner effect but have a spontaneously broken symmetry associated with the relative phases and relative densities of the components. We term these phases CP N −1neutral states. These states result from proliferation of non-topological vortices that eliminate the Meissner effect, yet cannot disorder the CP N −1 -neutral sector of the model. Given the importance of these non-topological vortices, we study these systems in external field and at finite temperature. These systems are not superconductors or superfluids in the usual sense because of the lack of a U (1) topological invariant. Nonetheless we demonstrated that, even at finite temperature, well defined lattices of non-topological vortices are formed in externally applied magnetic field.