Shape stability of pasta phases: Lasagna case

The stability of periodically placed slabs occurring in neutron stars (lasagna phase) is examined by exact geometrical methods for the first time. It appears that the slabs are stable against any shape perturbation modes for the whole range of volume fraction occupied by the slab. The calculations are done in the framework of the liquid drop model and obtained results are universal --they do not depend on model parameters like surface tension or charge density. The results shows that the transition to other pasta shapes requires crossing the finite energy barrier.


Introduction
The very long structures appearing in neutron star matter, called pasta phases, are a commonly accepted phenomenon. Following the seminal work [1], the pasta phases have been studied in many different manners. Many approaches are based on the compressible liquid drop model (CLDM) where the system is described by two homogeneous phases separated by a sharp boundary with nonzero surface tension [2][3][4]6]. The competition between the Coulomb and surface energy leads to different shapes which are usually described in the Wigner-Seitz approximation, where the geometry of the phase is imposed at the beginning. The spherical and cylindrical cells are not the correct unit cells as they are not able to fill the whole space by periodic placement. Structures obeying periodicity, in the form of gyroid and diamond-like shapes, were introduced in [7,8]. However, they must be treated only as an approximation of a true shape, because they do not satisfy the necessary condition for the cell energy extremum. It was shown in [9] that the condition relates the mean curvature H(x) of the cluster surface and electrostatic potential where σ and Δρ are surface tension and charge difference between the cluster and its surroundings. The constant C depends on pressure difference between neutron and proton phases and mean value of potential C = P N − P P + Δρ Φ P . It means the surface curvature depends on the potential distribution and is not constant in general. In fact, the phases considered in [7] represent surfaces with H = 0 and as such they cannot be the true a e-mail: skubis@pk.edu.pl solution of eq. (1). In works [7,8] it was also shown that they have larger energy than the flat slabs. Nevertheless, the energy difference is not large, so it seems that such triply periodic structures are likely to occur. Recent analysis of pasta phases, including the periodicity, based on quantum or classical Molecular Dynamics [10][11][12][13][14] and Time-Dependent Hartree-Fock [15,16] or Hartree-Fock with twist-averaged boundary conditions [17] have shown the existence of triply periodic ones in the form of gyroid, diamond, twisted spaghetti and waffles. Topological analysis of such structures was presented in [18].
The general motivation of our work is to confirm the existence of such non-trivial periodic structures also in the framework of CLDM. Therefore, we propose to examine the possibility of transition from lasagna phase to other shapes by its shape deformation. The inspection of various perturbation modes could indicate into what kind of structures the lasagna is going to transform. Such consideration corresponds to the stability analysis of the lasagna phase. The answer to the stability question is not obvious. Though the contribution from surface energy for flat slab is always positive, the Coulomb energy is not and could destabilize the lasagna.
Another type of considerations, connected to the stability issue, were presented in the works [19][20][21][22] where the hydrodynamic approach was used to describe the modes with the density perturbation in pasta phases. The obtained mode frequencies were always real, which means that the lasagna phase is stable for this kind of collective modes. Our approach concerns a different kind of excitations -it considers only the shape-changing modes while the nucleon density is kept constant. Such analysis makes the overall stability discussion more complete.
In the work [9] an overall discussion of the rigorous treatment of periodically placed proton clusters of any shape was presented. The stable cluster surface should satisfy not only eq. (1), which represents the necessary condition for the extremum, but also the condition for the minimum coming from the inspection of second variation of the energy functional. The second order analysis, in a limited sense, was already carried for cylinders and balls in the [5], where a particular deformation mode was examined in the isolated Wigner-Seitz cell. Periodically placed cylinders and balls were considered in [6] and later in [23], but we need to be aware that these structures do not represent the proper minimum determined by eq. (1), their shape was assumed a priori. Up to now, the only known true periodic solution of eq. (1) is lasagna (the boundaries of slabs coincide with the equipotential surfaces) and thus the stability analysis of these structure may be carried out Partial stability analysis of lasagna phase was done by Pethick and Potekhin in the context of elastic properties of the phase in the work [23]. In order to determine the elasticity coefficient, they considered one particular deformation mode and moreover the expression for the energy was valid only in the limit of wavelength going to infinity. Such analysis corresponds to the stability consideration, however, being limited to the only one type of deformation. Here, we present general, unconstrained stability analysis for any kind of deformation with finite wavelength, which, to the authors' knowledge, has never been carried out.

Energy variation for single slab
Let us consider a proton cluster P in the shape of one slab placed in the center of unit cell with size a, b, c and volume V C = abc. The slab is perpendicular to the xaxis and occupies a w fraction of the total cell volume, 0 < w < 1. The charge density contrast between phases is Δρ = ρ + − ρ − , so the slab is positively charged with the density ρ + = (1 − w) Δρ and immersed in negatively charged neutron gas ρ − = −w Δρ (we assume the homogeneous electron background). For such charge distribution, the unperturbed potential Φ 0 , fig. 1, is a function of x only and takes the form where x 1,2 = ± aw 2 correspond to the positions of slab edges. Such configuration fulfill the necessary conditions for minimum, i.e. equations coming from vanishing variation of total energy in the first order. The sufficient condition for minimum is expressed by the positive value of total energy variation in the second order δ 2ε where the constraints for baryon number and charge conservation are imposed. Such energy variation for any deformation of the proton cluster is expressed by the integral over proton cluster surface ∂P [9] where B 2 = κ 2 1 + κ 2 2 = 0 is the sum of squared principal curvatures, ∂ n Φ 0 is the normal derivative of unperturbed potential, δ Φ is the first order perturbation of the potential caused by the surface deformation . One must remember that here we consider only the deformation which preserves the volume of the cluster, which means Deformations of this kind may be expressed in terms of Fourier series on the y, z-plane. The deformations on each slab face are independent, so we get two series for each face located at x 1 and x 2 . For further calculations it is convenient to introduce an expansion based on complex amplitudes α j mk , where m and k are the mode indices and the superscript j corresponds to the face number where we introduce the 3-dimensional discrete wave vector and x = (x, y, z). The m and k indices take integer values from −∞ to +∞ except the case when both of them equal zero, which is indicated by an apostrophe in the sum sign. The lack of (0, 0)-mode is consistent with the volume conservation condition, eq. (4). The inclusion of such compression modes would require taking into account the particle density perturbation. Their stability is controlled mainly by the volume compressibility coefficient The slab deformation in the unit cell (the grey rectangle). The first face is only with sine mode whereas the second face with cosine mode.
-quantity being dependent on the details of nuclear interactions. The compression modes may be carried out separately as the (0, 0)-mode is orthogonal to the shape changing modes. As we are interested only in the shape stability, we postpone the compression modes for future work.
Since the deformation function j (y, z) must be real it means that the complex amplitudes should fulfill the following relations: It is more natural to introduce the cosine mk,C cos( my b + kz c ) and sine modes mk,S sin( my b + kz c ) in the expansion eq. (5). Then the relation between complex and real amplitudes is In fig. 2 an example of the slab deformation with real amplitudes taking the values 1 10,S = 2 10,C = 0 and 1 10,C = 2 10,S = 0 is shown. Below we present and discuss the subsequent contribution to the total energy variation δ 2ε , eq. (3). In terms of complex amplitudes the surface energy contribution to the 2nd order variation takes a simple form, It is always positive, as B 2 = 0 for the slab and due to relations (7) we get α j mk α j −m,−k = |α j mk | 2 > 0. That is an important fact, because in many cases the surface energy acts as a destabilizer, like, for example, δ 2ε S < 0 in the case of Rayleigh-Plateau instability.
The contribution coming from the electrostatic interaction includes two terms. The first one, determined by the normal derivative of potential Φ 0 , is given by The normal derivatives at the slab faces are ∂ n Φ 0 ( which is always negative. The second term of the Coulomb interaction part is associated with perturbation of the potential δ Φ where δ Φ is the first order potential perturbation calculated thanks to the periodic Green function For the three-dimensional unit cell with sizes a, b, c the periodic Green function can be expressed as the sum over discrete modes numbered by three indices n, m, k The prime sign in the summation means that n, m, k cannot vanish simultaneously, similarly as in eq. (5). Such Green function fulfills the Poisson equation in the unit cell with periodic boundary conditions [24,25] (we use CGS units) and the absence of (0, 0, 0) mode in the expansion (15) is a consequence of the unit cell neutrality. By joining eqs. (5), (14), (15) we obtain the change in the energy with respect to the potential variation where ξ ij corresponds to the difference between the faces location and χ mk is the norm of the dimensionless wave vector The vector χ mk describes the manner in which the face is vibrating -the indices m and k determine the number of wavelengths being placed in the face sizes b and c. The function F (ξ, χ) comes from the summation over the n index. As the n-numbered modes do not depend on x, the summation over n can be done separately and defines a function F The above series can be evaluated [26] to a closed form .
The function is positive for any ξ and χ. In some special cases the function simplifies to Taking together all energy terms we get the total energy variation expressed in terms of mode amplitudes α mk As the δ 2ε is the quadratic form of the amplitudes α mk the competition between the terms in curly brackets decides about the stability of the face surface. There are three characteristic terms: one from the surface energy being positive and two from Coulomb interactions: the first is negative and the second is always positive. It seems that the stability consideration depends on the values of surface tension σ or charge contrast Δρ but it appears that these parameters may be removed from our analysis. One of the conditions for the minimum of the energy for unit cell is the virial theorem [9]. The theorem takes the form of relation between the surface and Coulomb energy of the cell. For the cell with high symmetry, it has the simple form from which we may get the relation between σ and Δρ Finally the total energy variation may be written as the quadratic form of α i with its coefficients A ij mk given by The obtained general expression, eq. (27) for the second order variation of the total energy allows us to determine whether the slab is stable with respect to any deformations preserving its volume. It is worth noting that A ij mk coefficients, which decide about stability, do not depend on the strength of interactions being determined by surface tension σ and charge contrast Δρ. In this way, we obtained an interesting result that the stability of pasta depends only on the geometry of phase and mode under consideration and not on the details of strong or electromagnetic interactions.
Before general discussion of the stability of a single slab we show how the above results work in the stability analysis for particular class modes.

An example of stability analysis
Let us consider the simplest surface perturbation consisting of combination of sine and cosine modes going along the y-axis for each face, which corresponds to m = ±1 and k = 0. Then, for the i-th face, the only non-vanishing complex amplitudes α i mk are where we introduced the real amplitudes i S and i C corresponding to the functions sin( 2πy b ) and cos( 2πy b ). The deformation i for the i-th face is then a function of y only It is convenient to introduce the vector built of deformation amplitudes Then the total energy variation for such deformation may be written down in the matrix form where the dimensionless matrixM iŝ and its elements are where w is the volume fraction and χ determines the wavelength of the mode in comparison to the cell size χ = 2πa/b. The inspection of eigenvalues and eigenvectors ofM allows for a complete stability analysis for the deformation we have chosen. The matrixM given by eq. (34) possesses the two-fold degenerated two eigenvalues. Further, we call theM -eigenvalues λ l , l = 1 . . . 4 as the stability function. For our concrete form ofM the eigenvalues and their eigenvectors are The stability functions λ i do not depend on the details of interactions σ, Δρ but only on the geometry of our system which is described by volume fraction w and mode wavelength χ. In fig. 3 the stability functions are plotted in the w, χ parameter space for their eigenvectors. These vectors represent two classes of modes. We may call them as snaky and hourglass-shaped modes. The snaky modes occur when the deformations on both faces are in phase ( 1 , 3 ) and hourglass modes occur when the deformations on faces are out of phase ( 2 , 4 ). As we may notice, for all values of volume fraction and mode wavelengths the stability functions are positive, which means that for both classes of modes the slab is stable. However the stability is not too strong. Careful inspection of λ 1 and λ 2 shows that those functions go to 0 when χ and w approach to some specific values and λ 2 | χ→∞ → 0 for w → 0 or 1.
It means that snaky modes become unstable in the limit of very long waves regardless of slab thickness, whereas the hourglass modes become unstable for very thin slab and very short waves. To sum up, we may say the slab becomes asymptotically unstable for very long mode or for very short mode when the slab becomes very thin in comparison to unit cell width. One should note that the case when w ≈ 0 or w ≈ 1 must be treated with caution because in reality the cluster surface has finite thickness and for very thin slab the validity of the liquid drop model could be questioned. Nevertheless, the first case of asymptotic instability, eq. (37), is worthy of careful inspection as it is connected to the macroscopic deformation and allows for determination of elastic properties of lasagna phase, which is shown in the next section.

Elastic properties
Nuclear pastas share their elastic properties with liquid crystals. If the wavelength of the snaky mode is very large in comparison to the cell size, it corresponds to the socalled splay deformation of the liquid crystal. It allows to determine elastic constant K 1 for the lasagna phase what was shown by Pethick and Potekhin in [23]. By definition, the constant K 1 relates the deformation energy to the transverse derivative of the deformation field when the mode wavelength becomes very large. In our notation the relation takes the form where the brackets . . . mean the average over the slab surface. Taking the definition of K 1 in the limit of the very long mode, we get (here denote the mode amplitude only). The deformation energy for the snaky mode is Applying the limit in eq. (39) we get which exactly corresponds to the result of [23] if the Δρ is replaced by the unperturbed Coulomb energy of the cell ε C,0 = 1 6 πa 2 Δρ 2 (1 − w) 2 w 2 . In comparison to [23] our approach represents an improvement. Pethick and Potekhin used a triple sum for the deformation energy, we mean eq. (8) in [23], which gives correct result only in the limit of the very long mode χ → 0. In fact, that series represents an asymptotic expansion in the powers of mode wavelengths and the leading term gives correct result. However the expansion is divergent for finite χ. The detailed discussion of this divergence is shown in the appendix. The approach, presented here, allows to avoid any divergences and is not limited to the case χ → 0.

General discussion of stability
In sect. 3 the stability analysis for the simplest slab deformation was carried out. The full stability analysis would require the inclusion of modes for all multiplicities m and k. Writing down the matrixM for all modes, it appears to take the block-diagonal form where in the mk-th position we get the same matrix as in eq. (34) with A and B elements taking the same form as in eq. (35) but with the replacement χ → χ km . Such a form ofM means that the modes with different multiplicity m, k and the cosine-like and sine-like modes do not couple. So, taking fixed m, k one may repeat the discussion from sect. 3 and finally conclude that the single slab is stable for any mode keeping in mind the asymptotic cases, described by eqs. (37), (38).

Multi-slab modes
In previous sections the cell with only one proton layer was considered. However, one may consider many slabs which are placed in a cell with periodic boundary conditions. Such system allows to test a larger class of perturbations and makes our analysis of stability more complete. Let us suppose we have N slabs in one cell. The expressions derived earlier for energy variation (10), (11), (17) take the same form, except the summation over the surfaces which now takes the range i, j = 1 . . . 2N , as we have 2N surfaces for N slabs. Moreover, the value of the normal derivative of the electric potential scales with the number of slabs N according to the rule

Conclusions
In this work, by use of analytical methods, we have shown that proton clusters having the form of slabs placed periodically in space are stable for all values of volume fraction occupied by the cluster. It is a compelling result. All works based on CLDM (see for example [1,4,6,7]) show that different shapes of pasta are preferred for different values of w. Here, we have shown that the lasagna phase is stable in the whole range of w. One must remember that our analysis means that the lasagna phase represents merely a local minimum. The transition to another geometry is not totally blocked, but requires finite size deformation in order to exceed the energy barrier. It may be interpreted as the fact that, at least for some range of volume fraction, the lasagna phase is metastable. That could be quite interesting for the dynamics of pasta appearance during the neutron star formation. Our analysis is based on small deformations so it cannot state at which range of w the lasagna represents global minimum of the cell energy. First, the global minimum statement requires the knowledge of exact solutions of eq. (1) for other kind of shape than flat slab. So far, we have not known such solutions in the CLDM approach. The seeking of them marks out the direction of further research of pasta by differential geometry methods.
We are also conscious that the approach, based on the CLDM, has its limitations and the inclusion of such effects like finite thickness of the cluster surface or the temperature fluctuations could change the final conclusion concerning lasagna phase stability.