Shape stability of pasta phases: Lasagna case

TL;DR

This paper demonstrates that neutron star pasta phases, specifically lasagna slabs, are geometrically stable against shape perturbations across all volume fractions, with stability results being universal and independent of model parameters.

Contribution

First exact geometrical analysis showing the universal stability of lasagna pasta phases in neutron stars against shape perturbations.

Findings

Lasagna slabs are stable against all shape perturbations.

Stability results are universal, independent of surface tension or charge density.

Transition to other pasta shapes requires crossing a finite energy barrier.

Abstract

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

Figures

Equations

$$2\sigma H(\mbox{\boldmath$x$})={\cal C}+\Delta\rho\;\Phi(\mbox{\boldmath$x$})~{},$$

$$\displaystyle\Phi_{0}(\mbox{\boldmath$r$})=$$

$$=$$

$$\begin{array}[]{ll}\delta^{2}\tilde{\varepsilon}=&\\ {\displaystyle\frac{1}{2V_{\cal C}}}\int_{\partial\cal P}\left(\sigma((\nabla\epsilon)^{2}-B^{2}\epsilon^{2})+\Delta\rho\;(\partial_{n}\Phi_{0}\;\epsilon^{2}+\delta_{\epsilon}\Phi\;\epsilon)\right)dS,&\end{array}$$

$${\int }_{\partial \mathcal{P}}ϵdS=0.$$

$$\epsilon^{j}(y,z)=\sideset{}{{}^{\prime}}{\sum}_{m,k=-\infty}^{\infty}\alpha_{mk}^{j}\exp(\mathrm{i}\mbox{\boldmath$K$}_{0mk}\cdot\mbox{\boldmath$x$})~{},$$

$$\mbox{\boldmath$K$}_{nmk}=\left(\frac{2\pi n}{a},\frac{2\pi m}{b},\frac{2\pi k}{c}\right)$$

$$({\alpha }_{mk}^j{)}^{\ast }={\alpha }_{-m,-k}^j.$$

$${ϵ}_{mk,C}^j$$

$${ϵ}_{mk,S}^j$$

$${\delta }^2{\tilde{\epsilon }}_S$$

$${\delta }^2{\tilde{\epsilon }}_{norm}$$

$${\delta }^2{\tilde{\epsilon }}_{norm}=-\pi \Delta {\rho }^{2}w(1-w)\sum _{j=1}^2\text{sideset}{}^{\prime }{\sum }_{mk}{\alpha }_{mk}^j{\alpha }_{-m,-k}^j$$

$${\delta }^2{\tilde{\epsilon }}_{\Phi }=\frac{\Delta \rho }{2V_{\mathcal{C}}}{\int }_{\partial \mathcal{P}}{\delta }_{ϵ}\Phi ϵdS,$$

$$\delta_{\epsilon}\Phi(\mbox{\boldmath$x$})=\Delta\rho\int_{{\partial\cal P}}G_{P}(\mbox{\boldmath$x$},\mbox{\boldmath$x$}^{\prime})\epsilon(\mbox{\boldmath$x$}^{\prime})\;dS^{\prime}~{}.$$

$$G_{P}(\mbox{\boldmath$x$},\mbox{\boldmath$x$}^{\prime})=\frac{4\pi}{abc}\sideset{}{{}^{\prime}}{\sum}_{n,m,k=-\infty}^{+\infty}\frac{\exp({\mathrm{i}\mbox{\boldmath$K$}_{nmk}\cdot(\mbox{\boldmath$x$}-\mbox{\boldmath$x$}^{\prime})})}{K_{nmk}^{2}}~{}.$$

$$\nabla^{2}G_{P}(\mbox{\boldmath$x$},\mbox{\boldmath$x$}^{\prime})=-4\pi\left(\delta(\mbox{\boldmath$x$}-\mbox{\boldmath$x$}^{\prime})-\frac{1}{abc}\right)$$

$${\delta }^2{\tilde{\epsilon }}_{\Phi }=2\pi \Delta {\rho }^2\sum _{i,j=1}^2\text{sideset}{}^{\prime }{\sum }_{m,k}{\alpha }_{mk}^i{\alpha }_{-m,-k}^{j}F({\xi }_{ij},{\chi }_{mk}),$$

$${\xi }_{ij}=\frac{x_i-x_j}{a}$$

$${\chi }_{mk}=a{K}_{0mk}=\sqrt{{\left(\frac{2\pi ma}{b}\right)}^2+{\left(\frac{2\pi ka}{c}\right)}^2}.$$

$$F(\xi ,\chi )\stackrel{\mathrm{df}}{=}\sum _{n=-\infty }^{\infty }\frac{{\mathrm{e}}^{\mathrm{i}2\pi n\xi }}{(2\pi n{)}^2+{\chi }^2}.$$

$$F(\xi ,\chi )=\frac{\cosh \left(\frac{\chi }{2}(1-2\left|\xi \right|)\right)}{2\chi \sinh (\frac{\chi }{2})}.$$

$$F(0,\chi )$$

$$F(\frac{1}{2},\chi )$$

$${\delta }^2\tilde{\epsilon }=\sum _{i,j=1}^2\text{sideset}{}^{\prime }{\sum }_{m,k}$$

$$+2\pi \Delta {\rho }^{2}F({\xi }_{ij},{\chi }_{mk})\}{\alpha }_{mk}^i{\alpha }_{-m,-k}^j.$$

$$E_S=2E_{Coul},$$

$$\sigma =\frac{1}{6}\pi a^3\Delta {\rho }^2(w-1{)}^2{w}^2.$$

$${\delta }^2\tilde{\epsilon }=\pi \Delta {\rho }^2\sum _{i,j=1}^2\text{sideset}{}^{\prime }{\sum }_{m,k}{\mathcal{A}}_{mk}^{ij}{\alpha }_{mk}^i{\alpha }_{-m,-k}^j$$

$${\mathcal{A}}_{mk}^{ij}=$$

$$+2F({\xi }_{ij},{\chi }_{mk})\}.$$

$${\alpha }_{10}^i$$

$${\alpha }_{-10}^i$$

Full text

††journal: Eur. Phys. J. A

1 Introduction

The very long structures appearing in neutron star matter, called pasta phases, are a commonly accepted phenomenon. Following the seminal work ravenhall83 , 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 non-zero surface tension Baym:1971ax ; hashimoto ; Oyamatsu:1993zz ; Williams:1985prf . 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 Nakazato:2009ed ; Nakazato:2010nf . 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 Kubis:2016fmw that the condition relates the mean curvature H(\mbox{\boldmathx}) of the cluster surface and electrostatic potential \Phi(\mbox{\boldmathx})

[TABLE]

where $\sigma$ and $\Delta\rho$ are surface tension and charge difference between the cluster and its surroundings.The constant $\cal C$ depends on pressure difference between neutron and proton phases and mean value of potential ${\cal C}=P^{\cal N}-P^{\cal P}+\Delta\rho\;\langle\Phi\rangle_{\cal P}$. It means the surface curvature depends on the potential distribution and is not constant in general. In fact, the phases considered in Nakazato:2009ed represent surfaces with $H\!=\!0$ and as such they cannot be the true solution of Eq. (1). In works Nakazato:2009ed ; Nakazato:2010nf 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 triple-periodic structures are likely to occur.

Recent analysis of pasta phases, including the periodicity, based on quantum or classical Molecular Dynamics Watanabe:2004tr ; Schneider:2013dwa ; Schneider:2014lia ; Horowitz:2015gda ; Alcain:2014fma and Time-Dependent Hartree-Fock Schuetrumpf:2012cj ; Schuetrumpf:2014aea or Hartree-Fock with twist-averaged boundary conditions Schuetrumpf:2015nza 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 Kycia:2017ibr .

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 DiGallo:2011cr ; Urban:2015dna ; Kobyakov:2018ftm ; Durel:2018cxs 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 Kubis:2016fmw 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 the 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 Iida:2001xy , where a particular deformation mode was examined in the isolated Wigner-Seitz cell. Periodically placed cylinders and balls were considered in Williams:1985prf and later in Pethick:1998qv , but we need to be aware that these structures do not represent the proper minimum determined by the Eq. (1), their shape was assumed apriori. 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 Pethick:1998qv . 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.

2 Energy variation for single slab

Let us consider a proton cluster ${\cal P}$ in the shape of one slab placed in the center of unit cell with size $a,b,c$ and volume $V_{\cal C}\!=\!abc$. The slab is perpendicular to the $x$-axis and occupies a $w$ fraction of the total cell volume, $0<w<1$. The charge density contrast between phases is $\Delta\rho=\rho_{+}-\rho_{-}$, so the slab is positively charged with the density $\rho_{+}=(1-w)\,\Delta\rho$ and immersed in negatively charged neutron gas $\rho_{-}=-w\,\Delta\rho$ (we assume the homogeneous electron background). For such charge distribution, the unperturbed potential $\Phi_{0}$, Fig. 1, is a function of $x$ only and takes the form

[TABLE]

where $x_{1,2}=\pm\frac{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 $\delta^{2}\tilde{\varepsilon}$ where the constraints for baryon number and charge conservation are imposed. Such energy variation for any deformation of the proton cluster $\epsilon$ is expressed by the integral over proton cluster surface ${\partial\cal P}$ Kubis:2016fmw

[TABLE]

where $B^{2}=\kappa_{1}^{2}+\kappa_{2}^{2}=0$ is the sum of squared principal curvatures, $\partial_{n}\Phi_{0}$ is the normal derivative of unperturbed potential, $\delta_{\epsilon}\Phi$ is the first order perturbation of the potential caused by the surface deformation $\epsilon$. One must remember that here we consider only the deformation which preserves the volume of the cluster, which means

[TABLE]

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 $\alpha^{j}_{mk}$, where $m$ and $k$ are the mode indices and the superscript $j$ corresponds to the face number

[TABLE]

where we introduce the 3-dimensional discrete wave vector

[TABLE]

and \mbox{\boldmathx}=(x,y,z). The $m$ and $k$ indices take integer values from $-\infty$ to $+\infty$ 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 $K_{V}=\frac{\partial P}{\partial\rho}$ - 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 $\epsilon^{j}(y,z)$ must be real it means that the complex amplitudes should fulfill the following relations

[TABLE]

It is more natural to introduce the cosine $\epsilon_{mk,C}\cos(\frac{my}{b}+\frac{kz}{c})$ and sine modes $\epsilon_{mk,S}\sin(\frac{my}{b}+\frac{kz}{c})$ in the expansion Eq. 5. Then the relation between complex and real amplitudes is

[TABLE]

In the Fig. 2 an example of the slab deformation with real amplitudes taking the values $\epsilon_{10,S}^{1}=\epsilon_{10,C}^{2}\neq 0$ and $\epsilon_{10,C}^{1}=\epsilon_{10,S}^{2}=0$ is shown.

Below we present and discuss the subsequent contribution to the total energy variation $\delta^{2}\tilde{\varepsilon}$, Eq. (3). In terms of complex amplitudes the surface energy contribution to the 2nd order variation takes simple form

[TABLE]

It is always positive, as $B^{2}=0$ for the slab and due to relations (7) we get $\alpha_{mk}^{j}\alpha_{-m,-k}^{j}=|\alpha_{mk}^{j}|^{2}>0$. That is an important fact, because in many cases the surface energy acts as a destabilizer, like, for example, $\delta^{2}\tilde{\varepsilon}_{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 $\Phi_{0}$, is given by

[TABLE]

The normal derivatives at the slab faces are $\partial_{n}\Phi_{0}(x_{1})=\partial_{n}\Phi_{0}(x_{1})=-2\pi(1-w)wa\Delta\rho$. So, finally one gets

[TABLE]

which is always negative. The second term of the Coulomb interaction part is associated with perturbation of the potential $\delta_{\epsilon}\Phi$

[TABLE]

where $\delta_{\epsilon}\Phi$ is the first order potential perturbation calculated thanks to the periodic Green function

[TABLE]

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$

[TABLE]

The prime sign in the summation means that $n,m,k$ cannot vanish simultaneously, similarly as in the Eq. (5). Such Green function fulfills the Poisson equation in the unit cell with periodic boundary conditions marshall ; tyagi (we use CGS units)

[TABLE]

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

[TABLE]

where $\xi_{ij}$ corresponds to the difference between the faces location

[TABLE]

and $\chi_{mk}$ is the norm of the dimensionless wave vector $\bm{\chi}_{mk}=(\frac{2\pi ma}{b},\frac{2\pi ka}{c})$

[TABLE]

The vector $\bm{\chi}_{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(\xi,\chi)$ 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$

[TABLE]

The above series can be evaluated gradshteyn to a closed form

[TABLE]

The function is positive for any $\xi$ and $\chi$. In some special cases the function simplifies to

[TABLE]

Taking together all energy terms we get the total energy variation expressed in terms of mode amplitudes $\alpha_{mk}$

[TABLE]

As the $\delta^{2}\tilde{\varepsilon}$ is the quadratic form of the amplitudes $\alpha_{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 $\sigma$ or charge contrast $\Delta\rho$ 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 Kubis:2016fmw . 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

[TABLE]

from which we may get the relation between $\sigma$ and $\Delta\rho$

[TABLE]

Finally the total energy variation may be written as the quadratic form of $\alpha^{i}_{km}$

[TABLE]

with its coefficients ${\cal A}_{mk}^{ij}$ given by

[TABLE]

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 ${\cal A}_{mk}^{ij}$ coefficients, which decide about stability, do not depend on the strength of interactions being determined by surface tension $\sigma$ and charge contrast $\Delta\rho$. 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.

3 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\!=\!\pm 1$ and $k\!=\!0$. Then, for the i-th face, the only non-vanishing complex amplitudes $\alpha^{i}_{mk}$ are:

[TABLE]

where we introduced the real amplitudes $\epsilon^{i}_{S}$ and $\epsilon^{i}_{C}$ corresponding to the functions $\sin(\frac{2\pi y}{b})$ and $\cos(\frac{2\pi y}{b})$. The deformation $\epsilon^{i}$ for the i-th face is then a function of $y$ only

[TABLE]

It is convenient to introduce the vector $\epsilon$ built of deformation amplitudes

[TABLE]

Then the total energy variation for such deformation may be written down in the matrix form

[TABLE]

where the dimensionless matrix $\hat{M}$ is

[TABLE]

and its elements are:

[TABLE]

where $w$ is the volume fraction and $\chi$ determines the wavelength of the mode in comparison to the cell size $\chi=2\pi a/b$. The inspection of eigenvalues and eigenvectors of $\hat{M}$ allows for a complete stability analysis for the deformation we have chosen. The matrix $\hat{M}$ given by Eq. (34) possesses the two-fold degenerated two eigenvalues. Further, we call the $\hat{M}$-eigenvalues $\lambda_{l}~{},~{}l=1\dots 4$ as the stability function. For our concrete form of $\hat{M}$ the eigenvalues and their eigenvectors are

[TABLE]

The stability functions $\lambda_{i}$ do not depend on the details of interactions $\sigma,{\Delta\rho}$ but only on the geometry of our system which is described by volume fraction $w$ and mode wavelength $\chi$. In Fig. 3 the stability functions are plotted in the $w,\chi$ 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 (\mbox{\boldmath\epsilon}_{1},~{}\mbox{\boldmath\epsilon}_{3}) and hourglass modes occur when the deformations on faces are out of phase (\mbox{\boldmath\epsilon}_{2},~{}\mbox{\boldmath\epsilon}_{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 $\lambda_{1}$ and $\lambda_{2}$ shows that those functions go to 0 when $\chi$ and $w$ approach to some specific values

[TABLE]

and

[TABLE]

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\approx 0~{}{\rm or}~{}w\approx 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.

4 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 so-called 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 Pethick:1998qv . By definition, the constant $K_{1}$, relates the deformation energy with the transverse derivative of the deformation field when the mode wavelength becomes very large. In our notation the relation takes the form

[TABLE]

where the brackets $\langle..\rangle$ mean the average over the slab surface. Taking the definition of $K_{1}$ in the limit of the very long mode, we get

[TABLE]

(here $\epsilon$ denote the mode amplitude only). The deformation energy for the snaky mode is

[TABLE]

Applying the limit in Eq.(39) we get

[TABLE]

which exactly corresponds to the result of Pethick:1998qv if the $\Delta\rho$ is replaced by the unperturbed Coulomb energy of the cell $\varepsilon_{C,0}=\frac{1}{6}\pi a^{2}\Delta\rho^{2}(1-w)^{2}w^{2}$. In comparison to Pethick:1998qv our approach represents an improvement.

Pethick and Potekhin used a triple sum for the deformation energy, we mean Eq.(8) in Pethick:1998qv , which gives correct result only in the limit of the very long mode $\chi\rightarrow 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 $\chi$. 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 $\chi\rightarrow 0$.

5 General discussion of stability

In the Section 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 matrix $\hat{M}$ for all modes it appears to take the block-diagonal form

[TABLE]

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. (3) but with the replacement $\chi\rightarrow\chi_{km}$. Such a form of $\hat{M}$ 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 the Section 3 and finally conclude that the single slab is stable for any mode keeping in mind the asymptotic cases, described by the Eqs. (37,38).

6 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

[TABLE]

As was discussed in section 42 it is enough to test the only one type of mode per surface. Then the stability matrix $\hat{M}$ for one given mode has dimension $2N\times 2N$. As an example we show results for the case of two slabs, $N=2$ and cosine-mode with $m=1,k=0$. Then we have four distinct eigenvalues of stability matrix $\hat{M}$ :

[TABLE]

where $D$ is

[TABLE]

All of these eigenvalues are positive. Their dependence on volume fraction $w$ and mode wavelength $\chi$ and the corresponding slab deformations are shown in the Fig.4. As one may see, the eigenmodes of the deformation of multi-slab system are always the combination of snaky and hourglass modes.

7 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 ravenhall83 ; Oyamatsu:1993zz ; Williams:1985prf ; Nakazato:2009ed ), 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.

Appendix

The Coulomb energy in the work Pethick:1998qv , Eq.(8), was finally expressed by the power series

[TABLE]

here we keep the notation used by Pethick and Potekhin: $2r_{c}$ is the cell size, $k$ is the mode wavenumber and $\xi$ is the dimensionless amplitude of perturbation. Coefficients $a_{j}(\xi)$ are determined by the double sum

[TABLE]

The series given by $\cal S$ is divergent. It is enough to take only one term from the double sum, Eq.(A.2) for $m=1$ and $n=j/2$ and then every coefficient $a_{j}(\xi)$ is bounded from below by the expression

[TABLE]

The $\frac{j}{2}$-th Bessel function for small amplitude ($\xi<1$), again, may be bounded by its Taylor expansion

[TABLE]

So, finally, all the coefficients of the power series $\cal S$ are underestimated by the new ones $\tilde{a}_{j}$

[TABLE]

where we removed factorial by use of the Stirling formula. The power series with coefficients $\tilde{a}_{j}$ is however divergent because for large $j$ the $\tilde{a}_{j}$ behaves like $j^{j}$ and the convergence radius is zero

[TABLE]

This divergence comes from the introduction of the sum over $j$ in the expression for Coulomb energy. The Coulomb energy, was initially expressed by the convergent double sum over $m,n$ indices, Eq.(7) of Pethick:1998qv , then the third summation over $j$ was introduced in the following way (screening was neglected $k_{TF}=0$)

[TABLE]

and then the order of summation sequence was changed $\sum_{n,m}\sum_{j}\rightarrow\sum_{j}\sum_{n,m}$. The change in the summation order is allowed only if all sub-series are convergent, however, the expansion given by Eq.(A.6) is convergent only if

[TABLE]

That means that the indices are not independent: the summation sequence is not interchangeable. The authors of Pethick:1998qv got a correct result because they took the double limit: $\xi\rightarrow 0$ (small amplitude of deformation) and $\frac{\pi}{kr_{c}}\rightarrow\infty$ (very long mode). Summarizing, we want to emphasize that such expression for the Coulomb energy, Eq.(8) in Pethick:1998qv , is valid only in the limit of the very long mode, $kr_{c}\rightarrow 0$ .