Chiral pasta: Mixed phases at the chiral phase transition

Interiors of neutron stars are ultradense and may contain a core of deconfined quark matter. Such a core connects to the outer layers smoothly or through a sharp microscopic interface or through an intermediate macroscopic layer of inhomogeneous mixed phases, which is globally neutral but locally contains electrically charged domains. Here I employ a nucleon-meson model under neutron star conditions that shows a first-order chiral phase transition at large densities. In the vicinity of this chiral transition I calculate the free energies of various mixed phases — different “ pasta structures ”— in the Wigner-Seitz approximation. Crucially, chirally broken nuclear matter and the approximately chirally symmetric phase (loosely interpreted as quark matter) are treated on the same footing. This allows me to compute the interface profiles of bubbles, rods, and slabs fully consistently, taking into account electromagnetic screening effects and without needing the surface tension as an input. I find that the full numerical results tend to disfavor mixed phases compared to a simple steplike approximation used frequently in the literature and that the predominantly favored pasta structure consists of slabs with a surface tension Σ ≃ 6 MeV = fm 2 .


I. INTRODUCTION
Neutron stars probe a large range of baryon densities, from subsaturation densities in the outer layers up to several times nuclear saturation density in the core. Quantum chromodynamics (QCD) predicts the existence of deconfined quark matter at ultrahigh densities, but it is unknown from first principles at which density the transition from nuclear matter to quark matter occurs. While theoretically extremely challenging, input to this question can be obtained from astrophysical observations. Quark matter has different thermodynamic and transport properties compared to nuclear matter, and thus a sizable quark matter core may have observable consequences for masses, radii, cooling behavior, etc., of single neutron stars as well as for gravitational wave signals from neutron star mergers.
In this paper, I am interested not primarily in a bulk phase of quark matter, but in the interface that separates nuclear matter from a potential quark matter core. This interface itself, depending on its size and structure, might influence observable properties of the neutron star. Since a continuity between quark and hadronic matter is conceivable [1][2][3][4][5], there may not be an interface at all but rather a smooth transition. The other possibility, predicted by various phenomenological models, is a first-order phase transition. In this case, there can be a sharp interface; i.e., the transition from nuclear matter to quark matter occurs via a domain wall of microscopic thickness. Depending on the extent of the discontinuity at the transition, it can induce qualitative changes in the mass-radius curve of the star [6] and possibly in the signals from neutron star mergers [7,8]. Or, the transition can be made less abrupt by a macroscopic region of a mixed phase, i.e., an inhomogeneous, possibly crystalline, phase, where quark and hadronic phases are separated spatially and occupy different volume fractions as the density varies [9][10][11]. Here, quark and hadronic phases each carry nonzero electric charge, such that their overall charge is zero in order to maintain global charge neutrality. The quark-hadron mixed phase can affect observable properties of the star through its transport properties [12] or through its rigidity, which may result in a sustained ellipticity of a rotating neutron star and thus in the emission of gravitational waves [13,14]. A similar mixed phase is expected to occur in the crust-core transition region of the star, where a crystalline coexistence of nuclei and a neutron (super)fluid is predicted to turn into homogeneous nuclear matter via different geometric structures as one moves to higher densities [15,16]. These structures, for example, bubbles, rods, or slabs of one phase immersed in the other, are often referred to as "nuclear pasta" due to their resemblance with gnocchi, spaghetti, and lasagne [17][18][19][20]. Here I will study pasta phases at the quark-hadron transition, more precisely at the chiral phase transition. 1 Various levels of approximations have been employed in the literature for mixed phases at the quark-hadron transition. Usually, nuclear and quark matter are treated with two different models [21][22][23][24][25][26][27][28][29][30]. In this approach, the profile of a quark-hadron interface cannot be calculated microscopically since this requires knowledge of an effective potential that connects the two phases. Instead, one may approximate the free energies of different mixed phase geometries by assuming spatially constant profiles on either side of the interface, with the interface itself being steplike [10], and treating the surface tension as an external parameter, either varying this parameter freely or employing constraints from microscopic calculations [31][32][33][34][35][36][37][38]. One can improve on this simple estimate by including the charge screening effect [39][40][41][42][43]. This renders the profiles nonconstant, but, if two different models are being used, it still requires an external choice for the surface tension.
I will employ a nucleon-meson model with a small explicit chiral symmetry breaking term [44][45][46][47], which shows a first-order transition to an approximately chirally symmetric phase. This model was used in Ref. [38] to compute the surface tension for isospin-symmetric nuclear matter at the chiral phase transition. I shall use the extension of the model to isospin-asymmetric matter [48] to construct mixed phases explicitly. This model has the obvious advantage of treating the phases on both sides of the transition on the same footing. In particular, location and properties of the chiral phase transition are determined by fixing all parameters to match known properties of nuclear matter at saturation, and interfaces between the two phases can be computed fully consistently, with the surface tension being a side result, not an input. The disadvantage is that the chirally restored phase is actually made of very light nucleons, not of quarks, and thus can only be a toy version of deconfined quark matter. Therefore, the quark-hadron pasta within this model is more accurately thought of as "chiral pasta." Obviously, it is desirable to have a unified description of quarks and hadrons which is more realistic on both sides of the transition. However, this is currently not possible from first principles and difficult to achieve even on the level of phenomenological models. (For a recent attempt in this direction using the gauge-gravity duality see Refs. [49][50][51][52].) I will work at zero temperature throughout the paper and employ the mean-field and no-sea approximations for the evaluation of the model. Furthermore, for the inhomogeneous phases I use the Thomas-Fermi approximation to integrate out the fermions and construct the mixed phase in the Wigner-Seitz approximation. In this approach, for a given geometry (bubbles, rods, slabs), the energetically preferred size of a unit cell is determined from the calculation of the interface profiles; see Fig. 1 for an illustration. This profile calculation is done by solving simultaneously and numerically the Euler-Lagrange equations for the mesonic condensates together with the Poisson equation for the electrostatic potential. A similar calculation was performed for the interface between vacuum and nuclear matter [53], but, to the best of my knowledge, the present paper is the first such calculation in the context of the quark-hadron transition.
The paper is organized as follows. In Sec. II, I introduce the model, explain the various approximations, derive the Euler-Lagrange equations, and comment on the model parameters and the specific choices I make for them. Section III contains an analysis of the homogeneous phases and serves to identify the regime in which mixed phases can potentially be found. The main results are presented in Sec. IV. After recapitulating the simple steplike approximation in Sec. IVA and explaining the numerical evaluation in Sec. IV B, I show and interpret the results for the chiral pasta phases in Sec. IV C. A summary and outlook are given in Sec. V. I use natural units c ¼ ℏ ¼ k B ¼ 1 and Heaviside-Lorentz units for the gauge fields, where the elementary electric charge is e ¼ ffiffiffiffiffiffiffiffi 4πα p ≃ 0.3 with the fine structure constant α. The convention for the metric tensor in Minkowski space is g μν ¼ ð1; −1; −1; −1Þ. Illustration of the geometries (bubbles, rods, slabs) of the "pasta phases" considered in this paper. The inner regions (shaded) are occupied by one phase, and the rest of the unit cell (volume bounded by the dashed lines minus the shaded region) by the other phase. Together with the complementary structures for bubbles and rods, there are five different geometries in total. The profiles of the condensates and the electrostatic potential across the interfaces (here shown as sharp surfaces for illustrative purpose) are computed numerically, and the free energy is minimized with respect to the size of the unit cell L. The unit cell is taken to have the same geometry as the structure of the pasta, thus approximating the actual Wigner-Seitz cell.

A. Nucleon-meson Lagrangian
The starting point is the Lagrangian [48] L ¼ψðiγ μ ∂ μ þ γ 0μ Þψ þ where the traces are taken over isospin space. The nucleon field includes neutron and proton spinors with corresponding chemical potentials μ n , μ p , such that is the chemical potential matrix in isospin space, with the baryon chemical potential μ B and the isospin chemical potential μ I . The mesonic fields are σ, ω μ , π ¼ π a τ a , and ρ μ ¼ ρ a μ τ a , where τ a (a ¼ 1, 2, 3) are the Pauli matrices, normalized such that Tr½τ a ; τ b ¼ 2δ ab , ½τ a ; τ b ¼ 2iϵ abc τ c . The vector meson mass will be set to m v ≃ 782 MeV. The kinetic terms of the vector mesons are defined through and F μν ¼ ∂ μ A ν − ∂ ν A μ is the electromagnetic field strength tensor with the electromagnetic gauge field A μ . The potential for the sigma and pion fields is taken to be and with the pion decay constant f π ≃ 93 MeV and the model parameters a 1 , a 2 , a 3 , a 4 , and ϵ. By requiring that in the vacuum hχi ¼ f π , the pion mass term fixes a 1 ¼ m 2 π with the pion mass m π ≃ 139 MeV, as well as the explicit chiral symmetry breaking term ϵ ¼ m 2 π f π . In addition to these parameters, the model also contains the coupling constants g σ , g ω , g ρ , which determine the nucleon-meson coupling via Yukawa interactions. The choice for the six remaining free parameters will be discussed in Sec. II D.

B. Free energy density
To facilitate the numerical evaluation of the pasta phases, I make use of several approximations. Employing the mean-field approximation, I neglect all meson fluctuations, such that the meson contributions arise only through the mesonic condensatesσ,ω, andρ. Here,σ ≡ hσi plays the role of the chiral condensate, giving rise to the nucleon mass. The omega condensate is the expectation value of the temporal component of the omega vector field,ω ≡ hω 0 i, and the rho condensate corresponds to the temporal component of the third isospin component of the rho vector field,ρ ≡ hρ 3 0 i. The condensatesω andρ act as effective contributions to the baryon and isospin chemical potentials. The resulting mean-field Lagrangian is equivalent to a Lagrangian of noninteracting fermions, with all effects of the interactions absorbed in the effective mass and effective chemical potentials. The ansatz for the condensates does not include the possibility of a chiral density wave, where the chiral condensate oscillates spatially between the σ and π 3 directions with fixed modulus χ [54]. The chiral density wave or related, more complicated, inhomogeneous structures are expected to appear in the vicinity of a first-order chiral phase transition [55], even without relaxing the local neutrality constraint to a global one, which is necessary for the mixed phases discussed here. I leave the study of the interplay between the chiral density wave and chiral pasta for the future. For simplicity, I will also drop the vacuum contribution; i.e., I work in the so-called no-sea approximation. Since the model is renormalizable, it is straightforward to include this contribution, which would introduce the renormalization scale as an additional parameter. However, since the model is of a phenomenological nature to begin with, there is no guarantee that the vacuum terms would make the results more realistic, although there are examples where it does so [56]. There are cases where the sea contribution makes an obvious qualitative difference, for instance, in the presence of a magnetic field [57], but such an effect is not expected here.
In order to describe inhomogeneous mixed phases, the condensatesσ,ω, andρ must be allowed to vary in space. (There is no time dependence since I am only interested in equilibrium configurations.) Within the Thomas-Fermi approximation, this spatial dependence is neglected when integrating over the fermionic fields to compute the partition function and thus the free energy density Ω of the system. As a consequence, the resulting free energy density has the same form as in the homogeneous case with all mesonic condensates allowed to be inhomogeneous and with the addition of the mesonic gradient terms. This approximation is expected to be valid for slowly varying condensates.
Within these approximations, the calculation of the free energy density follows the standard procedure, and at zero temperature one obtains Here, the Coulomb energy density is given by the electric field E ¼ −∇ϕ, where the scalar electric potential ϕ is related to the electron chemical potential μ e by ϕ ¼ μ e =e.
In addition to the terms coming from the nucleon-meson Lagrangian (1), I have added a leptonic contribution from noninteracting, approximately massless electrons with chemical potential μ e and noninteracting muons with chemical potential μ μ and mass m μ ≃ 106 MeV, where I have introduced the zero-temperature pressure of a noninteracting spin-1 2 fermion species with chemical potential μ and mass m, with the Fermi momentum For m ¼ 0 the pressure reduces to the simple form used for the electrons in Eq. (8). The nucleonic contribution in Eq. (7) is where UðσÞ is the potential (5) with all fields replaced by their expectation values, i.e., χ ¼σ, where are the effective nucleon chemical potentials, and where is the effective nucleon mass. For low temperatures, where the neutrino mean-free path is much larger than the size of the neutron star, equilibrium with respect to the weak processes p þ l → n þ ν l , n → p þ l þν l where l ¼ e, μ is either an electron or a muon, and e → μ þν μ þ ν e , μ → e þν e þ ν μ yields the following relations between the nucleonic and leptonic chemical potentials 2 : I will use these constraints to eliminate μ p and μ μ , such that the only remaining chemical potentials are μ n and μ e .

C. Euler-Lagrange equations
From the free energy density (7) one can now straightforwardly derive the Euler-Lagrange equations, The last equation is the Poisson equation, which can also be written as where ρ is the electric charge density, and q ≡ n p − n e − n μ ; ð17Þ such that positive ρ and q correspond to an excess of proton over lepton charges. The scalar density is where k F;n and k F;p are the neutron and proton Fermi momenta. The baryon and isospin densities are 2 Notice that these relations hold for the actual chemical potentials μ n and μ p , not for the auxiliary quantities μ Ã n and μ Ã p . The reason is that the nucleonic single-particle energies at the Fermi momentum are μ n and μ p , not μ Ã n and μ Ã p . n B ¼ n n þ n p , n I ¼ n n − n p with the neutron and proton densities and the electron and muon densities are The four coupled second-order differential equations (15) will later be solved numerically for the four functionsσðrÞ, ωðrÞ,ρðrÞ,μ e ðrÞ. In the homogeneous limit, Eqs. (15) become coupled algebraic equations. Then, Eqs. (15a)-(15c) are the stationarity equations of the free energy with respect to the variablesσ,ω,ρ, and Eq. (15d) accounts for electric charge neutrality.

D. Parameter values
It remains to fix the parameters g σ , g ω , g ρ , a 2 , a 3 , a 4 . The matching procedure is almost the same as in Ref. [38], where only isospin-symmetric matter was discussed. We first set g σ ¼ 10.097, which follows from the nucleon mass in the vacuum m N ¼ 939 MeV and the relation m N ¼ g σ f π . The remaining five parameters depend, in a coupled way, on saturation properties of symmetric nuclear matter which are not all known to high accuracy. I use the saturation density n 0 ¼ 0.153 fm −3 and the binding energy E B ¼ −16.3 MeV, both known to good accuracy. Furthermore, I need the incompressibility K ≃ ð200-300Þ MeV [58,59], the effective Dirac mass at saturation M 0 ≃ ð0.7-0.8Þm N [10, [60][61][62][63][64], and the asymmetry energy S ≃ ð30.2-33.7Þ MeV [65,66]. For definiteness, I set S ¼ 32 MeV, as in Ref. [48], and M 0 ¼ 0.75m N . One can easily express the coupling constant g ω in an analytical form in terms of the experimental values, MeV is the baryon chemical potential at the onset of symmetric nuclear matter, and where is the Fermi momentum of symmetric nuclear matter at saturation. An analogous expression can be derived for the coupling constant g ρ , This relation is derived in Appendix A. The remaining parameters a 2 , a 3 , a 4 are then determined from the stationarity equation with respect toσ, the condition that nuclear matter at the onset has zero pressure, and that the incompressibility K at saturation have a certain value, where the expression for K is also derived in Appendix A; see Eq. (A15) for the final result. These three conditions are coupled and have to be solved numerically. Due to the experimental uncertainty in K, and in order to look for qualitatively different scenarios, I shall first treat K as a free parameter. For the actual numerical calculation of the mixed phases I will then focus on the single parameter value K ¼ 252 MeV, for which, with the given values for

III. IDENTIFYING MIXED PHASE REGIONS
Before calculating the profiles of the condensates in various mixed phase geometries, one needs to identify the regions in the phase diagram where the mixed phases are candidates to replace the homogeneous solutions. At zero temperature, the only independent thermodynamic variable is the neutron chemical potential μ n ; i.e., for a given set of model parameters, this amounts to identifying a range in μ n around the first-order phase transition.
The model has three different phases: the baryonic vacuum (V), chirally broken nuclear matter (N), and the chirally restored phase, which I loosely interpret as quark matter (Q), as discussed in the Introduction. Since chiral symmetry is explicitly broken, N and Q are not distinguished by symmetry and in principle there can be a smooth transition between them. However, in the parameter range of interest it turns out that they are separated by a first-order transition, and hence it does make sense to label them as distinct phases. For isospin-symmetric nuclear matter, there is a first-order transition from V to N, the "baryon onset," to whose properties the parameters of the model are fitted, as discussed above. Then, at larger chemical potentials, there is a transition to Q. This transition is of first order unless the incompressibility assumes unrealistically large values, in which case this transition becomes a smooth crossover [38].
For isospin-asymmetric, electrically neutral, betaequilibrated matter the overall phase structure is similar, but there are some differences, as Fig. 2 shows. The black curves in this figure correspond to the phase transitions between the homogeneous phases. There is a small region of unrealistically small K where a direct transition from V to Q occurs. In this region, nuclear matter is only metastable and never the state of lowest free energy for any neutron chemical potential. (For symmetric nuclear matter, this regime is larger, extending into the region of realistic values of K [38].) Then, there is a very small regime, hardly visible in the plot, where the onset to neutral baryonic matter is first order. For most of the parameter range-including all realistic values of K-the baryon onset is second order, shown by the dashed line. The solid line just next, and almost parallel, to the dashed line indicates a (weak) first-order phase transition within the nuclear matter phase. I am not aware of any other model showing such a first-order transition, and possibly it is an artifact of the model. In any case, this transition will not be relevant for the main results of this paper. A more detailed view of the homogeneous phases, for a particular value of K, can be found in Figs. 3 and 4, where the free energy density, the effective nucleon mass, and the electron chemical potential are plotted. All these results are obtained by solving the Euler-Lagrange equations (15a)-(15c) without gradient terms together with the local neutrality constraint from Eq. (15d). Figures 2-4 also contain results for the VQ and NQ mixed phases in the absence of surface and Coulomb effects, where Eq. (15d) is turned into a global neutrality condition, as I explain now.
Mixed phases can appear in the vicinities of the firstorder phase transitions if one allows for spatial regions that are electrically charged, with the overall charge of the system being zero. In the simplest approximation, Coulomb energy and surface tension are neglected, and the mixed phase is characterized solely by the volume fraction of one of the phases χ ∈ ½0; 1. In particular, in this approximation, the specific spatial structure of the system is completely irrelevant. Adding Coulomb and surface energies increases the free energy of the mixed phase, and thus this simple approximation overestimates the stability of the mixed phases and gives an upper limit for the range in μ n in FIG. 2. Mixed phase regions in the plane of incompressibility at saturation K and neutron chemical potential μ n , without Coulomb effects and surface tension for the nuclear-quark (NQ, red shaded region) and vacuum-quark (VQ, blue shaded region) mixed phases. Thick black curves are first-order (solid) and secondorder (dashed) phase transitions between homogeneous phases. The baryon onset is second order almost everywhere, and there is a first-order transition within the nuclear matter phase just after the onset. The grey band indicates the experimental uncertainty for K, and the horizontal thin dashed line shows the value of K used in all following results. which the mixed phase is favored. The calculation is done as follows. Suppose phases 1 and 2 occupy volume fractions χ and 1 − χ, respectively. Then, Eqs. (15a)-(15c) for each phase (six equations in total) together with the condition of global charge neutrality, with the charge densities for the two phases defined via Eq. (17), and the equality of the free energies densities, are solved simultaneously for the eight scalar variables σ 1 ;ω 1 ;ρ 1 ;σ 2 ;ω 2 ;ρ 2 ; μ e ; χ for a given neutron chemical potential μ n . As a result, one finds a range in μ n where χ decreases continuously from χ ¼ 1 at the lower boundary to χ ¼ 0 at the upper boundary. In other words, as χ goes from 1 to 0, one moves from a pure system in phase 1 through a mixture of both phases up to the pure phase 2. The resulting free energy densities for the VQ and NQ mixed phases are shown in Fig. 3. The VQ phase is shown for all μ n where it exists, even though in a large part of this regime it is not the state of lowest free energy. Here and throughout the paper the VN mixed phase is ignored for simplicity. The situation becomes even more complicated if the mixed phase in the vicinity of the first-order transition within the N phase is constructed, made of nuclear matter with two different densities. This mixed phase is not included in the analysis either. In any case, the focus of this paper is on the NQ mixed phase, which is stable throughout the regime where it is a solution. The corresponding curves for the effective masses and the electron chemical potential in Fig. 4 are therefore shown only for the NQ phase.
To get an idea of the potential importance of mixed phases in this model, I have determined the range for the mixed phases for all values of the incompressibility K. Of course, there are experimental constraints for K, and one might argue that large parts of the parameter space are therefore not interesting. However, one needs to keep in mind that the model used here is of phenomenological nature. Therefore, it is useful to scan the parameter space and look for interesting qualitative features regarding mixed phases. While certain features might not be realized in the present model for realistic values of K, these features might be realized in other models, or in improved versions of this model, or in QCD. In this spirit, one might also add variations in M 0 and S to extend the result presented here, as it was done, varying K and M 0 , for isospin-symmetric nuclear matter in Ref. [38]. The result for fixed M 0 and S, but K unconstrained, is shown in the form of the red and blue shaded areas in Fig. 2. In contrast to Fig. 3, the VQ mixed phase is now indicated only in the regime where it has lower free energy than all homogeneous phases and the NQ mixed phase. One observes that the NQ mixed phase region in the vicinity of the chiral phase transition becomes smaller as K is increased, eventually vanishing for very large K, where the transition becomes a smooth crossover. Interestingly, the results show the possibility of a direct transition from the VQ mixed phase to the NQ mixed phase. Such a behavior is suggested for K ≃ ð30-240Þ MeV. If realized in QCD, this might, for instance, be relevant for quark stars, which are usually thought of as having an outer layer where strangelets are immersed in an electron gas; i.e., they have a crystalline VQ crust and a homogeneous Q body. Figure 2 suggests that stars of the form VQ-NQ-Q are conceivable. In other words, nuclear matter might find its way into a quark star through a mixed phase that exists in a certain layer of the star. Moreover, the transition from the VQ to the NQ mixed phase raises the question whether a VQN mixed phase exists in the transition region; i.e., a phase in which baryonic vacuum, nuclear matter, and quark matter coexist in spatially separated regions. These intriguing questions are beyond the scope of the present paper but deserve further studies, be it in the present model or in another approach.
Here I proceed with the parameter set of Figs. 3 and 4, which shows the "standard" form of the NQ mixed phasesurrounded by pure nuclear matter N at lower density and by pure quark matter Q at higher density-and will no longer be concerned with the VQ mixed phase and any complications that are connected with it.

IV. PASTA PHASES
Whether chiral pasta phases are preferred or not is determined by the total free energy per unit volume, which is given by the spatial integral over the free energy density (7) where partial integration and the Euler-Lagrange equations (15a)-(15c) have been used, where the surface terms are dropped due to ∇ω ¼ ∇ρ ¼ ∇σ ¼ 0 at the boundaries, and where the Coulomb energy per unit volume has been separated, with the electric field E ¼ ∇μ e =e. One can also compute F by directly integrating over Ω from Eq. (7), without rewriting the gradient terms with the help of partial integration. It is a good numerical check, however, to use both expressions because Eq. (26) makes explicit use of the Euler-Lagrange equations. This check can thus be used to confirm that these equations are solved to good accuracy. For the determination of the preferred pasta structure, F is the only relevant quantity and there is no need to extract the surface tension Σ. Nevertheless, it is useful to compute Σ because, first, it can be used as an input for simple approximations that do not rely on the knowledge of the exact spatial profiles and, second, it is interesting to compare the result to the surface tension of isospin-symmetric matter. The surface tension is defined by the difference in free energies between a domain wall configuration and the corresponding homogeneous configuration, where x 01 and x 02 are the boundaries of the onedimensional structure and ðΩ N þΩ l Þ 0 ≡ðΩ N þΩ l Þ x¼x 01 ¼ ðΩ N þΩ l Þ x¼x 02 is the free energy of either phase at the boundary. Note that the definition (28) does not include the Coulomb energy w C . I will work within the Wigner-Seitz approximation, which reduces the problem to a single unit cell instead of the full crystalline structure. Moreover, as already shown in Fig. 1, the unit cells are assumed to have the form of the pasta structure itself, i.e., spherical for bubbles, cylindrical for rods, and rectangular for slabs, rather than using the actual Wigner-Seitz cells, which are the unit cells of the specific lattice structure (for instance, cubic, face-centered cubic). As a consequence, this approximation does not distinguish between different lattice structures of the bubbles or rods. The advantage is of course that the calculation becomes effectively one-dimensional for all three geometries.

A. Steplike approximation
The Wigner-Seitz approximation (with conveniently shaped unit cells) can be further simplified if, additionally, the profiles of the condensates are assumed to be steplike. This approximation is well known and frequently used in the literature. Nevertheless, I will briefly recapitulate it because it serves as an introduction to the basic concepts of the chiral pasta phases, which is useful for the more complete calculation explained subsequently. Also, I will use the results of the steplike approximation later as a comparison to the full numerical result.
In the steplike approximation, the sketches in Fig. 1 can be taken literally: the interfaces that separate the two phases from each other are assumed to be sharp surfaces. The condensates and the electron chemical potential assume the values calculated in the previous section, as shown in Fig. 4. Therefore, for each μ n for which a mixed phase is conceivable according to Figs. 3 and 4, there are two electric charge densities, one for each constituent phase, which I will denote by ρ 1 and ρ 2 , and a volume fraction χ. This is the input from the microscopic calculation. One can now choose a certain geometry, and a simple exercise in electrostatics gives the Coulomb energy per unit volume. By assumption, the electrostatic potential, which is generated by the steplike charge distribution and which adds to the electron chemical potential, does not backreact on the values of the condensates. The calculation of the Coulomb energy is presented in Appendix B, which leads to the wellknown result [10] where L 0 is the width of the inner region of the Wigner-Seitz cell (shaded regions in Fig. 1), i.e., the radius of the bubble or the radius of the rod or (half) the width of the slab, which is related to the volume fraction by where L is the width of the Wigner-Seitz cell (see also Fig. 1) and d is the codimension of the structure, i.e., d ¼ 3 for bubbles, d ¼ 2 for rods, and d ¼ 1 for slabs. Moreover, The Coulomb energy is obviously an energy cost for creating spatial regions that are electrically charged. As Eq. (29) confirms, this cost increases as the Wigner-Seitz cell becomes larger at a fixed volume fraction χ. (Increasing the width of the Wigner-Seitz cell at fixed χ is equivalent to increasing L 0 at fixed χ.) A counteracting effect comes from the surface energy per unit volume, which decreases as the size of the Wigner-Seitz cell is increased because a larger Wigner-Seitz cell means fewer surfaces per volume where the condensates have to interpolate between the different phases, which costs energy. For a given χ, the sum of Coulomb and surface energies always has a minimum value at a certain nonzero and finite L 0 because the Coulomb energy goes like L 2 0 and the surface energy goes like 1=L 0 . This will be different in the full calculation, for instance, due to screening effects.
Minimizing ΔF with respect to L 0 at fixed χ and inserting the result back into ΔF gives In order to determine the critical values of χ at which the geometry of the mixed phase changes, one simply has to evaluate the function d 2 f d ðχÞ for the five possible phases that arise from d ¼ 1, 2, 3 and from replacing χ → 1 − χ in the cases d ¼ 2, 3. The latter accounts for the fact that the bubbles or rods are made of phase 1 immersed in phase 2 or vice versa. For slabs the situation is symmetric and no additional information is obtained by the replacement χ → 1 − χ. One finds that the geometry changes from bubbles to rods to slabs to complementary rods to complementary bubbles as the volume fraction is decreased from 1 to 0. With the help of Eq. (31) one finds that the four critical values of χ are given by Gð1 − χÞ ¼ Hð1 − χÞ ¼ HðχÞ ¼ This results in the critical values These values are completely independent of the details of the underlying microscopic model. Only the translation from χ into μ n and the actual value of the free energy cost (34) depends on the microscopic physics. We show the result of the steplike approximation in comparison to the full result in Fig. 8 for two fixed values of the surface tension. (As the full calculation will show, the surface tension itself is a function of μ n .)

B. Numerical evaluation
A more complete and reliable result is obtained by solving numerically the coupled differential equations (15) without further approximations. To get a first idea of the form of the solution it is useful to discuss the behavior in the center and at the edge of the unit cell. This behavior is found by linearizing Eqs. (15) about the boundary values. A detailed derivation, presented in Appendix C, shows that the lowest-order behavior is quadratic for all condensates and the electron chemical potential, ρðxÞ ≃ρ 0 þ ðm 2 vρ − g ρ n I Þ 0 Here, I have denoted the relevant coordinate collectively by x for all three geometries; i.e., x is a Cartesian coordinate for slabs, the cylindrical radial coordinate for rods, and the spherical radial coordinate for bubbles. The point at the center or the edge of the unit cell is denoted by x 0 , andσ 0 ;ω 0 ;ρ 0 ; μ e0 are the boundary values. The subscripts 0 at the coefficients of the quadratic terms indicate that the expressions are evaluated at the boundary. Moreover, k ¼ 2 at the outer boundary of the unit cell for all geometries, while in the center k ¼ 2d. Importantly, the boundary values themselves are not known a priori, and they cannot be computed by a local analysis. They rather have to be computed dynamically by solving the coupled differential equations over the entire domain of the unit cell. The coefficients of the quadratic terms are nothing but the right-hand sides of the Euler-Lagrange equations (15). This shows that the profiles would be flat (the coefficient of the quadratic term would be zero) if the boundary values fulfilled the stationarity and neutrality equations. The more a certain unit cell deviates from these properties at its boundaries, the less flat the profiles become. Besides these general insights into the solution, Eqs. (37) are also useful as a check for the numerical solution, at least in the vicinity of the boundaries. I perform the full numerical calculation along the following steps: (1) Fix a neutron chemical potential μ n for which the mixed phase is preferred in the absence of Coulomb energy and surface tension, i.e., some μ n in Fig. 3 in the range of the (red) NQ curve. At this given μ n , two charged phases with the same μ e have the same pressure and form a globally neutral system with a volume fraction χ. For notational compactness, let v ≡ ðσ;ω;ρ; μ e Þ and let v 1 and v 2 be the values of v in these two charged phases. (2) Pick a geometry, i.e., bubbles, rods, or slabs. For bubbles and rods, decide which of the two phases occupies the interior of the cell. This yields five cases, and in all cases the Wigner-Seitz approximation ensures that the problem is effectively onedimensional. (3) Choose a domain, x ∈ ½x 01 ; x 02 , such that L ¼ jx 01 − x 02 j is the width of the unit cell. In the case of bubbles and rods x is a radial coordinate and thus x 01 ¼ 0. Construct four initial guess functions vðxÞ, which are as close as possible to the (yet unknown) solutions. This can be done by interpolating between the values of the homogeneous phases, vðx 01 Þ ¼ v 1 , vðx 02 Þ ¼ v 2 , for instance, with the tanh function in the case of the condensates and with the constant function μ e ðxÞ ¼ μ e in the case of the electron chemical potential. (4) Solve the four Euler-Lagrange equations (15) for v numerically with the successive over-relaxation method. 3 Since the initial guess functions contain a nontrivial electric charge distribution, the Poisson equation induces a nontrivial behavior μ e ðxÞ. In particular, μ e ðxÞ will differ at the boundaries from the (constant) starting value. This, in turn, also changes the condensates at the boundaries. Therefore, it is important to keep all boundary values dynamical in the relaxation procedure. (5) Enforce electric neutrality in each iteration step of the relaxation procedure: By integrating the charge density over the unit cell, the total charge is a functional of the profiles computed at a given iteration step, in particular of the electron chemical potential Q ¼ Q½μ e ðxÞ. Electric neutrality is enforced by adding a constant contribution φ that solves the algebraic equation Q½μ e ðxÞ þ φ ¼ 0. This defines the corrected electron chemical potential μ e ðxÞ þ φ used for the next iteration step. After sufficiently many iterations, the result is a stable interface for a given geometry and given L, where the boundary behavior of vðxÞ is given by Eqs. (37). (6) Repeat the procedure for different L. For instance, if L was chosen to be smaller than the (yet unknown) favorable L, successively increase x 02 by a small amount (while keeping x 01 fixed). It is useful to extend the solution from the previous, smaller, domain to the larger domain in order to use it as an initial guess, for instance, by adding constant pieces to the profiles. Increasing the domain leads to nontrivial changes, including small changes in the boundary values, and therefore the relaxation procedure has to be repeated for every L. (7) For each L, compute the free energy per unit volume F via Eq. (26) and determine the (local) minimum of F as a function of L if such a minimum exists. The final result, i.e., the preferred size of the unit cell together with the profile functions, only depends on μ n and the chosen geometry. Quantities such as L 0 , χ, and Σ, needed to construct the steplike approximation, play no role in the calculation, everything is determined dynamically.

C. Results and discussion
Using the slab geometry as an example, I will first explain some details of the results and afterwards present the main conclusions, bringing together the results of all three geometries.

Unit cells and profiles for slabs
In Fig. 5, results for three different neutron chemical potentials μ n are shown as a function of the width of the unit cell L. In the upper panel the free energy, calculated from Eq. (26), is plotted. One of the curves has a very pronounced minimum, as one would expect from the steplike approximation. However, moving to smaller and larger chemical potentials, i.e., toward more imbalanced volume fractions, the minimum becomes more shallow. For the (red and black) curves shown here the minimum barely The green line is the surface tension at the minimum for all μ n where a minimum exists. Here and in all following results, K ¼ 252 MeV. 3 The same method was used in the calculation of the surface tension of the VQ, VN, NQ transitions for symmetric nuclear matter [38] and for flux tube profiles of multicomponent superconductors [67,68]. I am working with a grid of about 5000 points on the domain, corresponding to a lattice spacing of about 10 −3 fm, and with up to about 10 6 iterations steps. These numbers were obtained by checking convergence for selected cases to reach an accuracy such that the results of all plots in the paper would be basically indistinguishable if more grid points or iterations were used. In the slab geometry convergence is much quicker, and a much coarser lattice and fewer iterations lead to the same results. exists and indeed the minimum ceases to exist for slightly smaller (red line) or larger (black line) chemical potentials. At least for the curve μ n ¼ 1016 MeV, it is obvious that the mixed phase with the slab geometry is metastable; i.e., the minimum is only local, not global. The system seems to be able to reduce its free energy by forming larger and larger unit cells. The corresponding profiles show that the volume partition between the two phases becomes more and more extreme; i.e., the tendency to form unit cells of large L is nothing but the tendency to create uniform matter. Since the numerics become challenging for very large L, it is difficult to continue the curve much further and show that it indeed asymptotes to the free energy of the uniform phase. But, it is easy to check that for the two cases μ n ¼ 1016 MeV and 1029 MeV the free energy of the uniform state (N in the case of the smaller μ n , Q for the larger one) is indeed smaller than the free energy of the local minimum, as shown explicitly in Fig. 8 below.
The lower panel of Fig. 5 shows the surface tension, computed from Eq. (28). First of all, not surprisingly, one sees that the surface tension is not a constant. For instance, it decreases as the unit cell gets larger for fixed μ n . The figure also shows the values of the surface tension at the energetic minimum for all μ n for which the slab solution is at least a local minimum. For these values ones reads off Σ ≃ ð5.2-6.2Þ MeV=fm 2 . These values are somewhat smaller than but very similar to the surface tension computed in the same model with isospin-symmetric matter from a domain wall configuration with semi-infinite phases at the chiral phase transition [38]. 4 Figure 6 shows the corresponding profiles across one wall of the slab for the same three neutron chemical potentials as in Fig. 5. In each case, the profiles are shown for the L at the local minimum. The coordinates are chosen such that x ¼ 0 corresponds to the center of the slab in all  4 With the chosen parameters M 0 ¼0.75m N and K ¼ 252 MeV, symmetric nuclear matter is only metastable in the given model, such that in Ref. [38] only the surface tension of the transition between the vacuum and the chirally restored phase was calculated at this parameter point; see, for instance, Fig. 6 in this reference. The same figure shows that for slightly larger values of K or M 0 , where symmetric nuclear matter is stable, one finds Σ ≃ 8 MeV=fm 2 , not far from the result obtained for asymmetric matter here. three cases, while the right end of the scale corresponds to the size of the largest of the three unit cells. For comparison, the figure also shows the constant values obtained by ignoring surface tension and Coulomb effects, i.e., the values that give the results for the NQ phase in Figs. 3 and 4. The plot confirms that the boundary values of the profiles deviate-sometimes significantly-from these constant values. The step that connects the constant values is chosen to reproduce the volume fractions from the NQ phases of Figs. 3 and 4, using the preferred unit cell size from the full calculation. 5 One can see from the curves for the largest and the smallest μ n that the full calculation tends to create a more pronounced imbalance between the two phases than predicted from the steplike curves.
The profiles for the effective mass M and the baryon number n B illustrate that the wall interpolates between the chirally broken (large M), low-density phase and the chirally restored (small M), high-density phase. One characteristic feature of the slab is the accumulation of isospin density at the interface, as the lower left panel shows. The profiles of the electric charge density q demonstrate the screening effect, i.e., the accumulation of negative and positive charge carriers close to the interface. This effect reduces the Coulomb energy cost and thus works in favor of large unit cells. One might wonder why the charge accumulation is not an even more pronounced effect. But of course one has to keep in mind that reducing the Coulomb energy is not the only effect at work. A more extreme charge accumulation is, for instance, prevented by the tendency of the condensates to remain as close as possible to the local minimum they would assume in the absence of electrostatic energy costs, in order to reduce the potential energy of the configuration.
The charge density profiles also show a cusp close to the interface at negative q. This cusp arises because there is a regime where there are no protons in the chirally broken phase: for all x up to the cusp for μ n ¼ 1023 MeV and 1029 MeV (blue and black curves) and for x ≃ 4 fm up to the cusp for μ n ¼ 1016 MeV (red curve). Thus, in this regime, the system decides to create baryon number in the chirally broken phase only from neutrons, with the lepton gas providing negative charge. This is possible because the chirally restored phase is positively charged. At this point, it is useful to recall that although in some sense the chirally restored phase can be interpreted as quark matter, there are important differences to realistic quark matter. First of all, this phase can, at best, be identified with two-flavor quark matter because strangeness is not included in the model. A priori, there is nothing wrong with this restriction because it is plausible that in the transition region under realistic neutron star conditions the strange quarks are sufficiently heavy to only play a very minor role [5]. However, in the context of locally charged mixed phases it is important to note that the chirally restored phase in the present model has different charge carriers than realistic two-flavor quark matter, which is made of up and down quarks with charges 2=3e and −1=3e. Most notably, in the present model there is no negative charge carrier that also carries baryon number, i.e., no analogue of the down quark. All negative charge is carried by the leptons. This has to be kept in mind in the interpretation of the results. In fact, this shortcoming of the model provides a strong motivation to include strangeness in future studies. The reason is that, even if hyperons might not play any role in the chirally broken phase, there will be negative baryonic charge carriers in the system, which, in the chirally restored phase, may provide a more realistic description of actual quark matter.

Comparison of all geometries
Figures 7 and 8 combine the results for all geometric structures considered in this paper. For the bubbles and rods one has to distinguish whether nuclear matter forms bubbles or rods immersed in quark matter ("N bubbles," "N rods") or vice versa ("Q bubbles," "Q rods"). Figure 7 shows the sizes L of the unit cell that minimize the free energy for a given neutron chemical potential μ n ; i.e., this is the result of the analysis illustrated by the upper panel of Fig. 5, now applied to all five geometries. In the case of the slabs, the three chemical potentials used for Figs. 5 and 6 are marked by diamonds. The curve for the slabs ends on both sides at a point, marked by a dot, beyond which the local minimum of the mixed phase configuration ceases to FIG. 7. Energetically preferred widths L of the unit cells as a function of the neutron chemical potential μ n for bubbles (green lines), rods (blue lines), and slabs (red line). Spherical dots mark the end of the line beyond which there is no stable mixed phase solution; at the other end of the curves for bubbles and rods the numerics start becoming too difficult. The (black) diamonds indicate the three points used for Figs. 5 and 6. exist. In the case of bubbles and rods the same disappearance of the local minimum occurs on one end of the curves. This is the end where the preferred uniform phase is the one that fills the bubbles or rods. In other words, these end points indicate the critical μ n 's above which Q bubbles and Q rods are no longer a minimum and the critical μ n 's below which N bubbles and N rods are no longer a minimum. The other end of the curves is in the regime where by reducing the size of the bubbles or rods to zero one arrives at the energetically preferred uniform phase. On this side, beyond a certain point the numerical relaxation algorithm does not yield any nontrivial profiles, the profiles rather relax to the homogeneous solution. This may well be a purely numerical problem; I cannot exclude that nontrivial profiles continue to exist beyond that point. Therefore, in Figs. 7 and 8 I have not marked these ends of the lines by a dot, suggesting that the lines potentially continue. Nevertheless, the range where I did find nontrivial solutions is sufficient in the sense that at the point where I have to stop one of the uniform phases is already favored, i.e., this numerical problem occurs in the regime where the pasta phases are metastable anyway. Figure 8 shows the main result of the paper by comparing the free energies of the uniform phases with the ones of the mixed phases from the full numerical calculation and from the steplike approximation. All free energies are plotted relative to the free energy of the mixed phase without any surface and Coulomb effects. This free energy difference is denoted by ΔF; i.e., ΔF ¼ 0 corresponds to the (red) NQ curve in Fig. 3. The steplike approximation shows the sequence of phases Q bubbles → Q rods → slabs → N rods → N bubbles, separated by black dots at the points given by the critical volume fractions (36). This sequence is generic for this approximation and does not depend on the details of the system. As discussed above, the steplike approximation requires a value for the surface tension as an input. For a rough comparison with the full result I have chosen the two constant values Σ ¼ 5.2 MeV=fm 2 and Σ ¼ 6.2 MeV=fm 2 , motivated by the minimal and maximal values obtained from the numerical calculation; see lower panel of Fig. 5. A density dependent surface tension between these two boundaries would then yield a curve within the band given by the two curves of constant surface tension. Considering only the most preferred configuration for each μ n , the steplike approximation with this choice of the surface tension predicts a sequence of phase transitions from uniform nuclear matter at low densities to Q rods, then to slabs, and then to uniform quark matter at high densities. I have checked that surface tensions Σ ≳ 13 MeV=fm 2 are needed to lift the entire dashed curve above the free energies of the uniform phases, i.e., to make the mixed phases disfavored for all μ n .
The solid lines show the full result, with the lower panel giving a more detailed view of the chiral pasta regime. The results indicate that the mixed phases are less favorable than predicted by the steplike approximation with the same surface tension. They exist in a range of neutron chemical potentials of Δμ n ≃ 2 MeV, compared to Δμ n ≃ 15 MeV for the steplike approximation and Δμ n ≃ 66 MeV if Coulomb and surface effects are entirely neglected. Interestingly, the preferred structures seem to be identical to the prediction of the steplike approximation: also in the full result most of the mixed phase regime consists of slabs, possibly with a small region of Q rods at the low-density end. All transitions between the different homogeneous and inhomogeneous structures are of first order, and in a more complete treatment more complicated geometries, or mixtures of different pasta structures, might occur in the vicinities of these first-order transitions. FIG. 8. Upper panel: Free energies per unit volume of homogeneous phases (black lines), chiral pasta phases from the full numerical calculation (solid red, green, and blue curves) and from the steplike approximation using the two constant surface tensions Σ ¼ 5.2 MeV=fm 2 and 6.2 MeV=fm 2 (thin dashed curves). All free energies are relative to the mixed phase without any Coulomb effects and surface tension, i.e., relative to the (red) NQ curve in Fig. 3. Lower panel: Magnification of the region where chiral pasta is stable or metastable, showing that slabs (and possibly Q rods in a tiny regime) are the preferred pasta configuration.
One might wonder whether it is a coincidence that the chiral pasta is just barely preferred and that predominantly slabs, not bubbles or rods, seem to appear. Is it possible that pasta is never preferred, i.e., that it is metastable for all densities? One way to check this would be to repeat the calculation for different model parameters. After all, there is some freedom in choosing the parameters due to the uncertainty in the saturation properties M 0 , S, and K. Due to the relatively tedious numerical calculation it is not easy to present a systematic survey of the parameter space, and I have not repeated the full-fledged calculation for different parameters. I did, however, repeat the calculation of the slabs for a larger value of the incompressibility at saturation K (still within the experimentally predicted range) while keeping M 0 and S fixed and found the free energy to be larger for any μ n than that of the homogeneous phases. This suggests that there is no pasta phase at all for large K. This is perhaps expected since it was already demonstrated in Fig. 2 that the potential range for a mixed phase, even in the absence of surface tension and Coulomb energy, becomes smaller as K is increased. Therefore, in this regime, the energy costs from surface and Coulomb effects can easily destroy the mixed phases. Consequently, the quantitative and even some of the qualitative details of Fig. 8 should not be taken too literally. Even within the chosen model they depend on the exact values of the model parameters.

V. SUMMARY AND OUTLOOK
I have computed pasta structures and their free energies consistently within a single model at the chiral phase transition. To this end, I have employed a phenomenological model where nucleons and mesons are the fundamental degrees of freedom, to which noninteracting electrons and muons are added. For model parameters that reproduce saturation properties of nuclear matter and with the constraints of electric charge neutrality and beta equilibrium, this model shows a first-order chiral phase transition at zero temperature and in the absence of mixed phases. I have first identified the region in the vicinity of this phase transition where globally, but not locally, neutral mixed phases are possible without taking into account surface and Coulomb effects. These effects have then been taken into account by calculating the profiles of the meson condensates and the electrostatic potential in a consistent way, solving the Euler-Lagrange equations for the condensates coupled with the Poisson equation for the electrostatic potential. Doing this for three different geometries-bubbles, rods, and slabs-and determining the preferred sizes of the unit cells dynamically for each case, this calculation yields the free energies of the different pasta structures. I have compared these free energies with each other, with the free energy of the homogeneous phases, and with the free energy obtained from a simple steplike approximation of the profiles, which is often used in the literature. As a result, I have found that for the chosen model parameters chiral pasta is favored in a vicinity of the first-order phase transition which is only about 2 MeV wide, measured in the neutron chemical potential, and that the predominant structure that appears consists of slabs. As a by-product, I have extracted the density dependent surface tension, computed from domain walls in the slab geometry, at the energetically preferred unit cell sizes. I have found values Σ≃ ð5.2-6.2Þ MeV=fm 2 , which are comparable to, although slightly smaller than the surface tension of isospinsymmetric nuclear matter at the chiral phase transition in the same model. All these results depend on the model parameters, which are not uniquely fixed due to the experimental uncertainty of some of the saturation properties of nuclear matter. Exploring this parameter space and the fate of the mixed phases systematically is left for future work, but I have pointed out that larger values of the incompressibility at saturation lead to even smaller and eventually vanishing regimes for the pasta phases.
I have used various approximations, which can be improved in future studies. For instance, one may go beyond the mean-field approximation, which becomes particularly relevant for nonzero temperatures, or improve on the Thomas-Fermi and Wigner-Seitz approximations that I have used for the inhomogeneous phases. Also, it should be emphasized that the high-density phase in the present model is only a very crude approximation to a realworld quark matter phase. In particular, there are no baryonic charge carriers with negative electric charge, which is relevant for the structure of the mixed phases. Therefore, one interesting extension would be to include strangeness, which introduces negative charge carriers in the form of hyperons and their chirally restored counterparts. Another extension would be to allow for an inhomogeneous chiral condensate in the form of a chiral density wave and study the interplay of the mixed phases with this inhomogeneous structure. Furthermore, the present model suggests that a vacuum-quark mixed phase might be adjacent in chemical potential to a quark-hadron mixed phase if the chiral phase transition is sufficiently close to the baryon onset. As a consequence, a layer containing nuclear matter within an inhomogeneous phase in quark stars or even the existence of a three-component mixed phase are conceivable. I have not addressed these intriguing possibilities in detail here, and it would be interesting to study them in the future. Finally it would, of course, be desirable to implement the results of this paper into a calculation of the structure of a neutron star or into a simulation of a neutron star merger, which potentially probes the fate of chiral pasta at nonzero temperatures.

APPENDIX A: ASYMMETRY ENERGY AND INCOMPRESSIBILITY AT SATURATION
The asymmetry energy is defined as [10] S ¼ 1 2 where ϵ is the energy density, where the derivatives are taken at fixed n B , and where the thermodynamic relation has been employed. The isospin chemical potential is μ I ¼ðμ n −μ p Þ=2, and thus with Eqs. (12), (15c), and (19) one computes Consequently, For the evaluation of this expression at the symmetric point, where μ Ã n ¼ μ Ã p , the derivative of M with respect to n I is not needed. With n I ¼ 0, evaluating the result at saturation and inserting it into the definition (A1) yields where This is in accordance with Eq. (4.218) in Ref. [10] (after g ρ → g ρ =2 due to the different definition of g ρ ) and gives Eq. (23) in the main text. The derivation of the incompressibility is similar. Since the incompressibility is needed for symmetric matter at saturation and no isospin derivatives are required, one can perform the entire calculation in the symmetric scenario.
The incompressibility is defined as again using the thermodynamic relation In analogy to Eq. (A3) one has where the omega condensate has been eliminated with the help of Eq. (15b). Consequently, and thus, with the definition (A7), Now the derivative of M does not drop out, and one has to compute it with the help of Eq. (15a). First note that at zero temperature Taking the derivative with respect to M on both sides of Eq. (15a) yields With Eq. (A12) one finds Inserting this into Eq. (A13) and the result into Eq. (A11) yields the incompressibility where the momentum integral can be done analytically, Z

APPENDIX B: COULOMB ENERGY IN STEPLIKE APPROXIMATION
In order to compute the Coulomb energy density (27) in the steplike approximation one first needs to compute the electric field. Due to the symmetry of the problem in all three cases (bubbles, rods, and slabs in unit cells of corresponding shapes), this is best done with the method of Gaussian surfaces. The integral form of Gauss's law is (in Heaviside-Lorentz units) where the integral is taken over a closed surface, which contains the total charge Q. With the charge densities ρ 1 (inner phase, shaded in Fig. 1) and ρ 2 (outer phase) and the volume fraction χ from Eq. (30), charge neutrality of the total cell can be written as

Bubbles
For the spherical geometry one has EðrÞ ¼ EðrÞe r , where r is the radial coordinate, and thus with a spherical Gaussian surface of radius r Eq. (B1) yields where Q is the charge contained in a sphere of radius r. The calculation of the electric field thus reduces to the calculation of that charge and one easily obtains Piecewise integration over the radial coordinate then yields the Coulomb energy (27) per unit volume, where the first result is general for any charge densities ρ 1 , ρ 2 , and in the second step the neutrality condition (B2) has been used.

Rods
In this case, using cylindrical coordinates whose radial component I also denote by r, the electric field has the form EðrÞ ¼ EðrÞe r . Using a cylindrical Gaussian surface with length L z and radius r, Gauss's law gives Calculating the enclosed electric charge Q yields the electric field EðrÞ ¼ 1 2 ( ρ 1 r for 0 < r < L 0 Therefore, the Coulomb energy per volume is again using charge neutrality in the second step.

Slabs
In this case, let the x direction be perpendicular to the slab, i.e., the slab is extended in the y and z directions, and let the point x ¼ 0 be such that the y − z plane is in the middle of the slab. Then one can focus on x > 0, and the electric field is EðrÞ ¼ EðxÞe x and Eð0Þ ¼ 0 for symmetry reasons. As a Gaussian surface one can, for instance, take a cylinder extended in the x direction with one end at x ¼ 0, such that the surface integral in Gauss's law only receives a contribution from the opposite end, where A is the circular area of the cylinder. The charge enclosed by the cylinder is easily computed, and one obtains This yields for the Coulomb energy one finds where ζ 0 0 η 00 0 − η 0 0 ζ 00 has been used, which can be checked explicitly for all three cases with the help of Eqs. (C12). One now has to remember that Eq. (C10) denotes only one of four components, and thus, returning to vector notation, the solution of Eq. (C7) close to x 0 is Multiplying both sides from the left with U and thus undoing the rotation that was necessary for the diagonalization, one obtains the simple result Note that only the lowest-order result allows for such a simple expression. The full result (for the linearized problem) is obtained by inserting C and D from Eq. (C14) into δṽ from Eq. (C11), and then undoing the rotation by multiplication with U. Since each of the four δṽ i depends on a different eigenvalue λ i , this leads to a linear combination with coefficients which are impossible to write down in a compact way (they can, of course, be evaluated numerically without problems). Only because the lowestorder result for δṽ i does not depend on the eigenvalue λ i it yields such a simple form for the unrotated functions δv i .
It remains to discuss the case x 0 ¼ 0. Obviously, in the case of slabs, there is no qualitative difference between the two ends of the domain, i.e., between the center of the slab and the boundary of the unit cell. Therefore, both ends are covered by the result (C19). In the center of the rod, where x 0 ¼ 0, the function ζðxÞ from Eq. (C12b) diverges, and thus regularity requires D ¼ 0. The boundary conditions (C13) are fulfilled with C ¼f 0 =λ, and thus (again omitting the index i) This can easily be expanded and again one is left with a result proportional tof 0 , such that after reinstating the vector index and undoing the rotation one obtains rod center : δvðxÞ ≃ x 2 4 f 0 : Finally, for the center of the bubble, regularity and the boundary conditions (C13) yield C ¼ −D ¼f 0 =ð2λ 3=2 Þ, such that In analogy to the center of the rod, one thus obtains bubble center : δvðxÞ ≃ x 2 6 f 0 : The results (C21) and (C23) differ from the result (C19) in the numerical prefactor. It is easy to see that this difference comes from the different form of the Laplacian (C2) and that the number in the numerator is simply 2d.
In summary, Eqs. (C19), (C21), and (C23) show the behavior of the profiles at the center and the edge of the unit cell for all geometries considered in this paper and are written in the form (37) in the main text.