{'fraction_non_alphanumeric': 0.03863227794694941, 'fraction_numerical': 0.013586723899671032, 'mean_word_length': 4.238101338799565, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Crystallography of Color Superconductivity

We develop the Ginzburg-Landau approach to comparing different possible crystal structures for the crystalline color superconducting phase of QCD, the QCD incarnation of the Larkin-Ovchinnikov-Fulde-Ferrell phase. In this phase, quarks of different flavor with differing Fermi momenta form Cooper pairs with nonzero total momentum, yielding a condensate that varies in space like a sum of plane waves. We work at zero temperature, as is relevant for compact star physics. The Ginzburg-Landau approach predicts a strong first-order phase transition (as a function of the chemical potential difference between quarks) and for this reason is not under quantitative control. Nevertheless, by organizing the comparison between different possible arrangements of plane waves (i.e. different crystal structures) it provides considerable qualitative insight into what makes a crystal structure favorable. Together, the qualitative insights and the quantitative, but not controlled, calculations make a compelling case that the favored pairing pattern yields a condensate which is a sum of eight plane waves forming a face-centered cubic structure. They also predict that the phase is quite robust, with gaps comparable in magnitude to the BCS gap that would form if the Fermi momenta were degenerate. These predictions may be tested in ultracold gases made of fermionic atoms. In a QCD context, our results lay the foundation for a calculation of vortex pinning in a crystalline color superconductor, and thus for the analysis of pulsar glitches that may originate within the core of a compact star.


Introduction
Since quarks that are antisymmetric in color attract, cold dense quark matter is unstable to the formation of a condensate of Cooper pairs, making it a color superconductor [1,2,3,4,5]. At asymptotic densities, the ground state of QCD with quarks of three flavors (u, d and s) is expected to be the color-flavor locked (CFL) phase [6,5]. This phase features a condensate of Cooper pairs of quarks that includes ud, us, and ds pairs. Quarks of all colors and all flavors participate in the pairing, and all excitations with quark quantum numbers are gapped. As in any BCS state, the Cooper pairing in the CFL state pairs quarks whose momenta are equal in magnitude and opposite in direction, and pairing is strongest between pairs of quarks whose momenta are both near their respective Fermi surfaces.
Pairing persists even in the face of a stress (such as a chemical potential difference or a mass difference) that seeks to push the quark Fermi surfaces apart, although a stress that is too strong will ultimately disrupt BCS pairing. The CFL phase is the ground state for real QCD, assumed to be in equilibrium with respect to the weak interactions, as long as the density is high enough. Now, imagine decreasing the quark number chemical potential µ from asymptotically large values. The quark matter at first remains color-flavor locked, although the CFL condensate may rotate in flavor space as terms of order m 4 s in the free energy become important [7]. Colorflavor locking is maintained until an "unlocking transition", which must be first order [8,9], occurs when [8,9,10,11,12] µ ≈ m 2 s /4∆ 0 . (1.1) In this expression, ∆ 0 is the BCS pairing gap, estimated in both models and asymptotic analyses to be of order tens to 100 MeV [5], and m s is the strange quark mass parameter. Note that m s includes the contribution from any ss condensate induced by the nonzero current strange quark mass, making it a densitydependent effective mass. At densities that may occur at the center of compact stars, corresponding to µ ∼ 400 − 500 MeV, m s is certainly significantly larger than the current quark mass, and its value is not well known. In fact, m s decreases discontinuously at the unlocking transition [13]. Thus, the criterion (1.1) can only be used as a rough guide to the location of the unlocking transition in nature [13,12]. Given this quantitative uncertainty, there remain two logical possibilities for what happens as a function of decreasing µ. One possibility is a first-order phase transition directly from color-flavor locked quark matter to hadronic matter, as explored in Ref. [11]. The second possibility is an unlocking transition, followed only at a lower µ by a transition to hadronic matter. We assume the second possibility here, and explore its consequences. One may think that ud BCS pairing could persist below the unlocking transition. However, this does not occur in electrically neutral bulk matter [12]. The unlocking transition is a transition from the CFL phase to "unpaired" quark matter in which electrical neutrality is enforced by an electrostatic potential µ e ∼ m 2 s /4µ, to lowest order in m s /µ, meaning that the Fermi momenta are related (to this order) by If this "unpaired" quark matter exists in some window of µ between hadronic matter and CFL quark matter, there will be pairing in this window: all that Ref. [12] shows is that there will be no BCS pairing between quarks of different flavors. One possible pattern of pairing is the formation of uu , dd and ss condensates. These must be either J = 1 or symmetric in color, whereas the QCD interaction favors the formation of color antisymmetric J = 0 pairs. The gaps in these phases may be as large as of order 1 MeV [14], or could be much smaller [3]. Another possibility for pairing in the "unpaired" quark matter is crystalline color superconductivity [15,16,17,18,19,20], which involves pairing between quarks whose momenta do not add to zero, as first considered in a terrestrial condensed matter physics context by Larkin, Ovchinnikov, Fulde and Ferrell (LOFF) [21,22]. Unpaired quark matter with Fermi momenta (1.2) is susceptible to the formation of a crystalline color superconducting condensate constructed from pairs of quarks that are antisymmetric in color and flavor, where both quarks have momenta near their respective unpaired Fermi surfaces. We shall argue in this paper that the crystalline color superconducting phase is more robust than previously thought, with gaps comparable to ∆ 0 . In fact, crystalline color superconducting pairing lowers the free energy so much that its pairing energy should be taken into account in the comparison (1.1). This comparison should be redone as a competition between the crystalline color superconducting and CFL phases: crystalline color superconductivity in the "unpaired" quark matter is sufficiently robust that it will delay the transition to the CFL phase to a significantly higher density than that in (1.1).
The idea behind crystalline color superconductivity is that once CFL pairing is disrupted, leaving quarks with differing Fermi momenta that are unable to participate in BCS pairing, it is natural to ask whether there is some generalization of the pairing ansatz in which pairing between two species of quarks persists even once their Fermi momenta differ. It may be favorable for quarks with differing Fermi momenta to form pairs whose momenta are not equal in magnitude and opposite in sign [21,22,15]. This generalization of the pairing ansatz (beyond BCS ansätze in which only quarks with momenta which add to zero pair) is favored because it gives rise to a region of phase space where both of the quarks in a pair are close to their respective Fermi surfaces, and such pairs can be created at low cost in free energy. Condensates of this sort spontaneously break translational and rotational invariance, leading to gaps that vary in a crystalline pattern.
Crystalline color superconductivity may occur within compact stars. With this context in mind, it is appropriate for us to work at zero temperature because compact stars that are more than a few minutes old are several orders of magnitude colder than the tens-of-MeV-scale critical temperatures for BCS or crystalline color superconductivity. As a function of increasing depth in a compact star, µ increases, m s decreases, and ∆ 0 changes also. This means that in some shell within the quark matter core of a neutron star (or within a strange quark star), m 2 s /µ∆ 0 may lie within the appropriate window where crystalline color superconductivity is favored. Because this phase is a (crystalline) superfluid, it will be threaded with vortices in a rotating compact star. Because these rotational vortices may pinned in place by features of the crystal structure, such a shell may be a locus for glitch phenomena [15].
To date, crystalline color superconductivity has only been studied in simplified models with pairing between two quark species whose Fermi momenta are pushed apart by a chemical potential difference [15,16,17,18,20] or a mass difference [19]. We suspect that in reality, in three-flavor quark matter whose unpaired Fermi momenta are split as in (1.2), the color-flavor pattern of pairing in the crystalline phase will involve ud, us and ds pairs, just as in the CFL phase. Our focus here, though, is on the form of the crystal structure. That is, we wish to focus on the dependence of the crystalline condensate on position and momenta. As in previous work, we therefore simplify the color-flavor pattern to one involving massless u and d quarks only, with Fermi momenta split by introducing chemical potentials In this toy model, we vary δµ by hand. In three-flavor quark matter, the analogue of δµ is controlled by the nonzero strange quark mass and the requirement of electrical neutrality and would be of order m 2 s /4µ as in (1.2). In our toy model, we shall take the interaction between quarks to be pointlike, with the quantum numbers of single-gluon exchange. This s-wave interaction is a reasonable starting point at accessible densities but is certainly inappropriate at asymptotically high density, where the interaction between quarks (by gluon exchange) is dominated by forward scattering. The crystalline color superconducting state has been analyzed at asymptotically high densities in Refs. [18,20]. As we shall explain, the analysis of Ref. [18] indicates that the crystal structure in this circumstance would be different from that we obtain.
Our toy model may turn out to be a better model for the analysis of LOFF pairing in atomic systems. (There, the phenomenon could be called "crystalline superfluidity".) Recently, ultracold gases of fermionic atoms have been cooled down to the degenerate regime, with temperatures less than the Fermi energy [23], and reaching the pairing transition (perhaps by increasing the atom-atom interaction rather than by further reducing the temperature) seems a reasonable possibility [24].
In such systems, there really are only two species of atoms (two spin states) that pair with each other, whereas in QCD our model is a toy model for a system with nine quarks. In the atomic systems, the interaction will be s-wave dominated whereas in QCD, it remains to be seen how good this approximation is at accessible densities. Furthermore, in the atomic physics context experimentalists can control the densities of the two different atoms that pair, and in particular can tune their density difference. This means that experimentalists wishing to search for crystalline superfluidity have the ability to dial the most relevant control parameter [25,26]. In QCD, in contrast, δµ is controlled by m 2 s /µ, meaning that it is up to nature whether, and if so at what depth in a compact star, crystalline color superconductivity occurs.
In the simplest LOFF state, each Cooper pair carries momentum 2q. Although the magnitude |q| is determined energetically, the directionq is chosen spontaneously. The condensate is dominated by those regions in momentum space in which a quark pair with total momentum 2q has both members of the pair within of order ∆ of their respective Fermi surfaces. These regions form circular bands on the two Fermi surfaces, as shown in Fig. 1. Making the ansatz that all Cooper pairs make the same choice of directionq corresponds to choosing a single circular band on each Fermi surface, as in Fig. 1. In position space, it corresponds to a condensate that varies in space like This ansatz is certainly not the best choice, because it only allows a small fraction of all the quarks near their Fermi surfaces to pair. If a single plane wave is favored, why not two? That is, if one choice ofq is favored, why not add a second q, with the same |q| but a differentq, to allow more quarks near their Fermi surfaces to pair? If two are favored, why not three? Why not eight? If eight plane waves are favored, how should theirq's be oriented? These are the sort of questions we seek to answer in this paper. Upon making the plane-wave ansatz (1.4), we know from previous work [21,22,27,15] that this simplest LOFF phase is favored over the BCS state and over no pairing at all in a range δµ 1 < δµ < δµ 2 , where δµ 1 = ∆ 0 / √ 2 ≃ 0.707∆ 0 and δµ 2 ≃ 0.754∆ 0 . At δµ 2 , there is a second-order phase transition at which ∆/∆ 0 → 0 and |q|/∆ 0 tends to a nonzero limit, which we shall denote q 0 /∆ 0 , where q 0 ≃ 0.90∆ 0 ≃ 1.2δµ 2 . 1 At the second-order phase transition, the quarks that participate in the crystalline pairing lie on circular rings on their Fermi surfaces that are characterized by an opening angle ψ 0 = 2 cos −1 (δµ/|q|) ≃ 67.1 • and an angular width that is of order ∆/δµ, which therefore tends to zero. At δµ 1 there is a first-order phase transition at which the LOFF solution with gap ∆ is superseded 1 In high density QCD in the limit of a large number of colors N c , the ground state is not a color superconductor. Instead, it features a chiral condensate that varies in some crystalline pattern [28], although at least at weak coupling it seems that this phase only occurs for N c 's of order thousands [29,30]. For this chiral crystal, which is qualitatively different from a crystalline color superconductor because |q| = µ rather than |q| ∼ δµ, several possible crystal structures have been analyzed in Ref. [31]. Figure 1: In the simplest LOFF state, in which every Cooper pair carries momentum 2q, an up quark with momentum p (or p ′ ) pairs with a down quark with momentum −p + 2q (or −p ′ + 2q). In the figure q has been taken to point upwards. The pairing is dominated by pairing between up quarks in a band centered on the dashed ring shown on the up Fermi surface, and down quarks in a band centered on the solid ring shown on the down Fermi surface. The vector sum of the momenta in any pair is 2q, and the pairing bands are defined as those pairs such that both quarks in a pair are within ∆ of their respective Fermi surface. The thickness of the pairing bands is thus of order ∆. The width of the pairing bands is of order ∆μ/δµ, meaning that their angular width is of order ∆/δµ. We take a weak coupling limit in which ∆ 0 → 0 while δµ/∆ 0 and |q|/∆ 0 are held fixed and nonzero. In this limit, ∆/∆ 0 also stays fixed, and therefore so does the angular width of the pairing bands. In this limit, the two Fermi surfaces in the figure get closer and closer together and the angles ψ u and ψ d become degenerate, and take on the value ψ 0 = 2 cos −1 (δµ/|q|) ≃ 67.1 • . If, in addition, we take the further limit in which δµ/∆ 0 → δµ 2 /∆ 0 , approaching the second-order phase transition, then in this Ginzburg-Landau limit the angular extent of the pairing bands shrinks to zero, reducing them to only the rings shown.
by the BCS solution with gap ∆ 0 . (The analogue in QCD would be a LOFF window in m 2 s /µ, with CFL below and unpaired quark matter above.) Solving the full gap equation for more general crystal structures is challenging. In order to proceed, we take advantage of the second-order phase transition at δµ = δµ 2 . Because ∆ tends to zero there, we can write the gap equation as a Ginzburg-Landau expansion, working order by order in ∆. In Section 2, we develop this expansion to order ∆ 6 for a crystal made up of the sum of arbitrarily many plane waves, all with the same |q| = q 0 but with different patterns ofq's. In Section 3 we present results for a large number of crystal structures.
The Ginzburg-Landau calculation finds many crystal structures that are much more favorable than the single plane wave (1.4). For many crystal structures, it predicts a strong first-order phase transition, at some δµ * ≫ δµ 2 , between unpaired quark matter and a crystalline phase with a ∆ that is comparable in magnitude to ∆ 0 . Once δµ is reduced to δµ 2 , where the single plane wave would just be beginning to develop, these more favorable solutions already have very robust condensation energies, perhaps even larger than that of the BCS phase. These results are exciting, because they suggest that the crystalline phase is much more robust than previously thought. However, they cannot be trusted quantitatively because the Ginzburg-Landau analysis is controlled in the limit ∆ → 0, and we find a first-order phase transition to a state with ∆ = 0.
Even though it is quite a different problem, we can look for inspiration to the Ginzburg-Landau analysis of the crystallization of a solid from a liquid [32]. There too, a Ginzburg-Landau analysis predicts a first-order phase transition, and thus predicts its own quantitative downfall. But, qualitatively it is correct: it predicts the formation of a body-centered-cubic crystal and experiment shows that most elementary solids are body-centered cubic near their first-order crystallization transition.
Thus inspired, let us ask what crystal structure our Ginzburg-Landau analysis predicts for the crystalline color superconducting phase. To order ∆ 2 , we learn that |q| ≃ 1.2δµ even when δµ = δµ 2 . We also learn that, apparently, the more plane waves the better. We learn nothing about the preferred arrangement of the set of q's. By extending the calculation to order ∆ 4 and ∆ 6 , we find: • Crystal structures with intersecting pairing rings are strongly disfavored.
Recall that eachq is associated with pairing among quarks that lie on one ring of opening angle ψ 0 ≃ 67.1 • on each Fermi surface. We find that any crystal structure in which such rings intersect pays a large free energy price. The favored crystal structures should therefore be those whose set ofq's are characterized by the maximal number of nonintersecting rings. The maximal number of nonintersecting rings with opening angle 67.1 • that fit on a sphere is nine [33,34].
• We also find that those crystal structures characterized by a set of q's within which there are a large number of combinations satisfying q 1 −q 2 +q 3 −q 4 = 0 and a large number of combinations satisfying q 1 − q 2 + q 3 − q 4 + q 5 − q 6 = 0 are favored. Speaking loosely, "regular" structures are favored over "irregular" structures. The only configurations of nine nonintersecting rings are rather irregular, whereas if we limit ourselves to eight rings, there is a regular choice which is favored by this criterion: choose eight q's pointing towards the corners of a cube. In fact, a deformed cube which is slightly taller or shorter than it is wide (a cuboid) is just as good.
These qualitative arguments are supported by the quantitative results of our Ginzburg-Landau analysis, which does indeed indicate that the most favored crystal structure is a cuboid that is very close to a cube. This crystal structure is so favorable that the coefficient of ∆ 6 in the Ginzburg-Landau expression for the free energy, which we call γ, is large and negative. (In fact, we find several crystal structures with negative γ, but the cube has by far the most negative γ.) We could go on, to ∆ 8 or higher, until we found a Ginzburg-Landau free energy for the cube which is bounded from below. However, we know that this free energy would give a strongly first-order phase transition, meaning that the Ginzburg-Landau analysis would anyway not be under quantitative control. A better strategy, then, is to use the Ginzburg-Landau analysis to understand the physics at a qualitative level, as we have done. With an understanding of why a crystal structure with eight plane waves whose wave vectors point to the corners of a cube is strongly favored in hand, the next step would be to make this ansatz and solve the gap equation without making a Ginzburg-Landau approximation. We leave this to future work. The eight q's that describe the crystal structure we have found to be most favorable are eight of the vectors in the reciprocal lattice of a face-centered-cubic crystal. Thus, as we describe and depict in Section 4, the crystal structure we propose is one in which the condensate forms a face-centered-cubic lattice in position space.

The gap equation
We study the crystalline superconducting phase in a toy model for QCD that has two massless flavors of quarks and a pointlike interaction. The Lagrange function is The τ 's are Pauli matrices in flavor space, so the up and down quarks have chemical potentials as in (1.2). The vertex is Γ A = γ µ T a so that our pointlike interaction mimics the spin, color, and flavor structure of one-gluon exchange. (The T a are color SU(3) generators normalized so that tr(T a T b ) = 2δ ab .) We denote the coupling constant in the model by λ. It is convenient to use a Nambu-Gorkov diagrammatic method to obtain the gap equation for the crystalline phase. Since we are investigating a phase with spatial inhomogeneity, we begin in position space. We introduce the two-component spinor Ψ(x) = (ψ(x),ψ T (x)) and the quark propagator iS(x, x ′ ) = Ψ(x)Ψ(x ′ ) , which has PSfrag replacements ∆ = 0000 0000 1111 1111 Figure 2: The Schwinger-Dyson graph for the LOFF gap parameter ∆. The black dot is the pointlike interaction vertex and the double line represents the full anomalous propagator F , which is given in terms of ∆ in Fig. 4.
"normal" and "anomalous" components G and F , respectively: . (2. 2) The conjugate propagatorsF andḠ satisfy The gap parameter ∆(x) that describes the diquark condensate is related to the anomalous propagator F by a Schwinger-Dyson equation illustrated diagrammatically in Fig. 2. In our toy model, we are neglecting quark masses and thus the normal part of the one-particle-irreducible self-energy is zero; the anomalous part of the 1PI self energy is just ∆(x). The crystal order parameter ∆(x) defined by (2.5) is a matrix in spin, flavor and color space. In the meanfield approximation, we can use the equations of motion for Ψ(x) to obtain a set of coupled equations that determine the propagator functions in the presence of the diquark condensate characterized by ∆(x): where∆(x) = γ 0 ∆(x) † γ 0 . Any function ∆(x) that solves equations (2.5) and (2.6) is a stationary point of the free energy functional Ω[∆(x)]; of these stationary points, the one with the lowest Ω describes the ground state of the system. Our task, then, is to invert (2.6), obtaining F in terms of ∆(x), substitute in (2.5), find solutions for ∆(x), and then evaluate Ω for all solutions we find. There are some instances where analytic solutions to equations (2.5) and (2.6) can be found. The simplest case is that of a spatially uniform condensate. Translational invariance then implies that the propagators are diagonal in momentum space: S(p, p ′ ) = S(p)(2π) 4 δ (4) (p − p ′ ). In this case, Eqs. (2.6) immediately yield which is easily inverted to obtain S, which can then be substituted on the right-hand side of Eq. (2.5) to obtain a self-consistency equation (i.e. a gap equation) for ∆. The solution of this gap equation describes the familiar "2SC" phase [1,2,3,4,5], a two-flavor, two-color BCS condensate, where T 2 , τ 2 , and Cγ 5 indicate that the condensate is a color antitriplet, flavor singlet, and Lorentz scalar, respectively. 2 The remaining factor ∆ 0 , which without loss of generality can be taken to be real, gives the magnitude of the condensate. In order to solve the resulting gap equation for ∆ 0 , we must complete the specification of our toy model by introducing a cutoff. In previous work [3,4,6,5], it has been shown that if, for a given cutoff, the coupling λ is chosen so that the model describes a reasonable vacuum chiral condensate, then at µ ∼ 400 − 500 MeV the model describes a diquark condensate which has ∆ 0 of order tens to 100 MeV. Ratios between physical observables depend only weakly on the cutoff, meaning that when λ is taken to vary with the cutoff such that one observable is held fixed, others depend only weakly on the cutoff. For this reason, we are free to make a convenient choice of cutoff so long as we then choose the value of λ that yields the "correct" ∆ 0 . Since we do not really know the correct value of ∆ 0 and since this is after all only a toy model, we simply think of ∆ 0 as the single free parameter in the model, specifying the strength of the interaction and thus the size of the BCS condensate. Because the quarks near the Fermi surface contribute most to pairing, it is convenient to introduce a cutoff ω defined so as to restrict the gap integral to momentum modes near the Fermi surface (||p| −μ| ≤ ω). In the weak coupling (small λ) limit, the explicit solution to the gap equation is then This is just the familiar BCS result for the gap. (Observe that the density of states at the Fermi surface is N 0 = 2μ 2 /π 2 .) We denote the gap for this BCS solution by ∆ 0 , reserving the symbol ∆ for the gap parameter in the crystalline phase. We shall see explicitly below that when we express our results for ∆ relative to ∆ 0 , they are completely independent the cutoff ω as long as ∆ 0 /µ is small. The BCS phase, with ∆ 0 given by (2.9), has a lower free energy than unpaired quark matter as long as δµ < δµ 1 = ∆ 0 / √ 2 [35]. The first-order unpairing transition at δµ = δµ 1 is the analogue in our two-flavor toy model of the unlocking transition in QCD. For δµ > δµ 1 , the free energy of any crystalline solution we find below must be compared to that of unpaired quark matter; for δµ < δµ 1 , crystalline solutions should be compared to the BCS phase. We shall work at δµ > δµ 1 .
The simplest example of a LOFF condensate is one that varies like a plane wave: . The condensate is static, meaning that q = (0, q). We shall denote |q| by q 0 . In this condensate, the momenta of two quarks in a Cooper pair is (p + q, −p + q) for some p, meaning that the total momentum of each and every pair is 2q. See Refs. [16,19] for an analysis of this condensate using the Nambu-Gorkov formalism. Here, we sketch the results. If we shift the definition of Ψ in momentum space to Ψ q (p) ≡ (ψ(p + q),ψ(−p + q)), then in this shifted basis the propagator is diagonal: and the inverse propagator is simply See Refs. [16,19] for details and to see how this equation can be inverted and substituted into Eq. (2.5) to obtain a gap equation for ∆. This gap equation has nonzero solutions for δµ < δµ 2 ≃ 0.7544∆ 0 , and has a second-order phase transition at δµ = δµ 2 with ∆ ∼ (δµ 2 − δµ) 1/2 . We rederive these results below.
If the system is unstable to the formation of a single plane-wave condensate, we might expect that a condensate of multiple plane waves is still more favorable. Again our goal is to find gap parameters ∆(x) that are self-consistent solutions of Eqs. (2.5) and (2.6). We use an ansatz that retains the Lorentz, flavor, and color structure of the 2SC phase: but now ∆(x) is a scalar function that characterizes the spatial structure of the crystal. We write this function as a superposition of plane waves: where, as before, q = (0, q). The {∆ q } constitute a set of order parameters for the crystalline phase. Our task is to determine for which set of q's the ∆ q 's are nonzero. Physically, for each ∆ q = 0 the condensate includes some Cooper pairs for which the total momentum of a pair is 2q. This is indicated by the structure of the anomalous propagator F in momentum space: Eqs. (2.5) and (2.13) together imply that Figure 3: A process whereby a quark with momentum p scatters by interactions with two plane-wave condensates and acquires a momentum p − 2q 1 + 2q 2 . and Eq. (2.15) yields an infinite set of coupled gap equations, one for each q. (Note that each F q depends on all the ∆ q 's.) It is not consistent to choose only a finite set of ∆ q to be nonzero because when multiple plane-wave condensates are present, these condensates induce an infinite "tower" (or lattice) of higher momentum condensates. This is easily understood by noting that a quark with momentum p can acquire an additional momentum 2q 2 − 2q 1 by interacting with two different plane-wave condensates as it propagates through the medium, as shown in Fig. 3. Note that this process cannot occur when there is only a single plane-wave condensate. The analysis of the single plane-wave condensate closes with only a single nonzero ∆ q , and is therefore much easier than the analysis of a generic crystal structure. Another way that this difficulty manifests itself is that once we move beyond the single plane-wave solution to a more generically nonuniform condensate, it is no longer possible to diagonalize the propagator in momentum space by a shift, as was possible in Eqs. (2.10, 2.11).

The Ginzburg-Landau approximation
The infinite system of equations (2.15) has been solved analytically only in one dimension, where it turns out that the gap parameter can be expressed as a Jacobi elliptic function that, as promised, is composed of an infinite number of plane waves [36]. In three dimensions, the crystal structure of the LOFF state remains unresolved [37,38,26]. In the vicinity of the second-order transition at δµ 2 , however, we can simplify the calculation considerably by utilizing the smallness of ∆ to make a controlled Ginzburg-Landau approximation. This has the advantage of providing a controlled truncation of the infinite series of plane waves, because near δµ 2 the system is unstable to the formation of plane-wave condensates only for q's that fall on a sphere of a certain radius q 0 , as we shall see below. This was in fact the technique employed by Larkin and Ovchinnikov in their original paper [21], and it has been further developed in Refs. [37,38,26]. As far as we know, though, no previous authors have done as complete a study of possible crystal structures in three dimensions as we attempt. Most have limited their attention to, at most, structures 1, 2, 5 and 9 from the 23 structures we describe in Fig. 10 and Table 1 below. As Sfrag replacements far as we know, no previous authors have investigated the crystal structure that we find to be most favorable.
The authors of Refs. [37,38,26] have focused on using the Ginzburg-Landau approximation at nonzero temperature, near the critical temperature at which the LOFF condensate vanishes. Motivated by our interest in compact stars, we follow Larkin and Ovchinnikov in staying at T = 0 while using the fact that, for a single plane-wave condensate, ∆ → 0 for δµ → δµ 2 to motivate the Ginzburg-Landau approximation. The down side of this is that, in agreement with previous authors, we find that the T = 0 phase transition becomes first order when we generalize beyond a single plane wave. In the end, therefore, the lessons of our Ginzburg-Landau approximation must be taken qualitatively. We nevertheless learn much that is of value.
To proceed with the Ginzburg-Landau expansion, we first integrate equations (2.6) to obtain Then we expand these equations in powers of the gap function ∆(x). ForF (x, x ′ ) we find (suppressing the various spatial coordinates and integrals for notational simplicity) as expressed diagrammatically in Fig. 4. We then substitute this expression forF into the right-hand side of the Schwinger-Dyson equation (actually the conjugate of equation (2.15)). After some spin, color, and flavor matrix manipulation, the result in momentum space is as shown in Fig. 5. The prefactors have been chosen for later convenience. The functions Π, J, and K corresponding to the three graphs in Fig. 5 are given by: (2.20) We shall see that δµ and |q| are both of order ∆ which in turn is of order ∆ 0 . This means that all these quantities are much less thanμ in the weak coupling limit. Thus, in the weak coupling limit we can choose the cutoff ω such that δµ, |q| ≪ ω ≪μ. In this limit, J and K are independent of the cutoff ω, as we shall see in the appendix where we present their explicit evaluation. In this limit, where we have used the explicit solution to the BCS gap equation (2.9) to eliminate the cutoff ω in favor of the BCS gap ∆ 0 , and where the last equation serves to define α. Note that α depends on the cutoff ω only through ∆ 0 , and depends only on the ratios |q|/∆ 0 and δµ/∆ 0 . It will prove convenient to use the definition of α to rewrite the Ginzburg-Landau equation Eq. (2.19) as To learn how to interpret α, consider the single plane-wave condensate in which ∆ q = 0 only for a single q. If we divide equation (2.19) by ∆ * q , we see that the equation Π = −π 2 /2λμ 2 , which is to say α = 0, defines a curve in the space of (|q|, δµ) where we can find a solution to the gap equation with ∆ q → 0, with |q| on the curve and for anyq. This curve is shown in Fig. 6. We shall see below that when only one ∆ q is nonzero, the "J sum" and "K sum" in (2.19) are both positive. This means that wherever α < 0, i.e. below the solid curve in Fig. 6, there are solutions with ∆ q = 0 for these values of |q|, and wherever α > 0, i.e. above the solid curve, there are no single plane-wave solutions to the gap equation. The solid curve in Fig. 6 therefore marks the boundary of the instability towards the formation of a single plane-wave condensate. The highest point on this curve is special, as it denotes the maximum value of δµ for which a single plane-wave LOFF condensate can arise. This second-order critical point occurs at (|q|, δµ) = (q 0 , δµ 2 ) with δµ 2 ≃ 0.7544∆ 0 and q 0 /δµ 2 ≃ 1.1997, where ∆ 0 is the BCS gap of Eq. (2.9).
As δµ → δµ 2 from above, only those plane waves lying on a sphere in momentum space with |q| = q 0 are becoming unstable to condensation. If we analyze them one by one, all these plane waves are equally unstable. That is, in the vicinity of the critical point δµ 2 , the LOFF gap equation admits plane-wave condensates with ∆ q = 0 for a single q lying somewhere on the sphere |q| = q 0 . For each such plane wave, the paired quarks occupy a ring with opening angle ψ 0 = 2 cos −1 (δµ/|q|) ≃ 67.1 • on each Fermi surface, as shown in Fig. 1.

The free energy
In order to compare different crystal structures, with (2.22) in hand, we can now derive a Ginzburg-Landau free energy functional Ω[∆(x)] which characterizes the system in the vicinity of δµ 2 , where ∆ → 0. This is most readily obtained by noting that the gap equations (2.22) must be equivalent to because solutions to the gap equations are stationary points of the free energy. This determines the free energy up to an overall multiplicative constant, which can be found by comparison with the single plane-wave solution previously known. The result is where N 0 = 2μ 2 /π 2 and where we have restricted our attention to modes with |q| = q 0 , as we now explain. Note that in the vicinity of |q| = q 0 and δµ = δµ 2 , We see that for α > 0 (that is, for δµ > δµ 2 ) ∆ q = 0 is stable whereas for α < 0 (that is, δµ < δµ 2 ), the LOFF instability sets in. In the limit α → 0 − , only those plane waves on the sphere |q| = q 0 are unstable. For this reason we only include these plane waves in the expression (2.24) for the free energy. We shall do most of our analysis in the vicinity of δµ = δµ 2 , where we choose |q| = q 0 = 1.1997δµ 2 as just described. However, we shall also want to apply our results at δµ > δµ 2 . At these values of δµ, we shall choose |q| in such a way as to minimize α(|q|), because this minimizes the quadratic term in the free energy Ω and thus minimizes the free energy in the vicinity of ∆ → 0, which is where the Ginzburg-Landau analysis is reliable. As Fig. 6 indicates, for any given δµ the minimum value of α is to be found at |q| = 1.1997δµ. Therefore, when we apply Eq. (2.24) away from δµ 2 , we shall set q 0 = 1.1997δµ, just as at δµ = δµ 2 . As a consequence, the opening angle of the pairing rings, ψ 0 = 2 cos −1 (δµ/|q|), is unchanged when we move away from δµ = δµ 2 .

Generalities
All of the modes on the sphere |q| = q 0 become unstable at δµ = δµ 2 . The quadratic term in the free energy includes no interaction between modes with different q's, and so predicts that ∆ q = 0 for all modes on the sphere. Each plane-wave mode corresponds to a ring of paired quarks on each Fermi surface, so we would obtain a cacophony of multiple overlapping rings, favored by the quadratic term because this allows more and more of the quarks near their respective Fermi surfaces to pair. Moving beyond lowest order, our task is to evaluate the quartic and sextic terms in the free energy. These higher order terms characterize the effects of interactions between ∆ q 's with differing q's (between the different pairing rings) and thus determine how condensation in one mode enhances or deters condensation in other modes. The results we shall present rely on our ability to evaluate J and K, defined in Eqs. (2.20). We describe the methods we use to evaluate these expressions in an appendix and focus here on describing and understanding the results. We shall see, for example, that although the quadratic term favors adding more rings, the higher order terms strongly disfavor configurations in which ∆ q 's corresponding to rings that intersect are nonzero. Evaluating the quartic and sextic terms in the free energy will enable us to evaluate the free energy of condensates with various configurations of several plane waves and thereby discriminate between candidate crystal structures.  A given crystal structure can be described by a set of vectors Q = {q a , q b , · · ·}, specifying which plane wave modes are present in the condensate, and a set of gap parameters {∆ qa , ∆ q b , · · ·}, indicating the amplitude of condensation in each of the modes. Let us define G as the group of proper and improper rotations that preserve the set Q. We make the assumption that G is also the point group of the crystal itself; this implies that ∆ q = ∆ q ′ if q ′ is in the orbit of q under the group action. For most (but not all) of the structures we investigate, Q has only one orbit and therefore all of the ∆ q 's are equal.
For a given set Q, the quartic term in the free energy (2.24) is a sum over all combinations of four q's that form closed "rhombuses", as shown in Fig. 7. The four q's are chosen from the set Q and they need not be distinct. By a rhombus we mean a closed figure composed of four equal length vectors which will in general be nonplanar. A rhombus is therefore characterized by two internal angles (ψ, χ) with the constraint 0 ≤ ψ + χ ≤ π. Each shape corresponds to a value of the J function (as defined in equations (2.20)); the rotational invariance of the J function implies that congruent shapes give the same value and therefore J(q 1 q 2 q 3 q 4 ) = J(ψ, χ). So, each unique rhombic combination of q's in the set Q that characterizes a given crystal structure yields a unique contribution to the quartic coefficient in the Ginzburg-Landau free energy of that crystal structure. The continuation to next order is straightforward: the sextic term in the free energy (2.24) is a sum over all combinations of six q's that form closed "hexagons", as shown in Fig. 7. Again these shapes are generally nonplanar and each unique hexagonal combination of q's yields a unique value of the K function and a unique contribution to the sextic coefficient in the Ginzburg-Landau free energy of the crystal.
When all of the ∆ q 's are equal, we can evaluate aggregate quartic and sextic coefficients β and γ, respectively, as sums over all rhombic and hexagonal combinations of the q's in the set Q: Then, for a crystal with P plane waves, the free energy has the simple form and we can analyze a candidate crystal structure by calculating the coefficients β and γ and studying the resultant form of the free energy function.
If β and γ are both positive, a second-order phase transition occurs at α = 0; near the critical point the value of γ is irrelevant and the minimum energy solution is ∆ = P |α| β 1 2 , for δµ ≤ δµ 2 (i.e. α ≤ 0). If β is negative and γ is positive, the phase transition is in fact first order and occurs at a new critical point defined by In order to find the δµ * corresponding to α * , we need to solve Since α * is positive, the critical point δµ * at which the first-order phase transition occurs is larger than δµ 2 . If α * is small, then δµ * ≃ (1 + α * )δµ 2 . Thus, a crystalline color superconducting state whose crystal structure yields a negative β and positive γ persists as a possible ground state even above δµ 2 , the maximum δµ at which the plane-wave state is possible. At the first-order critical point (3.4), the free energy has degenerate minima at If we reduce δµ below δµ * , the minimum with ∆ = 0 deepens. Once δµ is reduced to the point at which the single plane wave would just be starting to form with a free energy infinitesimally below zero, the free energy of the crystal structure with negative β and positive γ has ∆ = |β| γ 1 2 , PSfrag replacements ψ ψ ψ Figure 8: Two different "hexagonal" shapes (as in Fig. 7) that can be constructed from two vectors q a and q b . These shapes correspond to the functions K 1 (ψ) and K 2 (ψ) in Eq. (3.11).
Finally, if γ is negative, the order ∆ 6 Ginzburg-Landau free energy is unbounded from below. In this circumstance, we know that we have found a first order phase transition but we do not know at what δµ * it occurs, because the stabilization of the Ginzburg-Landau free energy at large ∆ must come about at order ∆ 8 or higher.

One wave
With these general considerations in mind we now proceed to look at specific examples of crystal structures. We begin with the single plane-wave condensate (P = 1). The quartic coefficient of the free energy is 8) and the sextic coefficient is yielding a second-order phase transition at α = 0. These coefficients agree with those obtained by expanding the all-orders-in-∆ solution for the single plane wave which can be obtained by variational methods [22,27,15] or by starting from (2.11), as in Refs. [16,19]. The coefficient β in (3.8) was first found by Larkin and Ovchinnikov [21].

Two waves
Our next example is a condensate of two plane waves (P = 2) with wave vectors q a and q b and equal gaps ∆ qa = ∆ q b = ∆. The most symmetrical arrangement is an antipodal pair (q b = −q a ), which yields a cosine spatial variation ∆(x) ∼ cos(2q a · x). We will find it useful, however, to study the generic case where q a  : β(ψ) and γ(ψ), the quartic and sextic coefficients in the Ginzburg-Landau free energy for a condensate consisting of two plane waves whose wave vectors define an angle ψ. and q b have the same magnitude but define an arbitrary angle ψ. We find that the quartic coefficient is β(ψ) = 2J 0 + 4J(ψ, 0) (3.10) and the sextic coefficient is where K 1 (ψ) = K(q a q a q a q a q b q b ) and K 2 (ψ) = K(q a q a q b q a q a q b ). (K 1 and K 2 arise from the "hexagonal" shapes shown in Fig. 8.) The functions β(ψ) and γ(ψ) are plotted in Fig. 9. These functions manifest a number of interesting features. Notice that the functions are singular and discontinuous at a critical angle ψ = ψ 0 ≃ 67.1 • , where ψ 0 is the opening angle of a LOFF pairing ring on the Fermi surface. For the two-wave condensate we have two such rings, and the two rings are mutually tangent when ψ = ψ 0 . For ψ < ψ 0 , both β and γ are large and positive, implying that an intersecting ring configuration is energetically unfavorable. For ψ > ψ 0 the functions are relatively flat and small, indicating some indifference towards any particular arrangement of the nonintersecting rings. There is a range of angles for which β is negative and a first-order transition occurs (note that γ is always positive). The favored arrangement is a pair of adjacent rings that nearly intersect (ψ = ψ 0 + ǫ).
It is unusual to find coefficients in a Ginzburg-Landau free energy that behave discontinuously as a function of parameters describing the state, as seen in Fig. 9. These discontinuities arise because, as we described in the caption of Fig. 1, we are taking two limits. We first take a weak coupling limit in which ∆ 0 , ∆, δµ, |q| ≪ ω ≪ µ while δµ/∆ 0 , |q|/∆ 0 and ∆/δµ (and thus the angular width of the pairing bands) are held fixed. Then, we take the Ginzburg-Landau limit in which δµ/∆ 0 → δµ 2 /∆ 0 and ∆/∆ 0 → 0 and the pairing bands shrink to rings of zero angular width. In the Ginzburg-Landau limit, there is a sharp distinction between ψ < 67.1 • where the rings intersect and ψ > 67.1 • where they do not. Without taking the weak coupling limit, the plots of β(ψ) and γ(ψ) would nevertheless look like smoothed versions of those in Fig. 9, smoothed on angular scales of order δµ 2 /ω 2 (for β) and δµ 4 /ω 4 (for γ). However, in the weak coupling limit these small angular scales are taken to zero. Thus, the double limit sharpens what would otherwise be distinctive but continuous features of the coefficients in the Ginzburg-Landau free energy into discontinuities.

Crystals
From our analysis of the two-wave condensate, we can infer that for a general multiple-wave condensate it is unfavorable to allow the pairing rings to intersect on the Fermi surface. For nonintersecting rings, the free energy should be relatively insensitive to how the rings are arranged on the Fermi surface. However, Eqs. (3.1) suggest that a combinatorial advantage is obtained for exceptional structures that permit a large number of rhombic and hexagonal combinations of wave vectors. That is, if there are many ways of picking four (not necessarily different) wave vectors from the set of wave vectors that specify the crystal structure for which q 1 − q 2 + q 3 − q 4 = 0, or if there are many ways of picking six wave vectors for which q 1 −q 2 + q 3 −q 4 + q 5 −q 6 = 0, such a crystal structure enjoys a combinatorial advantage that will tend to make the magnitudes of β or γ large. For a rhombic combination q 1 − q 2 + q 3 − q 4 = 0, the four q's must be the four vertices of a rectangle that is inscribed in a circle on the sphere |q| = q 0 . (The circle need not be a great circle, and the rectangle can degenerate to a line or a point if the four q's are not distinct). For a hexagonal combination q 1 − q 2 + q 3 − q 4 + q 5 − q 6 = 0, the triplets (q 1 q 3 q 5 ) and (q 2 q 4 q 6 ) are vertices of two inscribed triangles that have a common centroid. In the degenerate case where only four of the six q's are distinct, the four distinct q's must be the vertices of an inscribed rectangle or an inscribed isoceles trapezoid for which one parallel edge is twice the length of the other. When five of the six q's are distinct, they can be arranged as a rectangle plus any fifth point, or as five vertices of an inscribed cuboid arranged as one antipodal pair plus the three corners adjacent to one of the antipodes.
We have investigated a large number of different multiple-wave configurations depicted in Fig. 10 and the results are compiled in Table 1. The name of each configuration is the name of a polygon or polyhedron that is inscribed in a sphere of radius q 0 ; the P vertices of the given polygon or polyhedron then correspond to the P wave vectors in the set Q. With this choice of nomenclature, keep in mind that what we call the "cube" has a different meaning than in much of the previous literature. We refer to an eight plane-wave configuration with the eight wave vectors directed at the eight corners of a cube. Because this is equivalent to eight vectors directed at the eight faces of an octahedron -the cube and the octahedron are dual polyhedra -in the nomenclature of previous literature this eight-wave crystal would have been called an octahedron, rather than a cube. Similarly, the crystal that we call the "octahedron" (six plane waves whose wave vectors point at the six corners of an octahedron) is the structure that has been called a cube in the previous literature, because its wave vectors point at the faces of a cube.
Because the LOFF pairing rings have an opening angle ψ 0 ≃ 67.1 • , no more than nine rings can be arranged on the Fermi surface without any intersection [33,34]. For this reason we have focussed on crystal structures with nine or fewer waves, but we have included several structures with more waves in order to verify that such structures are not favored. We have tried to analyze a fairly exhaustive list of candidate structures. All five Platonic solids are included in Table 1, as is the simplest Archimedean solid, the cuboctahedron. (All other Archimedean solids have even more vertices.) We have analyzed many dihedral polyhedra and polygons: regular polygons, bipyramids, prisms, 3 antiprisms, 4 and various capped or augmented polyhedra 5 . For each crystal structure we list the crystal point group G and the Föppl configuration of the polyhedron or polygon. The Föppl configuration is a list of the number of vertices on circles formed by intersections of the sphere with consecutive planes perpendicular to the principal symmetry axis of the polyhedron or polygon. We use a modified notation where a orā indicates that the points on a given circle are respectively eclipsed or staggered relative to the circle above. Note that polyhedra with several different principal symmetry axes, namely those with T , O, or I symmetry, have several different Föppl descriptions: for example, a cube is (44) along a fourfold symmetry axis or (1331) along a threefold symmetry axis. (That is, the cube can equally be described as a square prism or a bicapped trigonal 3 A trigonal prism is two triangles, one above the other. A cube is an example of a square prism. 4 An antiprism is a prism with a twist. For example, a square antiprism is two squares, one above the other, rotated relative to each other by 45 • . The octahedron is an example of a trigonal antiprism. 5 Capping a polyhedron adds a single vertex on the principal symmetry axis of the polyhedron (or polygon). Thus, a capped square is a square pyramid and an octahedron could be called a bicapped square. Augmenting a polyhedron means adding vertices on the equatorial plane, centered outside each vertical facet. Thus, augmenting a trigonal prism adds three new vertices.  Candidate crystal structures with P plane waves, specified by their symmetry group G and Föppl configuration. Bars denote dimensionless equivalents: β = β δµ 2 ,γ = γ δµ 4 ,Ω = Ω/(δµ 2 2 N 0 ) with N 0 = 2μ 2 /π 2 .Ω min is the (dimensionless) minimum free energy at δµ = δµ 2 , obtained from (3.7). The phase transition (first order forβ < 0 andγ > 0, second order forβ > 0 andγ > 0) occurs at δµ * . antiprism. This should make clear that the singly capped trigonal antiprism of Fig. 10 and Table 1 is a cube with one vertex removed.) We do not claim to have analyzed all possible crystal structures, since that is an infinite task. However, there are several classic mathematical problems regarding extremal arrangements of points on a sphere and, although we do not know that our problem is related to one of these, we have made sure to include solutions to these problems. For example, many of the structures that we have evaluated correspond to solutions of Thomson's problem [39,33] (lowest energy arrangement of P point charges on the surface of a sphere) or Tammes's problem [34,33] (best packing of P equal circles on the surface of a sphere without any overlap). In fact, we include all solutions to the Thomson and Tammes problems for P ≤ 9. Our list also includes all "balanced" configurations [40] that are possible for nonintersecting rings: a balanced configuration is a set Q with a rotational symmetry about every q ∈ Q; this corresponds to an arrangement of particles on a sphere for which the particles are in equilibrium for any two-particle force law.
For each crystal structure, we have calculated the quartic and sextic coefficients β and γ according to Eqs. (3.1), using methods described in the appendix to calculate all the J and K integrals. To further discriminate among the various candidate structures, we also list the minimum free energy Ω min evaluated at the plane-wave instability point δµ = δµ 2 where α = 0. To set the scale, note that the BCS state at δµ = 0 has Ω BCS = −µ 2 ∆ 2 0 /π 2 , corresponding toΩ BCS = −0.879 in the units of Table 1. For those configurations with β > 0 and γ > 0, Ω min = 0 at δµ = δµ 2 and Ω min < 0 for δµ < δµ 2 , where α < 0. Thus, we find a second-order phase transition at δµ = δµ 2 . For those configurations with β < 0 and γ > 0, at δµ = δµ 2 the minimum free energy occurs at a nonzero ∆ with Ω min < 0. (The value of ∆ at which this minimum occurs can be obtained from (3.7).) Because Ω min < 0 at δµ = δµ 2 , if we go to δµ > δµ 2 , where α > 0, we lift this minimum until at some δµ * it has Ω = 0 and becomes degenerate with the ∆ = 0 minimum. At δµ = δµ * , a first-order phase transition occurs. For a very weak first-order phase transition, δµ * ≃ δµ 2 ≃ 0.754∆ 0 . For a strong first-order phase transition, δµ * ≫ δµ 2 and the crystalline color superconducting phase prevails as the favored ground state over a wider range of δµ.

Crystal structures with intersecting rings lose
There are seven configurations in Table 1 with very large positive values for γ. These are precisely the seven configurations that have intersecting pairing rings: the hexagon, hexagonal prism, capped square prism, bicapped square antiprism, icosahedron, cuboctahedron, and dodecahedron. The first two of these include hexagons, and since ψ 0 > 60 • the rings intersect. The last four of these crystal structures have more than nine rings, meaning that intersections between rings are also inevitable. The capped square prism is an example of a nine-wave structure with intersecting rings. It has a γ which is almost two orders of magnitude larger than that of the augmented trigonal prism and the capped square antiprism which, in contrast, are nine wave structures with no intersecting rings. Because of their very large γ's all the structures with intersecting rings have either second-order phase transitions or very weak first-order phase transitions occurring at a δµ * ≃ δµ 2 . At δµ = δµ 2 , all these crystal structures have Ω min very close to zero. Thus, as our analysis of two plane waves led us to expect, we conclude that these crowded configurations with intersecting rings are disfavored.

"Regular" crystal structures rule, and the cube rules them all
At the opposite extreme, we see that there are several structures that have negative values of γ: the square, square pyramid, square antiprism, trigonal prism, capped trigonal antiprism, and cube. Our analysis demonstrates that the transition to all these crystal structures (as to those with β < 0 and γ > 0) is first order. But, we cannot evaluate Ω min or δµ * because, to the order we are working, Ω is unbounded from below. For each of these crystal structures, we could formulate a well-posed (but difficult) variational problem in which we make a variational ansatz corresponding to the structure, vary, and find Ω min without making a Ginzburg-Landau approximation. It is likely, therefore, that within the Ginzburg-Landau approximation Ω will be stabilized at a higher order than the sextic order to which we have worked.
Of the sixteen crystal structures with no intersecting rings, there are seven that are particularly favored by the combinatorics of Eqs. (3.1). It turns out that these seven crystal structures are precisely the six that we have found with γ < 0, plus the octahedron, which is the most favored crystal structure among those with γ > 0. As discussed earlier, more terms contribute to the rhombic and hexagonal sums in Eqs. (3.1) when the q's are arranged in such away that their vertices form rectangles, trapezoids, and cuboids inscribed in the sphere |q| = q 0 . Thus the square itself fares well, as do the square pyramid, square antiprism, and trigonal prism which contain one, two, and three rectangular faces, respectively. The octahedron has three square cross sections. However, the cube is the outstanding winner because it has six rectangular faces, six rectangular cross sections, and also allows the five-corner arrangements described previously. The capped trigonal antiprism in Table 1 is a cube with one vertex removed. This seven-wave crystal has almost as many waves as the cube, and almost as many combinatorial advantages as the cube, and it turns out to have the second most negative γ.
With eight rings, the cube is close to the maximum packing for nonintersecting rings with opening angle ψ 0 ≃ 67.1 • . Although nine rings of this size can be packed on the sphere, adding a ninth ring to the cube and deforming the eight rings into a cuboid, as we have done with the capped square prism, necessarily results in intersecting rings and the ensuing cost overwhelms the benefits of the cuboidal structure. To form a nine-ring structure with no intersections requires rearranging the eight rings, spoiling the favorable regularities of the cuboid. Therefore a nine-ring arrangement is actually less favorable than the cuboid, even though it allows one more plane wave. We see from Table 1 that γ for the cube is much more negative than that for any of the other combinatorially favored structures. The cube is our winner, and we understand why.
To explore the extent to which the cube is favored, we can compare it to the octahedron, which is the crystal structure with γ > 0 for which we found the strongest first-order phase transition, with the largest δµ * and the deepest Ω min . The order ∆ 6 free energy we have calculated for the cube is far below that for the octahedron at all values of ∆. To take an extreme example, at δµ = δµ 2 the octahedron hasΩ min = −13.365 at ∆ = 1.263δµ 2 = 0.953∆ 0 whereas for the cube we find thatΩ = −2151.5 at this ∆. As another example, suppose that we arbitrarily add + 1 4 800∆ 8 /δµ 8 to theΩ of the cube. In this case, at δµ = δµ 2 we find that the cube hasΩ min = −32.5 at ∆ = 0.656∆ 0 and is thus still favored over the octahedron, even though we have not added any ∆ 8 term to the free energy of the octahedron. These numerical exercises demonstrate the extraordinary robustness of the cube, but should not be taken as more than qualitative. We do not know at what ∆ and at what value Ω min the true free energy for the cube finds its minimum. However, because the qualitative features of the cube are so favorable we expect that it will have a deeper Ω min and a larger δµ * than the octahedron. Within the Ginzburg-Landau approximation, the octahedron already has ∆ = 0.953∆ 0 and a deepΩ min = −13.365, about fifteen times deeper thanΩ BCS = −0.879 for the BCS state at δµ = 0.
Even if we were to push the Ginzburg-Landau analysis of the cube to higher order and find a stable Ω min , we would not be able to trust such a result quantitatively. Because it predicts a strong first-order phase transition, the Ginzburg-Landau approximation predicts its own quantitative demise. What we have learned from it, however, is that there are qualitative reasons that make the cube the most favored crystal structure of them all. And, to the extent that we can trust the quantitative calculations qualitatively, they indicate that the first-order phase transition results in a state with ∆ comparable to or bigger than that in the BCS phase, with Ω min comparable to or deeper than that of the BCS phase, and occurs at a δµ * ≫ δµ 2 .

Varying continuous degrees of freedom
None of the regularities of the cube which make it so favorable are lost if it is deformed continuously into a cuboid, slightly shorter or taller than it is wide, as long as it is not deformed so much as to cause rings to cross. Next, we investigate this and some of the other possible continuous degrees of freedom present in a number of the crystal structures we have described above.
So far we have neglected the fact, mentioned at the start of this section, that some of the candidate structures have multiple orbits under the action of the point group G. These structures include the square pyramid, the four bipyramids, and the five capped or augmented structures listed in Table 1; all have two orbits except for the three singly capped crystal structures, which have three orbits. For these multiple-orbit structures each orbit should have a different gap parameter but in Table 1 we have assumed that all the gaps are equal. We have, however, analyzed each of these structures upon assuming different gaps, searching for a minimum of the free energy in the two-or three-dimensional parameter space of gaps. In most cases, the deepest minimum is actually obtained by simply eliminating one of the orbits from the configuration (i.e. let ∆ = 0 for that orbit); the resultant structure with one less orbit appears as another structure in Table 1. For example, the bicapped square antiprism has two orbits: the first is the set of eight q's forming a square antiprism, the second is the antipodal pair of q's forming the two "caps" of the structure. Denote the gaps corresponding to these two orbits as ∆ 1 and ∆ 2 . This structure is overcrowded with intersecting rings, so it is not surprising to find that a lower-energy configuration is obtained by simply letting ∆ 2 = 0, which gives the "uncapped" square antiprism. Configurations with fewer orbits are generally more favorable, with only three exceptions known to us: the trigonal bipyramid is favored over the triangle or the antipodal pair; the square pyramid is favored over the square or the point; and the capped trigonal antiprism is favored over any of the structures that can be obtained from it by removing one or two orbits. For these configurations, Table 1 lists the results for ∆ 1 = ∆ 2 (= ∆ 3 ); the numbers can be slightly improved with ∆ 1 = ∆ 2 ( = ∆ 3 ) but the difference is unimportant.
For some configurations in Table 1 the positions of the points are completely fixed by symmetry while for others the positions of the points can be varied continuously, while still maintaining the point-group symmetry of the structure. For example, with the square pyramid we can vary the latitude of the plane that contains the inscribed pyramid base. Similarly, with the various polygonal prism and antiprism structures (and associated cappings and augmentations), we can vary the latitudes of the inscribed polygons (equivalently, we can vary the heights of these structures along the principal symmetry axis). For each structure that has such degrees of freedom, we have scanned the allowed continuous parameter space to find the favored configuration. Table 1 then shows the results for this favored configuration, and the latitude angles describing the favored configuration are given as footnotes. However, if the structure always has overlapping rings regardless of its deformation, then either no favorite configuration exists or the favorite configuration is a degenerate one that removes the overlaps by changing the structure. There are two instances where this occurs: the capped square prism can be deformed into a square pyramid by shrinking the height of the square prism to zero, and the bicapped square antiprism can be deformed into an antipodal pair by moving the top and bottom square faces of the antiprism to the north and south poles, respectively. For these structures, Table 1 just lists results for an arbitrarily chosen nondegenerate configuration.
A typical parameter scan is shown in Fig, 11, where we have plotted γ for the square antiprism as a function of the polar angle θ of the top square facet (the polar angle of the bottom square facet is π − θ). As we expect, γ is very large in regions where any rings intersect, and we search for a minimum of γ in the region where no rings intersect. The plot has a rather complicated structure of singularities and discontinuities; these features are analogous to those of Fig. 9, and as there they arise as a result of the double limit we are taking. Primary singularities occur at critical angles where pairing rings are mutually tangent on the Fermi surface. Secondary singularities occur where rings corresponding to harmonic q's, obtained by taking sums and differences of the fundamental q's that define the crystal structure, are mutually tangent. Such q's arise in the calculation of J and K because these calculations involve momenta corresponding to various diagonals of the rhombus and hexagon in Fig. 7.
In addition to varying the latitudes of the Föppl planes in various structures, we varied "twist angles". For example, we explored the continuous degree of freedom that turns a cube into a square antiprism, by twisting the top square relative to the bottom square by an angle φ ranging from 0 • to 45 • . In Fig. 12, we show a parameter scan in which we simultaneously vary the twist angle φ and the latitudes of the square planes in such a way that the scan interpolates linearly from the cube to the most favorable square antiprism of Table 1. In this parameter scan, we find a collection of secondary singularities and one striking fact: γ is much more negative when the twist angle is zero (i.e. for the cube itself) than for any nonzero value. For the cube, γ = −459.2/δµ 4 , whereas the best one can do with a nonzero twist is γ = −64.2/δµ 4 , which is the result for an infinitesimal twist angle. Thus, any nonzero twist spoils the regularities of the cube that contribute to its combinatorial advantage, and this has a dramatic and unfavorable effect on the free energy.
Finally, we have scanned the parameter space of a generic cuboid to see how this compares to the special case of a cube. That is, we vary the height of the cuboid relative to its width, without introducing any twist. This continuous variation does not reduce the combinatorial advantage of the crystal structure. As expected, therefore, we find that as long as the cuboid has no intersecting rings, it has a free energy that is very similar to that of the cube itself. Any cuboid with intersecting rings is very unfavorable. In the restricted parameter space of nonintersecting cuboidal arrangements, the cuboid with the most negative free energy is a square prism with a polar angle of 51.4 • for the top square face. (For this polar angle the pairing rings corresponding to the four corners of the square are almost mutually tangent.) This prism is slightly taller than a perfect cube, which has a polar angle of 54.7 • . The free energy coefficients of the best cuboid are β = −111.563/(δµ) 2 , γ = −463.878/(δµ) 4 . These coefficients differ by less than 1% from those for the cube, given in Table 1. There is no significant difference between the cube and this very slightly more favorable cuboid: all the qualitative arguments that favor the cube favor any cuboid with no intersecting rings equally well. We therefore expect that if we could determine the exact (rather than Ginzburg-Landau) free energy, we would find that the favored crystal structure is a cuboid with a polar angle somewhere between 51.4 • and 56.5 • , as this is the range for which no rings intersect. We expect no important distinction between the free energy of whichever cuboid in this narrow range happens to be favored and that of the cube itself.

Conclusions and Open Questions
We have argued that the cube crystal structure is the favored ground state at zero temperature near the plane-wave instability point δµ = δµ 2 . By the cube we mean a crystal structure constructed as the sum of eight plane waves with wave vectors pointing towards the corners of a cube. The qualitative points (which we have demonstrated in explicit detail via the analysis of many different crystal structures) that lead us to conclude that the cube is the winner are: • The quadratic term in the Ginzburg-Landau free energy wants a |q| such that the pairing associated with any single choice ofq occurs on a ring with opening angle ψ 0 ≃ 67.1 • on each Fermi surface.
• The quadratic term in the Ginzburg-Landau free energy favors condensation with many different wave vectors, and thus many different pairing rings on the Fermi surfaces. However, the quartic and sextic terms in the free energy strenuously prohibit the intersection of pairing rings. No more than nine rings with opening angle 67.1 • can be placed on the sphere without overlap.
• The quartic and sextic terms favor regular crystal structures, for example those that include many different sets of wave vectors whose tips form rectangles. None of the nine-wave structures with no intersections between pairing rings are regular in the required sense. The cube is a very regular eight-wave crystal structure.
Quantitatively, we find that a cube (actually, a cuboid that is only slightly taller than it is wide) has by far the most negative Ginzburg-Landau free energy, to sextic order, of all the many crystal structures we have investigated.
Let us now see what the cubic crystal structure looks like in position space. The eight q vectors are the eight shortest vectors in the reciprocal lattice of a facecentered-cubic crystal. Therefore, we find that ∆(x) exhibits face-centered-cubic symmetry. Explicitly, where the lattice constant (i.e. the edge length of the unit cube) is PSfrag replacements where the last equality is valid at δµ = δµ 2 and where ∆ 0 is the gap of the BCS phase that would occur at δµ = 0. A unit cell of the crystal is shown in Fig. 13. We have taken ∆(x) to be real by convention: we have the freedom to multiply ∆(x) by an overall x-independent phase. This freedom corresponds to the fact that the condensate spontaneously breaks the U(1) symmetry associated with conservation of quark number. Our Ginzburg-Landau analysis predicts a first-order phase transition to the cubic crystalline color superconductor at some δµ = δµ * . The fact that we predict a firstorder phase transition means that the Ginzburg-Landau analysis cannot be trusted quantitatively. Furthermore, at order ∆ 6 , which is as far as we have gone, the Ginzburg-Landau free energy for the cube is unbounded from below. We therefore have no quantitative prediction of δµ * or the magnitude of ∆. The best we can do is to note that the cube is significantly favored over the octahedron, for which the order ∆ 6 Ginzburg-Landau analysis predicts δµ * ≃ 3.6∆ 0 and predicts that at δµ = δµ 2 = 0.754∆ 0 , the gap is ∆ ≃ 0.95∆ 0 and the condensation energy is larger than that in the BCS state by a factor of about fifteen. As we have warned repeatedly, these numbers should not be trusted quantitatively: because the Ginzburg-Landau approximation predicts a strong first-order phase transition, it predicts its own breakdown. We have learned several qualitative lessons from it, however: • We have understood the qualitative reasons that make the cube the most favored crystal structure of them all.
• The Ginzburg-Landau analysis indicates that ∆, the gap parameter in the crystalline phase, is comparable to ∆ 0 , that in the BCS phase.
• We learn that the condensation energy by which the crystalline phase is favored over unpaired quark matter at δµ = 0 may even be larger than that for the BCS phase at δµ = 0. (Note that the spatial average of ∆(x) 2 in (4.1) is 8∆ 2 . It is not surprising that ∆ ∼ ∆ 0 corresponds to an |Ω min | in the crystalline phase which is larger than |Ω BCS |.) • We learn that the crystalline color superconductivity window δµ 1 < δµ < δµ * is large. Because δµ 2 is not much larger than δµ 1 , the window δµ 1 < δµ < δµ 2 wherein the single plane-wave condensate is possible is narrow. We have learned, however, that δµ * ≫ δµ 2 . Furthermore, because the condensation energy of the crystalline phase is so robust, much greater than that for the single plane wave and likely comparable to that for the BCS phase, the value of δµ 1 , the location of the transition between the crystalline phase and the BCS phase, will be significantly depressed.
Much remains to be done: • Our Ginzburg-Landau analysis provides a compelling argument that the crystalline color superconductor is face-centered-cubic. Given that, and given the prediction of a strong first-order phase transition, the Ginzburg-Landau approximation should now be discarded. What should be taken from our work is the prediction that the structure (4.1) is favored. Although ∆ could be estimated by going to higher order in the Ginzburg-Landau approximation, a much better strategy is to do the calculation of ∆ upon assuming the crystal structure (4.1) but without requiring ∆ to be small. One possible method of analysis, along the lines of that applied in Refs. [22,27,15] to the single plane wave, would be to find a variational wave function that incorporates the structure (4.1), vary within this ansatz, and find ∆ and Ω min . Another possibility is to assume (4.1) and then truncate the infinite set of coupled Nambu-Gorkov gap equations without assuming ∆ is small. This strategy would be the analogue of that applied to the chiral crystal in Ref. [31]. A third possibility is to analyze the state with crystal structure (4.1) using methods like those applied to the two plane-wave (antipodal pair) crystal structure in Ref. [41]. Regardless of what method is used, because ∆ = 0 a complete calculation would require the inclusion of condensates with q's corresponding to higher harmonics of the fundamental q's that define the crystal structure.
• The cubic crystal structure (4.1) breaks x-, y-and z-translational invariance, and therefore has three phonon modes. With the crystal structure (4.1) in hand, it would be very interesting to use the methods developed in Refs. [17] to derive the dispersion relations for these phonons and compute the parameters of the low energy effective Lagrangian describing these phonons.
• How could the "crystalline superfluid" state be detected in an ultracold gas of fermionic atoms? The simplest idea would be to look for the periodic modulation of the atom density, which goes like ∆(x) 2 µ 2 . This means that the ratio of the magnitude of the spatial modulation of the density to the density itself is of order ∆(x) 2 /µ 2 . Since the density modulation is only sensitive to ∆(x) 2 , it will not see the sign of ∆ in Fig. 13. That is, the crystal structure for ∆ 2 is simple cubic, with a unit cell of size a/2.
• In the QCD context, we need to do a three-flavor analysis, beginning with the unpaired quark matter of (1.2) and letting ud, us and us crystalline condensates form. Given the robustness of the crystalline color superconducting condensates we have found, we expect the three-flavor crystalline color superconductor gap, critical temperature and condensation energy to be comparable to or even larger than those of the CFL phase. This means that just as the crystalline color superconductor will push δµ 1 down in our model, in three flavor QCD it will push the unlocking transition between the CFL and crystalline phases up to a higher density than has been estimated to date, since present estimates have assumed the condensation energy in the crystalline phase is negligible. It also means that just as in our model the crystalline window δµ 1 < δµ < δµ * is in no sense narrow, crystalline color superconductivity will be a generic feature of nature's QCD phase diagram, occurring wherever quark matter that is not color-flavor locked is to be found.
• What effects could the presence of a layer of crystalline color superconducting quark matter have on observable properties of a compact star? Could it be the place where some pulsar glitches originate, as suggested in Ref. [15]? The next step toward answering this question and toward studying the properties of glitches that may originate within crystalline color superconducting quark matter is to take the crystal structure (4.1) and rotate it. What is the structure of a rotational vortex in this phase? Are vortices pinned? It is reasonable to guess that the vortex core will prefer to live where the condensate is already weakest, thus along the line of intersection of two of the nodal planes in Fig. 13. If this intuition is correct, the vortices would be pinned. Constructing the vortices explicitly (and then calculating the pinning force) will be a challenge. Vortices are usually constructed beginning with a Ginzburg-Landau free energy functional written in terms of ∆(x) and ∇∆(x). Instead, we have constructed a Ginzburg-Landau functional written in terms of the ∆ q 's. In principle, this contains the same information. But, it is not well-suited to the analysis of a localized object like a vortex. This means that it will take some thought to come up with a sensible model within which one can construct a vortex in the crystalline condensate (4.1).
The greatest difficulties in finding the explicit form of a vortex will be in determining the details of how the vortex core interacts with the crystal structure. Unfortunately, therefore, calculating the pinning force will likely be a hard problem. As always, it should be easier to understand the physics of a vortex far from its core. The natural expectation is that far from its core, a vortex will be described simply by multiplying the ∆(x) of (4.1) by exp[iθ(x)], where θ(x) is a slowly varying function of x that winds once from 0 to 2π as you follow a loop encircling the vortex at a large distance. In a uniform superfluid, this slowly varying phase describes a particle-number current flowing around the vortex. Here, the nodal planes make this interpretation difficult. Thus, even the long distance vortex physics will be interesting to work out.
• Our analysis does not apply to QCD at asymptotically high density, where the QCD coupling becomes weak. In this regime, quark-quark scattering is dominated by gluon exchange and because the gluon propagator is almost unscreened, the scattering is dominated by forward scattering. This works in favor of crystalline color superconductivity [18], but it also has the consequence of reducing q 0 /δµ and hence reducing ψ 0 . The authors of Ref. [18] find q 0 /δµ reduced almost to 1, meaning ψ 0 reduced almost to zero. However, the authors of Ref. [20] find q 0 /δµ ≃ 1.16 at asymptotically high density, meaning that ψ 0 ≃ 61 • . If the opening angle of the pairing rings on the Fermi surface does become very small at asymptotic densities, then the crystal structure there is certain to be qualitatively different from that which we have found. At present, the crystal structure at asymptotic densities is unresolved. This is worth pursuing, since it should ultimately be possible to begin with asymptotically free QCD (rather than a model thereof) and calculate the crystal structure at asymptotic density from first principles. (At these densities, the strange quark mass is irrelevant and a suitable δµ would have to be introduced by hand.) Although such a first-principles analysis of the crystalline color superconducting state has a certain appeal, it should be noted that the asymptotic analysis of the CFL state seems to be quantitatively reliable only at densities that are more than fifteen orders of magnitude larger than those reached in compact stars [42]. At accessible densities, models like the one we have employed are at least as likely to be a good starting point.
Notice that these vectors are the coordinates of vertices in the rhombus and hexagon shapes of Fig. 7. In particular, (k 1 k 2 ) and (ℓ 1 ℓ 2 ) are the pairs of endpoints for the solid and dashed diagonals of the rhombus figure, while (k 1 k 2 k 3 ) and (ℓ 1 ℓ 2 ℓ 3 ) are the triplets of vertices of the solid and dashed triangles in the hexagon figure.
Finally, we must integrate the Feynmann parameters. For the J integral, two of the integrals can be done using the delta functions, a third can be done analytically, and the final integration is done numerically, using an integration contour that avoids the singularity at |k − ℓ| = δµ. For the K integral, we do the x 3 and y 3 integrals using the delta functions, and then make a linear transformation of the remaining integration variables x 1 , x 2 , y 1 , y 2 , introducing new variables with a ij , b ij and c i chosen such that |k − ℓ| 2 = r 2 1 + r 2 2 + r 2 3 . (A.10) While such a transformation puts the integrand in a convenient simple form, it complicates the description of the integration region considerably. Therefore we use a Fourier-Motzkin elimination procedure [43] to express the four-dimensional integration region as a sum of subregions, for each of which we have an iterated integral with "affine" limits of integration: For each subregion A, we can immediately do the r 4 integration since the integrand is independent of r 4 . We are left with a three-dimensional integral over the volume of a polyhedron with six quadrilateral faces. We then apply the divergence theorem to turn the three-dimensional integral into a sum of surface integrals over the faces. For each surface integral, we convert to plane polar coordinates (ρ, φ) so that |k − ℓ| 2 = ρ 2 + d 2 , where d is the distance from the origin (r 1 , r 2 , r 3 ) = (0, 0, 0) to the plane of integration. Now the φ integration can be done because the integrand is independent of φ. Finally, the ρ integration is done numerically, using a deformed integration contour that avoids the singularity at |k − ℓ| = δµ.