Chiral Soliton Lattice turns into 3D crystal

Chiral perturbation theory predicts the chiral anomaly to induce a so-called Chiral Soliton Lattice at sufficiently large magnetic fields and baryon chemical potentials. This state breaks translational invariance in the direction of the magnetic field and was shown to be unstable with respect to charged pion condensation. Improving on previous work by considering a realistic pion mass, we employ methods from type-II superconductivity and construct a three-dimensional pion (and baryon) crystal perturbatively, close to the instability curve of the Chiral Soliton Lattice. We find an analogue of the usual type-I/type-II transition in superconductivity: Along the instability curve for magnetic fields $eB>0.12\, {\rm GeV}^2$ and chemical potentials $\mu<910\, {\rm MeV}$, this crystal can continuously supersede the Chiral Soliton Lattice. For smaller magnetic fields the instability curve must be preceded by a discontinuous transition.

A Fourier transform of the 2D Abrikosov lattice 29 1 Introduction and conclusion

Context
Type-II superconductors can admit magnetic flux in the form of a flux tube lattice.Magnetic flux tube lattices have been predicted a long time ago by Abrikosov within Ginzburg-Landau theory [1] and are routinely observed in superconductors in the laboratory [2,3].They may also play an important role in high-energy systems governed by Quantum Chromodynamics (QCD).For instance, they have been suggested to exist in the interior of neutron stars in the form of a proton superconductor in nuclear matter [4][5][6] or a color superconductor in quark matter [7][8][9], and may be found in the QCD phase diagram at nonzero isospin chemical potential in the form of a charged pion condensate [10].It was pointed out in our previous work [11] that elements of type-II superconductivity also apply to charged pion condensation without isospin chemical potential, but in the presence of a baryon chemical potential, which couples to the magnetic field and the pions through the chiral anomaly.In that study, we restricted ourselves for simplicity to the limit of zero pion mass and found a two-dimensional Abrikosov lattice, as in ordinary superconductors.In this paper, we generalize this study to the case of a physical pion mass and construct the resulting three-dimensional lattice.

Background and main idea
The background for our study is the so-called Chiral Soliton Lattice (CSL).It was pointed out within chiral perturbation theory that the CSL phase is the ground state of QCD in a certain region of the phase diagram spanned by the magnetic field B and the baryon chemical potential µ [12] (see also studies in the chiral limit using the gauge-gravity duality [13][14][15]).The crucial ingredient for this observation is the chiral anomaly [16,17], implemented via a Wess-Zumino-Witten (WZW) term [18,19].This term results in an energy gain from a spatially varying neutral pion condensate.Fig. 1 shows the phase transition line between the vacuum and CSL.At this phase transition it becomes energetically favorable to create a single domain wall with nonzero baryon number [20].Close to this transition line, CSL can be considered as a stack of well-separated domain walls.When they start to overlap they can no longer be described as independent domain walls -the exact CSL solution needs to be considered.In the CSL phase, the magnetic field is homogeneous, say ⃗ B = B⃗ e z .In Fig. 1, we have generalized this to the average magnetic flux per area ⟨B⟩ through the transverse x-y plane, kept constant along z.This allows for inhomogeneous magnetic fields and is necessary for our 3D lattice.The discussion of the CSL phase has been extended to nonzero temperatures [21], next-to-leading order contributions have been included [22], and its dynamical formation was studied [23].Here we do not consider any of these extensions but focus on the instability of CSL with respect to charged pion condensation.This instability was already pointed out in the original work [12], without constructing the state that supersedes CSL beyond this instability.In Fig. 1, the corresponding instability curve (blue) is labelled by B c2 , borrowing the terminology from the second critical magnetic field in type-II superconductivity.In the chiral limit, i.e., setting the pion mass to zero, m π = 0, it was shown that the system forms magnetic flux tubes arranged in a hexagonal lattice in the transverse plane with respect to the magnetic field [11].This case is particularly simple because the system does not break translational invariance in the z-direction.Besides the magnetic field, also baryon number is concentrated in the flux tubes, which renders the system a 2D baryon crystal.The construction of this crystal was performed by solving the coupled equations of motion for the electromagnetic gauge field and the neutral and charged pion condensates in an expansion close to the critical field B c2 , following the Ginzburg-Landau analysis of type-II superconductivity [1,24,25].In this paper, we shall use the same idea.However, the calculation is more involved because in the presence of a nonzero pion mass already the CSL phase itself breaks translational invariance in the direction of the magnetic field.Together For large magnetic fields our perturbative analysis predicts a continuous transition to a three-dimensional charged pion lattice, indicated schematically by the shaded region.We find that this behavior changes at a point reminiscent of the type-I/type-II transition in ordinary superconductors (blue star): As one leaves the CSL region across the dashed segment of the B c2 curve, our charged pion lattice is not the solution to the CSL instability.This includes the limit of a single domain wall (small blue dot at e⟨B⟩ = 3m 2 π ).The diagram is obtained for realistic pion mass m π = 140 MeV and pion decay constant f π = 92 MeV.In this case, the type-I/type-II transition point sits at µ ≃ 910 MeV and e⟨B⟩ ≃ 0.12 GeV 2 ≃ 6.3 m 2 π ≃ 6.1 × 10 18 G.For decreasing (increasing) pion mass, this point moves towards larger (smaller) µ.The result suggests additional phase transitions not captured by our approximation, see conjectured phase structure in Fig. 2.
with an Abrikosov lattice of charged pions in the transverse directions, the resulting phase can be expected to be a 3D crystal.Indeed, this is what we shall find: We will solve the (expanded) equations of motion in Fourier space to construct the crystal, where magnetic field, baryon number, and all pion condensates vary periodically in all three dimensions.

Main results and discussion
The properties and precise structure of our 3D crystal will be discussed in the main part of the paper.Here we focus on the main results for the phase diagram.
We find that our calculation provides us with a consistent solution only on one side of the instability curve B c2 , either outside or inside the CSL phase.If the solution exists inside CSL, we find that our 3D lattice is energetically disfavored compared to CSL, if it is outside we find that it is favored.Importantly, there is a flip from outside (large e⟨B⟩, solid segment of B c2 ) to inside (small e⟨B⟩, dashed segment of B c2 ), marked by a star in Fig. 1.As a consequence, as one moves across the solid segment of the B c2 line, a continuous transition from CSL to the 3D crystal is predicted, as indicated by the shaded region in the phase diagram.This is analogous to a transition from the normal-conducting state to the flux tube lattice in a type-II superconductor.In contrast, our lattice does not provide a solution for a continuous transition across the dashed segment.
At the endpoint of the B c2 curve (small blue dot) CSL is identical to a single domain wall.Therefore, our solution along B c2 -although not preferred everywhere -interpolates from a 2D lattice with no dependence on z (chiral limit of Ref. [11], valid for asymptotically large B c2 ) to a lattice enclosed by a single domain wall localized in the z-direction.In particular, our perturbative solution does not yield a stable lattice in the domain wall.
The point marked by the star is reminiscent of the transition from type-II to type-I superconductivity.However, there is an important difference to ordinary superconductivity to be kept in mind.In our pionic system there is no equivalent of a Meissner phase.The would-be Meissner phase is a homogeneous charged pion condensate that completely expels the magnetic field.In the absence of a magnetic field, however, there is no CSL phase and thus no neutral pion condensate.Therefore, since we consider zero isospin chemical potential, there is no effective potential for the charged pions to condense in the first place.
We cannot rigorously claim to have found the ground state of chiral perturbation theory; let alone the ground state of QCD, due to the limitations of chiral perturbation theory.Nevertheless, it is tempting to speculate what our findings suggest for the QCD phase structure.In Fig. 2 we present a phase diagram with conjectured phase transition lines (red), based on the following observations.Firstly, if one reaches the dashed segment of B c2 from within CSL, one cannot transition to our charged pion lattice and yet we know that CSL becomes unstable there.Therefore, there must be an earlier transition either to the pion lattice or some other phase not considered in our calculation.This transition is necessarily discontinuous since it is not indicated by a CSL instability.It can connect to the continuous transition at any point on the solid B c2 segment.Secondly, we know that at zero magnetic field there is a first-order transition from the vacuum to homogeneous nuclear matter at µ 0 ≃ 922.7 MeV 1 .Consequently, there must be at least one more transition, namely from homogeneous nuclear matter to the pion (and baryon) crystal.While the existence of such a liquid-solid transition seems inevitable, its exact location in our sketch is purely speculative.Our conjecture shows the simplest possible scenario to accommodate our findings with the existence of homogeneous nuclear matter at zero magnetic field.The phase structure may, however, be more complicated.For instance, we have ignored superfluidity and superconductivity in homogeneous nuclear matter, which may break down as the magnetic field is increased or which may turn via a single phase transition into the superconducting baryon crystal.(The critical magnetic field for neutron and proton Cooper pairing depends on density with maximal values of the order of 10 17 G ∼ m 2 π [28,29].)

Outlook
To put our results into context we should first recall that they are obtained within leadingorder chiral perturbation theory.Going beyond leading order requires the calculation of 1 Interestingly, by only slightly varying the pion mass and the pion decay constant, using mπ = 139 MeV and fπ = 92.4MeV, we find that the chemical potential for our type-I/type-II transition point is µ ≃ 922 MeV, strikingly close to µ0.However, there seems to be no reason why our transition point should reproduce the baryon onset µ0, mainly because it is located at a nonzero magnetic field, where the baryon onset is likely to be different from µ0 [14,26].Also, by decreasing mπ at fixed fπ, our transition point moves towards larger chemical potentials.Since we do not expect the baryon mass to increase accordingly [27], the transition point does not seem to be coupled to the baryon onset for all pion masses.  1 plus the fact that at zero magnetic field there is a transition from the vacuum to homogeneous nuclear matter at µ 0 ≃ 922.7 MeV.Our type-I/type-II transition point (pale blue star) sets a limit for where the (red) discontinuous transition can attach to the (blue) continuous transition between CSL and the 3D baryon crystal.Pale segments are taken over from Fig. 1 but do not correspond to actual transitions in our conjecture.
the excitations of our lattice, which would be an interesting problem for future studies.In the regime where our calculation predicts a 3D pion crystal, the energy scales set by the magnetic field and/or the chemical potential become comparable to 4πf π and thus are sufficiently large for chiral perturbation theory itself to approach its limit of validity.Hence we do not interpret our results as rigorous QCD predictions; ultimately, it would be interesting to compare our predictions with a first-principle QCD calculation.However, due to the presence of the chemical potential this seems to be out of reach for current methods in lattice gauge theory.Therefore, it would also be interesting to search for our 3D crystal within suitable models, such as nucleon-meson models or employing the gaugegravity duality.
To explore our conjectured discontinuous phase transition lines, it would be necessary to go beyond our perturbative approach, which is only valid close to a continuous transition.Possibly this can be done along the lines of previous studies [30][31][32][33][34][35], which have constructed baryon crystals without using an expansion at the CSL instability, but which have made other simplifications compared to our study.Reference [35], which we learned of after completing the work reported here, predicts a discontinuous transition from the CSL phase to a phase made of domain-wall Skyrmion chains, however without allowing the magnetic field to be dynamical and without interaction of the chains in the transverse plane.It would be interesting to compare and combine these results with ours -which are fully dynamical in longitudinal and transverse directions, including the magnetic field, but only valid at a continuous transition.
Since there exist various extensions and generalizations of the CSL phase, one can also think about applying these extensions to our crystal, for instance including a nonzero temperature [21] or the η meson [36], or replacing the magnetic field with an externally imposed rotation [37][38][39].Especially going beyond zero temperature might be crucial to connect our results to real-world systems.Large magnetic fields are reached in neutron stars and mergers thereof, as well as in heavy-ion collisions.In neutron star mergers and heavy-ion collisions, however, the zero-temperature approximation is not applicable.Also the magnetic field strengths where our crystal is stable are most likely pushing the limits of these systems.Therefore, it is of particular interest to find the critical magnetic field for our conjectured transition of liquid nuclear matter to the baryon crystal.
Including an isospin chemical potential is another such extension.An isospin chemical potential gives rise to charged pion condensation more directly, without the detour via the chiral anomaly and neutral pion condensation.In the absence of a baryon chemical potential, a standard 2D vortex lattice can thus be induced by a magnetic field [10,40].It was later pointed out that the chiral anomaly allows for a CSL phase in this case as well, and its competition with the vortex lattice was studied [41] (however, without taking into account the possibility of a 3D lattice of coexisting neutral and charged pion condensates).It would be interesting to connect those findings to our results, i.e., to ask whether and how our anomaly-induced superconductivity relates to the more conventional superconductivity already present there.As a first step, this can be done by switching on an isospin chemical potential in our perturbative solution of the 3D crystal, or by searching for an instability of the isospin analogue of the CSL phase towards charged pion condensation.

Structure of the paper
The preparatory sections 2 and 3 serve to introduce our setup within chiral perturbation theory (Sec.2) and to discuss the CSL phase and its instability (Sec.3), thereby establishing our notation and collecting the necessary equations for the subsequent sections.Our main calculation, leading to the main results discussed in Sec.1.3, is laid out in Secs. 4 and 5.This contains the explanation of our expansion in Sec.4.1 and the construction and discussion of the 3D lattice in Secs.4.2 -4.5, as well as the calculation of the free energy in Secs.5.1 -5.4.The results that feed directly into the phase diagram in Fig. 1 can be found at the end of Sec.5.4.We work in natural units ℏ = c = k B = 1 and use Heaviside-Lorentz units for the electromagnetic gauge field, in which the elementary charge is e ≃ 0.30282.

Lagrangian
In this section, we write down our Lagrangian and collect the most important definitions.Since the starting point and notation are the same as in Ref. [11], we refer the reader to this reference for additional details.We work within leading-order, two-flavor chiral perturbation theory with a dynamical electromagnetic field and a WZW term, such that the Lagrangian can be written as The electromagnetic part is where is the electromagnetic field strength tensor, with the electromagnetic U (1) gauge field A µ .Following the parameterization of Ref. [11] (see also Ref. [12]), we write the pionic contribution as where is the covariant derivative.The real scalar field α and the complex scalar field φ parameterize the neutral and charged pions in a convenient way.In the most conventional parameterization, the chiral SU (2) field is written as Σ = (σ + iπ a τ a )/f π , where a = 1, 2, 3, and τ a are the Pauli matrices.In the (σ, π 3 ) sector we introduce α as the polar angle and eliminate the modulus with the help of the constraint f 2 π = σ 2 + π a π a .In the (π 1 , π 2 ) sector the complex field φ is obtained from a rotation by the angle α of the field (π 1 + iπ 2 )/ √ 2. The WZW term can be written as where A B µ = (µ, 0, 0, 0) contains the baryon chemical potential µ, and the Goldstone-Wilczek baryon current [12,20,42] can be written as [11] with the non-anomalous contribution to the charge current (2.7) Since we will only consider equilibrium configurations, we ignore any time dependence.Also, we only consider a magnetic field, ⃗ B = ∇ × ⃗ A, neglecting any electric field contribution.This requires the assumption that the system is locally neutral due to the addition of electrons or positrons that will not appear explicitly in our calculation (analogous to the usual Ginzburg-Landau study of superconductivity).We also perform our analysis on the classical level without taking into account fluctuations of the pions and the photon.Then, the local free energy density is simply given by Ω = −L and can be written as where is the local baryon number density.The free energy density (2.8) can be viewed as an effective potential and we can use it to derive the equations of motion for the dynamical fields ⃗ A, α, φ.

Chiral Soliton Lattice and its instability
Next, let us recapitulate the CSL phase and its instability [12] to set the stage for our construction of the 3D pion lattice.

CSL solution
In the absence of charged pions, the effective potential is where we have subtracted the constant B = α = 0 contribution to normalize the vacuum free energy to zero.Without charged pions, α is identical to the neutral pion field π 3 , and we use the notation α 0 to distinguish the CSL solution from our full solution in the subsequent sections.The equation of motion for α 0 is Aligning the z axis with the magnetic field, ⃗ B = B⃗ e z , the physically relevant solution only depends on z.We write the solution as where sn denotes the Jacobi elliptic function, and p ∈ [0, 1] is a dimensionless number (the "elliptic modulus") whose value has to be determined by minimizing the free energy for given µ and B.Moreover, we have introduced the dimensionless variable The boundary conditions for the solution (3.3) are chosen such that α 0 winds from 0 to 2π in the unit cell z ∈ [−K(p 2 ), K(p 2 )], where K is the complete elliptic integral of the first kind.For p → 1, the unit cell becomes infinitely large; this limit corresponds to an isolated domain wall, where the phase α 0 winds around once in a spatial domain of order m −1 π , domain wall limit: α 0 (z) = 4 arctan e mπz . (3.5) In the following we denote spatial averaging over the z direction of any field f (⃗ r) by Since the CSL configuration is constant in the transverse directions, the average free energy density is given by where E is the complete elliptic integral of the second kind.Minimizing F 0 with respect to p gives which is an implicit equation for the elliptic modulus p.This yields the free energy at the minimum, We can use Eqs.(3.8) and (3.9) to obtain a simple expression for the average baryon number density, We see from Eq. (3.9) that the CSL state is favored over the vacuum for p < 1 and the continuous transition occurs at p = 1, the domain wall limit.With Eq. (3.8) and E(1) = 1, the critical magnetic field for the onset of CSL is given by the black curve in Fig. 1.We recover the chiral limit with m π → 0, p → 0 and, using Eq.(3.8) with In this limit, the CSL solution (3.3) reduces to chiral limit: Up to the irrelevant constant shift by π, this is the m π = 0 solution used in Ref. [11].While in Ref. [11] the chiral limit was employed for all magnetic fields for simplicity, it is a valid approximation for magnetic fields much larger than the pion mass squared.Our full solution will also, just like α 0 in this section, interpolate between the domain wall limit and the chiral limit.

CSL instability
The CSL phase is unstable with respect to charged pion condensation [12].To observe this instability we go back to the original Lagrangian (2.1), including the complex field φ.The equation of motion for φ * is We insert the CSL solution α 0 into this equation and linearize in the field φ.Then, within our assumptions of a locally charge neutral system, setting A 0 = 0 and e ⃗ A + ∇α 0 = eBx⃗ e y and using the ansatz φ(⃗ r, t) = e −iωt e ikyy h(x, z), with the yet to be determined function h(x, z), Eq. (3.14) becomes where we have used We factorize h(x, z) = f (z)g(x) and divide by f (z) to obtain an equation for g(x) that has the form of the Schrödinger equation for the harmonic oscillator.We can thus immediately read off the (z-dependent) energy eigenvalues labeled by l = 0, 1, 2, . .., Rearranging this equation gives with the Lamé operator where, in our case, m = 2, and The lowest eigenvalue of the m = 2 Lamé equation is [43] with corresponding eigenfunction where This factor is introduced to ensure that f 0 is normalized, which is convenient, although not strictly necessary, for our calculation.With Eq. (3.20) and setting l = 0, the lowest-energy excitation is The instability with respect to pion condensation sets in at ω 2 = 0, i.e., at a critical field given by Our notation B c2 anticipates that this critical field is reminiscent of the second critical magnetic field of a type-II superconductor.For the following it is useful to keep in mind that Eq.
Again, it is useful to check the domain wall and chiral limits.For the domain wall, p = 1, we easily obtain from Eq. (3.26), domain wall: in accordance with Ref. [20].In the chiral limit, p → 0, we find with Eq. (3.12), chiral limit: in accordance with Refs.[11,12].The instability curve given by Eq. (3.27) is shown as the blue line (solid and dashed segments together) in Fig. 1 for the physical values of m π and f π .The solution to the instability, at least close to the instability curve, is the main topic of this paper, to which we turn now.
4 Calculation of the 3D pion crystal

Expansion at the critical magnetic field
At the second critical magnetic field of a type-II superconductor, there is a continuous transition from the normal-conducting state, where the magnetic field penetrates homogeneously, to the superconducting state, where the magnetic field is confined into flux tubes.Magnetic flux tubes arise from a periodically varying charged condensate, which is zero along the central axes of the flux tubes.At the transition, while the lattice constant of the flux tube lattice is nonzero and finite, the amplitude of the condensate is infinitesimally small.This allows for a small-field expansion, which we will also employ at our critical field B c2 .For a consistent power counting, let us introduce the smallness parameter Here, as already introduced in the phase diagram in Fig. 1, ⟨B⟩ denotes the average magnetic flux per area through any z = const plane.We will require this average flux to be the same at each z and use ⟨B⟩ as our externally given thermodynamic variable2 .For the CSL phase of the previous section, ⟨B⟩ = B. We expand in ϵ at fixed baryon chemical potential µ.Therefore, due to the two-valuedness of B c2 as a function of µ in a certain regime (see Fig. 1) the region where we expect an Abrikosov lattice can either be at ⟨B⟩ > B c2 or at ⟨B⟩ < B c2 , hence the modulus in Eq. (4.1).We denote the expansion of our dynamical fields by where The expansion of the gauge field implies for the magnetic field with ⃗ B 0 = B c2 ⃗ e z of order 1 and δ ⃗ B ∼ ϵ 2 .We can now write down the equations of motion derived from the effective potential (2.8) order by order in ϵ.The leading and next-toleading order equations for each field are [11] D 0 φ 0 = 0 , (4.4a) ) ) where we have abbreviated The equation of motion (4.4c) is solved by the CSL solution, as discussed in Sec.3.1, while Eq.(4.4e) is trivially solved by the constant leading-order field ⃗ B 0 = B c2 ⃗ e z .The other equations will now be discussed step by step, starting with Eq. (4.4a).

Compute φ 0
With e ⃗ A 0 + ∇α 0 = eB c2 x⃗ e y and Eq.(3.16), we write Eq. (4.4a) as This equation is solved by a product ansatz, where the z-dependence is given by the wave function of the state that becomes unstable, with f 0 (z) from Eq. (3.22).With B c2 from Eq. (3.26) this yields This is the usual differential equation obtained at the second critical field of a type-II superconductor [1,11,24], such that we can immediately write down the solution where C n are complex coefficients, q is the wave number, and is the coherence length.We shall only consider solutions that have discrete translational invariance and thus produce a regular 2D Abrikosov lattice, namely with C ∈ C and |C| to be determined.The unit cell in the x-y-plane is given by the rectangle with side lengths 12) The lattice structure can be varied continuously by changing q.It is convenient to work instead with the dimensionless parameter Restricting ourselves to a ∈ [0, 1], we have a quadratic lattice for a = 1 and a hexagonal lattice for a = 1/ √ 3. (The range a ∈ [1, ∞] gives the same structures with x and y exchanged.)The hexagonal structure minimizes the free energy in a standard type-II superconductor [24].We shall see later that this is still the case in our system, despite the additional lattice structure in the z direction.For the following we introduce the notation for averaging any field f (⃗ r) over the x-y plane, In order to solve the remaining equations of motion and to compute the free energy, it is convenient to apply a Fourier decomposition.Eventually, we will need to go to Fourier space with respect to all three spatial dimensions.Here we start by a two-dimensional Fourier series for the 2D Abrikosov lattice given by ϕ 0 .We will need the Fourier series of We denote the Fourier transform by ω(n, m) = ω( ⃗ k ⊥ ) with where ⃗ k ⊥ = (k x , k y ) with Hence, the sum over ⃗ k ⊥ stands for a double sum over n, m and we have the orthogonality relation

.19)
In Appendix A we compute the Fourier transform of ω(x, y) explicitly, with the result (obtained already in Ref. [24] for the hexagonal lattice) Here, the ± originates from the two possible sign choices in Eq. (4.11).This choice plays no role for any of the physics and thus we may simply continue with the upper sign, in which case we have ω( ⃗ k ⊥ ) = ω(− ⃗ k ⊥ ).We can use this Fourier transform to compute where we have abbreviated

.22)
The quantity β was introduced in Abrikosov's original paper [1].It only depends on a and plays an important role in determining the preferred lattice structure.With the help of the Fourier transform of ω, the calculation of β is particularly simple.For a direct calculation in position space see for instance Appendix A in Ref. [11], where the result is expressed in terms of the Jacobi theta function instead of the double sum of Eq. (4.21).We will use β also in our final result for the free energy but will find that, in contrast to the standard Ginzburg-Landau superconductor, additional functions of a occur.
In the following, we will need the full 3D Fourier transform of the charged condensate squared |φ 0 (⃗ r)| 2 .Having discussed the transverse directions, we now turn to the longitudinal direction and denote s(z) ≡ f 2 0 (z) .(4.23) We write the Fourier series in the z direction as where such that the sum over where we have introduced the dimensionless wave number kz = pk z /m π , and the orthogonality relation is ⟨e Since s is even in z we can work with the somewhat simpler form Here we have used ŝ(0) = 1 due to the normalization (3.24).All other Fourier components have to be computed numerically from Eq. (4.26).

Compute δB
The goal of this section is to calculate the correction to the magnetic field due to the charged pion condensate.Since the charged pions couple to the neutral pion field, which, in turn, is modulated in the z direction already in the pure CSL phase, we can expect a modulation of the magnetic field in the z direction as well.Because of the Maxwell equation ∇• ⃗ B = 0, this inevitably creates nonzero (and non-constant) x and y components of the magnetic field.These expectations are borne out in the following solution to Eq. (4.4f).The components of this vector equation are We solve these coupled partial differential equations in Fourier space.Separating a term linear in x to account for the correct boundary conditions of the magnetic field and combining the longitudinal and transverse Fourier series discussed above to a 3D Fourier series, we write where c is a constant, to be determined later, and where ⃗ k = ( ⃗ k ⊥ , k z ).Employing Coulomb gauge, ∇ • δ ⃗ A = 0, Eqs.(4.29) can now easily be solved in Fourier space to obtain As a consistency check, one confirms that the gauge condition is obeyed in Fourier space, k x δ Âx + k y δ Ây + k z δ Âz = 0.With the help of Eq. (4.30) we can write the magnetic field as where Consequently, the magnetic field components in Fourier space are Inserting this expression for c, we can write the magnetic field components in position space as From top left to bottom right, the elliptic modulus is p = 0.1, 0.5, 0.9, 0.99, which gives a sequence from (near) the chiral limit to (near) the limit of a single domain wall.For all plots, we have set ⟨B⟩/(eω 0 ) = 1, which exaggerates the effect of deviations from a constant magnetic field for illustrative purposes.Our expansion assumes these deviations to be small.
In contrast to the lowest-order contribution to the charged pion condensate φ 0 (⃗ r), the magnetic field correction δ ⃗ B(⃗ r) does not factorize into fields depending on longitudinal and transverse directions separately.With the help of Eq. (4.35) we can evaluate ⃗ B/⟨B⟩ numerically after choosing values for a, ⟨B⟩/(eω 0 ), and p.In Fig. 3 we plot the result for the hexagonal lattice a = 1/ √ 3 and ⟨B⟩/(eω 0 ) = 1.Since ⃗ B/(eω 0 ) is of order 1/ϵ 2 , the chosen value is somewhat small, and our approximation is strictly speaking questionable for this value.The benefit of this choice is that the deviations of ⃗ B from the constant field ⃗ B 0 become more visible.The final parameter, the elliptic modulus p, is used to show different scenarios in the range between the chiral limit and the domain wall.To this end, we have chosen four different values of p.In the upper left panel we show a case close to the chiral limit, which reproduces the result of Ref. [11]: The magnetic field is confined in flux tubes and (nearly) constant in the z direction, as in an ordinary type-II superconductor.The distance between the flux tubes is set by the scale ξ.As we increase the elliptic modulus we see that the z-dependence becomes non-trivial, with the magnetic field enhanced in bubble-like regions.The lower right panel is close to the domain wall configuration, where the layers in the z direction are well separated and the magnetic field between them is nearly constant.

Compute δα
Next, we turn to Eq. (4.4d) to compute δα.We first separate a topological (= nonzero winding) contribution from the contribution that we will write as a Fourier series.For fixed µ, the CSL solution α 0 with elliptic modulus p is valid at B c2 .We are interested in the new α 0 + δα at ⟨B⟩.Therefore, we introduce the elliptic modulus p + δp, which reproduces the CSL solution at ⟨B⟩.Using Eq. (3.8) for p at B c2 and for p + δp at ⟨B⟩, we find The topological contribution to δα can now be obtained by expanding α 0 at elliptic modulus p + δp for small δp.This gives the O(1) contribution α 0 at elliptic modulus p plus an O(ϵ 2 ) contribution that we denote by α 1 δp.Together with a non-topological contribution, the total O(ϵ 2 ) correction is thus where the dimensionless factor ω0 /f 2 π of order ϵ 2 has been introduced for convenience, such that δ α is of order 1, and where with the Jacobi epsilon function E. (The derivative with respect to p is taken at fixed z, not fixed z.)With we easily verify that α 1 fulfills the differential equation We are now prepared to solve Eq. (4.4d).To this end, we first compute where This function is antisymmetric, t(z) = −t(−z), and thus its Fourier components can be computed numerically via In particular, the zero mode vanishes, t(0) = 0. Furthermore, we abbreviate the function whose Fourier components û(ℓ) we compute numerically.Inserting all this into Eq.( 4.4d), we obtain in Fourier space, Importantly, the topological contribution has dropped out due to Eq. (4.40).We have assumed the non-topological part to be odd in z, i.e., δ α(ℓ) = −δ α(−ℓ), which gives a consistent solution to the differential equation.For the practical calculation, it is useful to note that only depends on the geometric parameter a and the elliptic modulus p (besides the integers n, m, ℓ), and not explicitly on m π .Therefore, for given a and p, Eq. (4.45) yields the Fourier modes δ α(n, m, ℓ) and thus the non-topological contribution δα 1 .Due to the term with the convolution, for each transverse mode (n, m) we have to solve a coupled set of linear equations for the longitudinal Fourier modes ℓ.In the practical calculation, we must, of course, restrict ourselves to a truncation, say n, m ∈ [−n ⊥ , n ⊥ ] and ℓ ∈ [−n z , n z ].It turns out that convergence to the final result is already achieved to a very good accuracy for n ⊥ of about 3 or 4 for all quantities we consider, while it is usually sufficient to choose n z of the order of 10 (except for elliptic moduli p very close to 1, where we need up to n z ∼ 100).
We have also verified our numerical results by inserting δα 1 (⃗ r), obtained from its Fourier components, back into the differential equation (4.4d).In Fig. 4 we show δα 1 for a = 1/ √ 3 and different values of p at a fixed point (x, y) as a function of z.This plot is shown mostly for replicability of our calculation, δα 1 is one of the ingredients needed later for our physical results.
We have now solved 5 of the 6 equations of motion (4.4).The equation for δφ (4.4b) does not need to be solved explicitly for our purposes.We will make use of it implicitly in the calculation of the free energy in Sec. 5. ), corresponding to a local extremum of the 2D hexagonal lattice.In each case the curve is plotted over one period −K(p 2 ) < z < K(p 2 ).

Crystalline baryon density
With the results of the previous subsections we can now calculate the local baryon number density n B (2.9).We shall evaluate n B consistently up to order ϵ 2 .As we see from Eq. (2.9), n B has a magnetic and a vorticity contribution (this structure from the chiral anomaly reflects the well-known chiral magnetic and chiral vortical effects [44]).The magnetic contribution contains the leading-order CSL part plus ϵ 2 corrections, where the Fourier transforms δ Bz and δ α are obtained from Eqs. (4.33) and (4.45).The vorticity contribution is of order ϵ 2 , and we first note that φ, α, and ⃗ A in the three-current ⃗ j can all be replaced by their lowest-order contributions.With the Maxwell law (4.4f)we can therefore use ⃗ j = ∇ × δ ⃗ B and thus ∇ × ⃗ j = −∆δ ⃗ B in the calculation of the baryon density, which yields For a quantitative plot of the baryon density it is necessary to anticipate a result of Sec.5.3, namely Eq. (5.17), which allows us to evaluate n B for given values of a, p, m π , and the smallness parameter ϵ (4.1).
The results are shown in Fig. 5 for the physical pion mass and a hexagonal lattice.We show the baryon density in units of saturation density of symmetric nuclear matter, Here and in all following numerical results we have used the value of the pion decay constant f π = 92 MeV.The vertical direction is distinguished by the magnetic field, which points (predominantly) into the z direction, see Fig. 3.The variation of n B in the transverse direction is due to charged pions, whose vorticity contributes to the topological baryon number, in addition to a magnetic field contribution.n 0 = 0.15 fm −3 , and use physical units for the spatial directions, in particular employing Eq. (4.10) for the coherence length.The average value of the baryon number for each plot is determined by the CSL phase at the specific point (µ, eB c2 ) on the instability curve, see Eq. (3.10).Also, the CSL phase already gives rise to variations in the baryon number in the z direction.The variations in the transverse direction, however, only arise due to the presence of the charged pions as we move away from the CSL instability curve by an amount given by ϵ.We have set ϵ = 0.1 for all plots except for the bottom right one, where a larger value, ϵ = 0.3, is needed to make the lattice structure clearly visible.In the upper left plot we see that for very large average magnetic field and thus very large average baryon density, baryon number is localized in tube-like structures, elongated in the direction of the magnetic field.(For eB c2 → ∞ the result approaches the chiral limit, where the baryons are localized in tubes without any variation in z.)As we move to smaller magnetic fields, the lattice becomes coarser and the size of the baryon lumps becomes larger, approaching the size of actual baryons as we approach densities n B ∼ n 0 .The bottom right panel shows a sheet-like structure, where baryon number is confined in almost spherical bubbles within a domain wall that is well separated in the z direction from the next sheet.As the results of Sec.5.4 will show, only the first three panels of Fig. 5 depict configurations that are stable in our perturbative approach.
5 Free energy of the 3D pion crystal

General expression
After having constructed the 3D pion crystal, we need to compute its free energy to determine whether it provides a solution to the CSL instability.Before inserting our solution into the effective potential, we can bring the free energy into a convenient form in exactly the same way as for the ordinary Abrikosov lattice and for the chiral limit of Ref. [11].Starting from (2.8), we derive [11] where V is the volume of the system, where we have subtracted the vacuum contribution, as in Eq. (3.1), and where we have introduced the effective coupling Here, ⟨−⟩ denotes averaging over all three spatial directions, combining the integrals (3.6) and (4.14).In the textbook treatment of a Ginzburg-Landau superconductor, this coupling corresponds to a constant in front of the quartic term of the complex scalar field.As written, the free energy (5.1) contains terms of all orders in ϵ through α and B, while the term λ * |φ 0 | 4 is of order ϵ 4 .This expression will turn out to be useful because we are eventually interested in the O(ϵ 4 ) free energy difference between our crystal and the CSL phase.
The final ingredient for the calculation is the orthogonality relation which, again following the steps of Ref. [11], can be recast in the form where Eqs.(4.4b) and (4.4f) have been used.In the following subsections we will compute λ * , evaluate Eq. (5.4) to extract an expression for eω 0 (which corresponds to determining the coefficient of the Abrikosov lattice |C|), and then use these results to compute the free energy F (5.1).

Compute eω 0
The evaluation of Eq. (5.4) will give us an expression for eω 0 , from which the yet unknown coefficient |C| 2 can be read off.We consider the various terms in Eq.With the help of our Fourier decomposition it is then easy to see that ⟨∇δα where, in the first step, we have dropped the boundary contributions, and in the second step we have used the Fourier decomposition.For the next term we recall that δα contains two contributions, see Eq. (4.37), and we compute The second term comes from the non-topological contribution δα 1 , which vanishes at the boundaries z = ±K(p 2 ), and we have used partial integration before inserting the Fourier decomposition to recover the function t(z) from Eq. (4.42).The first term originates from the topological contribution α 1 δp, and one can perform all z integrals explicitly with the help of Eqs.(4.38) and (4.39).The result is expressed in terms of the function with N (p) from Eq. (3.23), and This function has the limits G(0) = 1 (chiral limit) and G(1) = −1 (domain wall).Finally, with Eqs.(5.7) and (5.10) we immediately obtain We now insert everything into Eq.(5.4), use the expression for c in Eq. (4.34), solve the resulting equation for eω 0 , separate the ℓ = 0 mode in the Fourier series to bring the result into a convenient form, and use the definition of β (4.21).This yields where we have introduced These two Fourier sums need to be computed numerically as functions of a, b, m π .Both of them vanish in the chiral limit.While this is obvious for H 2 due to the prefactor m 2 π , for H 1 we note that in the chiral limit, where s(z) = 1, we have ŝ(ℓ) = δ ℓ0 and thus all terms in the sum are zero.
The result (5.17) is important for the interpretation of our pion crystal.First of all, we recall ω0 = |C| 2 / √ a (4.22).The right-hand side is a function of a (characterizing the transverse lattice structure), b (dimensionless version of eB c2 ), and the pion mass m π .Moreover, the right-hand side is proportional to the difference ⟨B⟩ − B c2 , which is by definition of order ϵ 2 , while all other quantities on the right-hand side are of order 1.Hence, Eq. (5.17) gives an expression for |C| 2 and we confirm as a consistency check that |C| 2 (and ω0 ) is of order ϵ 2 .As another check we observe that in the chiral limit, with G going to 1 and both H 1 and H 2 going to zero, we exactly reproduce the result of Ref. [11], which has the same form as for a standard superconductor.
Since the left-hand side is greater than zero, our configuration is only a physical solution if the right-hand side is greater than zero as well.One easily checks numerically that G is positive if and only if p < 0.91792.This is close to the point where the CSL instability curve turns around (reaches its maximal µ, see Fig. 1).With the help of Eq. (3.27) one finds that this turning point is at p ≃ 0.92751, with p being smaller than this value towards larger eB c2 .The discrepancy between these two numerical values is due to our expansion and we will ignore the small region between them 3 .Hence, the function G takes care of the twovaluedness of the CSL instability curve.As a consequence, along the entire instability curve, our solution is valid outside the stable CSL region (⟨B⟩ > B c2 on the upper branch and ⟨B⟩ < B c2 on the lower branch) if and only if the denominator in Eq. (5.17) is positive.This denominator is not only a function of p, but also depends on a and m π .We shall encounter the same denominator in the calculation of the free energy and discuss its behavior at the end of the following subsection.

When the 3D pion crystal is preferred
We are now prepared to compute the free energy (5.1).Again, let us discuss the various terms separately.The pure magnetic field contribution is first rewritten with the help of Eq. (4.3) as where we have used ⟨δB x ⟩ x,y = ⟨δB y ⟩ x,y = 0, ⟨ ⃗ B⟩ = ⟨B⟩⃗ e z .Apart from the ⟨B⟩ 2 contribution, all terms are of order ϵ 4 .This result would not further change if we allowed for higher-order corrections to ⃗ B. With Eqs.(4.33) we find (5.20) In the chiral limit, the last term vanishes because of ŝ(ℓ) = δ ℓ0 .The terms including α require a separate discussion.So far, we expanded α = α 0 + δα and included a topological correction of order ϵ 2 in δα, together with the non-topological contribution δα 1 that was computed via its Fourier modes.Since we need the free energy up to order ϵ 4 , this is not sufficient.One can easily show that any non-topological ϵ 4 contribution that vanishes at the boundaries of the period in z does not contribute.For the topological contribution, however, higher orders do matter.Fortunately, there is an easy way to include them: In the following, α 0 is understood to be evaluated at p + δp (4.36), which is the correct elliptic modulus at ⟨B⟩.This creates contributions to all orders in ϵ.These orders will eventually be absorbed in the CSL free energy, to which we need to compare the free energy of our 3D crystal.With this interpretation of α 0 in mind and using the equations of motion for α 0 (4.4c) and δα  and including a transverse lattice confined in a domain wall.Since this is not the case, the single domain wall does not turn into a charged pion lattice, at least not continuously 4 .The type-I/type-II transition point moves towards smaller critical magnetic fields (and thus to larger µ) as we decrease the pion mass.This is demonstrated by the red and blue curves in the right panel.In particular, the blue curve is the chiral limit from Ref. [11].
The dependence of the type-I/type-II transition point on the pion mass is shown in Fig. 7.The two curves give the coordinates of the point (µ, e⟨B⟩), shown for the physical pion mass as a star in the phase diagram of Fig. 1.We see that this point approaches a finite limit for vanishing pion mass (which sits at a chemical potential far beyond the validity of chiral perturbation theory) and that for large pion masses it moves up the CSL instability line without any finite boundary.
Absorbing for now the irrelevant details in the function F (n 1 ), we compute by splitting the sum into even and odd contributions, for even n where, in the second step, we have used the form of the coefficients C n given in Eq. (4.11).

2 ]Figure 1 .
Figure1.Phase diagram in the plane of magnetic field and baryon chemical potential, showing our main results along the CSL instability curve B c2 (blue).For large magnetic fields our perturbative analysis predicts a continuous transition to a three-dimensional charged pion lattice, indicated schematically by the shaded region.We find that this behavior changes at a point reminiscent of the type-I/type-II transition in ordinary superconductors (blue star): As one leaves the CSL region across the dashed segment of the B c2 curve, our charged pion lattice is not the solution to the CSL instability.This includes the limit of a single domain wall (small blue dot at e⟨B⟩ = 3m 2 π ).The diagram is obtained for realistic pion mass m π = 140 MeV and pion decay constant f π = 92 MeV.In this case, the type-I/type-II transition point sits at µ ≃ 910 MeV and e⟨B⟩ ≃ 0.12 GeV 2 ≃ 6.3 m 2 π ≃ 6.1 × 10 18 G.For decreasing (increasing) pion mass, this point moves towards larger (smaller) µ.The result suggests additional phase transitions not captured by our approximation, see conjectured phase structure in Fig.2.

2 ]Figure 2 .
Figure 2. Conjectured phase diagram, based on the results of Fig.1plus the fact that at zero magnetic field there is a transition from the vacuum to homogeneous nuclear matter at µ 0 ≃ 922.7 MeV.Our type-I/type-II transition point (pale blue star) sets a limit for where the (red) discontinuous transition can attach to the (blue) continuous transition between CSL and the 3D baryon crystal.Pale segments are taken over from Fig.1but do not correspond to actual transitions in our conjecture.
(3.26)  gives us a one-to-one relation between the critical magnetic field and the elliptic modulus p ∈ [0, 1].The CSL instability line is parameterized by p, from p = 0, corresponding to infinitely large b, to p = 1, corresponding to the end of the instability line at b = 3, where CSL has turned into a single domain wall.The corresponding chemical potential at the onset of the instability is then obtained from Eq. (3.8) as a function of b (or p), .33c) To implement the boundary condition, we use Eq.(4.3) to compute ⟨B⟩ ≡ ⟨B z ⟩ x,y = c + B c2 − eω 0 , i.e., c = ⟨B⟩ − B c2 + eω 0 .(4.34)

Figure 3 .
Figure 3. Dimensionless magnetic field ⃗ B/⟨B⟩ from Eq. (4.35) for the hexagonal lattice a = 1/ √ 3.From top left to bottom right, the elliptic modulus is p = 0.1, 0.5, 0.9, 0.99, which gives a sequence from (near) the chiral limit to (near) the limit of a single domain wall.For all plots, we have set ⟨B⟩/(eω 0 ) = 1, which exaggerates the effect of deviations from a constant magnetic field for illustrative purposes.Our expansion assumes these deviations to be small.

Figure 5 . 3 ,
Figure 5. Baryon density in units of nuclear saturation density n 0 for the physical pion mass m π = 140 MeV, hexagonal lattice in the transverse plane, a = 1/ √ 3, and elliptic moduli p = 0.3, 0.5, 0.7, 0.999, from top left to bottom right, which correspond to four points on the CSL instability line (µ [MeV], eB c2 [GeV 2 ]) ≃ (366, 0.83), (628, 0.28), (889, 0.13), (1016, 0.059).The lattices are computed away from the instability line with ϵ = 0.1 (first three plots) and ϵ = 0.3 (bottom right plot).Here and in all following numerical results we have used the value of the pion decay constant f π = 92 MeV.The vertical direction is distinguished by the magnetic field, which points (predominantly) into the z direction, see Fig.3.The variation of n B in the transverse direction is due to charged pions, whose vorticity contributes to the topological baryon number, in addition to a magnetic field contribution.

Figure 6 . 3 ,e B c 2 [
Figure 6.Left panel: Free energy difference (5.24) at the physical pion mass as a function of the geometric parameter a, normalized to the value at a = 1, for three different values of the elliptic modulus p.The vertical dashed line indicates a = 1/ √ 3, where all curves have their minimum.Right panel: Free energy difference for the hexagonal lattice a = 1/ √ 3 along the critical magnetic field, for three values of the pion mass.The dashed vertical lines indicate the smallest critical field of the CSL instability, where CSL approaches a single domain wall (p = 1).The solid vertical lines indicate the type-I/type-II transition point, where ∆f changes sign; for critical fields larger than this point, the pion crystal is preferred over CSL.

Figure 7 .
Figure 7.Chemical potential and magnetic field of the type-I/type-II transition point as a function of the pion mass.The dashed vertical line indicates the physical point, where (µ, eB c2 ) ≃ (907.7 MeV, 0.1234 GeV 2 ), as shown in the phase diagram in Fig. 1.