Liquid-phase epitaxy of neutron star crusts and white dwarf cores

Near-equilibrium bottom-up crystallization of fully-ionized neutron star crusts or white dwarf cores is considered. We argue that this process is similar to liquid-phase epitaxial (i.e. preserving order of previous layers) crystal growth or crystal pulling from melt in Earth laboratories whereby lateral positions of newly crystallizing ions are anchored by already solidified layers. Their vertical positions are set by charge neutrality. Consequently, interplane spacing of a growing crystal either gradually increases, tracing $n_\mathrm{e}$ decrease, as the crystallization front moves away from the stellar center, or decreases, tracing decrease of $\langle Z \rangle$, when the crystallization front crosses a boundary between layers of different compositions. This results in a formation of stretched Coulomb crystals, in contrast to the standard assumption of cubic crystal formation, which is based on energetics arguments but does not take into account growth kinetics. Overstretched crystals break, which limits the vertical sizes of growing crystallites. We study breaking shear strain and effective shear modulus of stretched matter and discuss possibility of macrocrystallite formation. The latter has interesting astrophysical implications, for instance, appearance of weak crustal layers, whose strength may increase by a few orders of magnitude upon breaking and refreezing at a late-time event. We also analyze interaction of adjacent Coulomb crystals, having different ion compositions, and estimate the strength of such interfaces.


INTRODUCTION
Neutron stars (NS) and white dwarfs (WD) are magnificent astrophysical objects renowned for unparalleled breadth and significance of observational manifestations as well as for absolutely extreme physical conditions in their interior (e.g.Fontaine & Brassard 2008;Winget & Kepler 2008;Althaus et al. 2010;Kaspi 2010;Kaspi & Kramer 2015;Mereghetti, Pons & Melatos 2015;Özel & Freire 2016;Kaspi & Beloborodov 2017;Córsico et al. 2019;Saumon, Blouin & Tremblay 2022, and references therein).These stars, with masses in excess of 2 M⊙ for NS and up to about 1.4 M⊙ for WD, with thermonuclear reactions turned off, are kept from gravitational collapse by pressure of degenerate fermions, neutrons and electrons, respectively (e.g.Haensel, Potekhin & Yakovlev 2007).Their central densities are expected to be ∼ 10 15 g/cc for NS and up to ∼ 10 10 g/cc for WD.At densities ρ below the nuclear saturation density, ρ0 = 2.8 × 10 14 g/cc, i.e. in the outer layers of NS called crust and in the entire WD interior, structure of matter remotely resembles that of ⋆ E-mail:baiko@astro.ioffe.ruEarth metals.There are fully-ionized ions, i.e. atomic nuclei, which, in deeper layers of NS, become neutron-rich, and strongly degenerate nearly uniform electron gas.The electron gas is ultrarelativistic at ρ ≫ 10 6 g/cc.In NS at ρ > ρ d = 4.3 × 10 11 g/cc, there is also Fermi-liquid of neutrons dripped from the nuclei.In the very outer layers at ρ 10 4 g/cc (depending on ion sort), conditions of full ionization and electron degeneracy gradually cease to exist.
Ion composition of such dense matter is a complicated topic (e.g.Beznogov, Potekhin & Yakovlev 2021;Shchechilin, Gusakov & Chugunov 2023, and references therein).In general, the composition varies strongly with density.For WD, one typically considers less massive stars made of helium, stars of intermediate mass which contain a C/O mixture with a pronounced transition from deeper O-enriched regions to C-enriched matter in the outer core, and finally, the most massive stars made of an O/Ne mixture with a possibility of Mg, Si etc. Also there are smaller admixtures of various other elements and isotopes.In NS, though lighter elements may be present at relatively low densities, one typically considers heavier nuclei, for instance, iron (which represents matter ground state at the lowest densities) and many other.The composition in this case is calculated by minimizing thermodynamic potential of matter or by following a nuclear reaction network, taking place in a freshly accreted fuel.The end result of such calculations is a sequence of shells, each corresponding to a specific range of mass densities with composition dominated by a particular atomic nucleus.Such a structure is sometimes referred to as the "onion" structure of the NS crust.Predicted isotopes oftentimes have different charge-to-mass ratios and tend to gravitationally separate making the "onion" structure even more pronounced.
As an NS or a WD cools, their internal temperature drops below local crystallization temperature, which results in a gradual freezing of matter from the deeper layers up.Formed crystals, which are called Coulomb crystals to emphasize long-range nature and simplicity of forces, binding them together, are typically assumed to have (the most tightly-bound) body-centred cubic (bcc) lattice.In what follows, we shall neglect electron screening, i.e. slight deviations of electron density from uniformity.This is a well justified approximation in compact degenerate stars at not too low densities, i.e. down to densities, where full-ionization approximation starts breaking down (e.g.Haensel et al. 2007).Hence in our approach, electrons will be treated as an ideal constant density background.Besides that, electrons are responsible for thermal and electric conductivities of matter, which are believed to be very high (e.g.Potekhin, Pons & Page 2015) and thus provide for an exclusively uniform thermal environment.
In this paper, we take a fresh look at the bottomup crystallization process in dense matter of compact stars and critically reappraise the bcc lattice assumption, breaking strain, and shear modulus for these objects.Also, using methods of lattice dynamics, we study an idealized model of an interface between layers with different compositions, e.g.shells of the "onion" structure in NS crust.In particular, we analyze interaction between different layers and try to answer the question as to the effect of such interfaces on the overall crustal strength.A related problem of spherical Coulomb crystal energy, which is relevant for dusty plasma crystallization and ion crystallization in traps, has been investigated e.g. by Hasse & Avilov (1991, and references therein).
In section 2.1, we deduce a formula for the interaction energy of a charge with a Coulomb half-crystal, consisting of ordered planar layers of identical ions immersed into chargeneutralizing electron background.In section 2.2, this formula is manipulated to produce an alternative expression for the Madelung energy, which reproduces previously known result for the bcc lattice.In section 2.3, we consider interaction of two Coulomb crystals composed of ions of two different sorts, propose a simplified model of the "onion" structure boundary, and analyze its breaking properties.In section 3.1, which is central for the whole work, we establish the possibility of epitaxial crystal growth in dense matter of degenerate stars and predict crystal stretching effect, associated with it.In section 3.2, breaking properties of stretched matter are studied quantitatively.In particular, we consider NS crust or WD core composed of crystallites with random orientations of crystallographic planes with respect to the stretch direction (polycrystalline model).By contrast, in section 3.2.1,we address the possibility of large-scale crystallite formation and possible relations between their properties and astrophysical phenomena (macrocrystallite model).In section 3.3, shear modulus of stretched matter is evaluated.Finally, in section 3.4, we list a few other effects, which are expected to accompany elongation and contraction of Coulomb crystals in NS and WD interior.

Interaction of a charge with a crystal
We would like to begin by studying the interaction energy of a charge Q with a plane, containing a 2D lattice of ions with charges Ze (e is the positron charge), spanned by basis lattice vectors a1 and a2, and a slab of neutralizing electrons of constant density ne, enveloping symmetrically the ion lattice.If Q is outside of the electron slab, the latter can be shrunk to the same plane with surface density σ = neη, where η is the slab width.
Let us pick an ion of the lattice and define its position as the origin.Let the position vector of the charge Q be d = d ⊥ +d .In this case, index ⊥ indicates vectors perpendicular to the plane whereas index marks in-plane vectors.Then the interaction energy U1 is given by where R = R ⊥ +R is the difference between an ion position vector and the charge and an analogous representation can be written for the second term in equation (1).Using the fact that ν,µ e iq(νa where G are 2D reciprocal lattice vectors for direct vectors spanned by a1 and a2, we deduce In this case, r = −d .The ρ-integral can be evaluated explicitly with the result Suppose now that there is not one ion plane but half a space filled with such planes separated by distance η.Ion positions in the κ-th plane (κ = 0, 1, 2, . ..) are obtained from ion positions in the original plane by adding an out of plane vector κa3 (a3 = a 3⊥ + a 3 , |a 3⊥ | = η).We note in passing that any simple 3D lattice can be specified in this way.
Interaction energy of the charge Q, located outside of the half-space, with each of these planes is given by the same equation ( 6) with R ⊥ = κa 3⊥ − d ⊥ and r = κa 3 − d .Moreover, one can sum over κ and find interaction energy U of the charge with the whole half-space as

Madelung energy of a crystal
As a first application of this formula, we shall reproduce the Madelung energy of a 3D ion lattice based on it.We suppose, that Q = Ze, d = −a3, and, moreover, that the charge Q belongs to a plane of identical charges, adjacent to the filled half-space and positioned and oriented with respect to the original plane and the half-space as it would be in an ideal 3D crystal.Let this plane be also immersed into its own electron slab of the same density ne and thickness η.Then, it is necessary to find the interaction energy of Q with other charges in its plane.Using Ewald transformation, we can write where A > 0 is an arbitrary real number, R = νa1 + µa2, R = ν,µ , and G is the same as in the equation ( 5).
Madelung energy is the crystal electrostatic energy per ion.Only half of U0 thus contributes to the Madelung energy.This is pretty obvious for the 1/R term, but the situation with the ion-background term in the first line of equation (8) is slightly more subtle.For Madelung energy, it should be replaced by two other terms, an ion-background contribution equal to and a background-background term equal to (10) where r → ∞, and the area S = Z/neη.Such a replacement results in a correction of U bg = πZe2 neη 2 /6.
For the U contribution, equation ( 7), factor 1/2 is not needed because U describes interaction with a half-space, but the second half-space produces equal contribution.It is easy to verify, that, for this term, there are no corrections associated with treatment of the background.Thus where ζ is the Madelung constant, ai = (4πni/3) −1/3 is the ion sphere radius, and ni = ne/Z is the ion 3D density.

Interaction of two crystals
As another application of equation ( 7), we shall consider a simple model of an interface in the "onion" structure predicted to form in NS crust.Accordingly, we suppose that there are two semi-infinite one-component crystals composed of ions, one with charges Z1e, the other with charges Z2e, Z2 > Z1, touching along some plane.This plane can belong to different crystallographic families in the two crystals.Without loss of generality, we shall assume that the Z2phase is "below" the Z1-phase.The Z2-phase can be identified with the half-space discussed in section 2.1, whereas charges Z1e can be identified with Q.The Z1-lattice is then spanned by basis vectors b1 and b2, belonging to the same plane as the vectors a1 and a2, and by an out of plane vector b3.If we pick a Z1-ion in the bottom plane of the Z1-structure and denote its position as d (the origin is still at one of the Z2-ions of the top plane of the Z2-half-space), positions of all the other Z1-ions in this plane are given by d + nb1 + mb2.Ion positions in the k-th plane are obtained by adding vector kb3.In this case, n, m = 0, ±1, ±2, . .., k = 0, 1, 2, . ...
It is crucial for our argument, that the phases have essentially the same 1 electron density ne, which, in the astrophysical context, is predominantly responsible for pressure and hydrostatic equilibrium.Then, if Z1-lattice is also bcc, its cube edge b l = (Z1/Z2) 1/3 a l .The interplane spacing (or the electron slab thickness) is h = |b 3⊥ |, and the distance from the top plane of the Z2-crystal to the bottom plane of the Z1-crystal is |d ⊥ | = 0.5(η + h).If both cubes are aligned with the boundary plane, i.e. this plane is a {100} plane for both crystals (more on this notation in section 3.2), then The interaction energy of each Z1-ion with the Z2lattice is given by equation ( 7), in which one has to replace |d ⊥ | by |d ⊥ | + k|b 3⊥ | and d by d + nb1 + mb2 + kb 3 .The sum over Z1-ions yields the interaction energy 2 between the two crystals.The main question then becomes, whether the interaction energy depends on d , i.e. on the Z1-lattice lateral 3 position with respect to the Z2-lattice.
Summation of the factor exp [−iG(d + nb1 + mb2 + kb 3 )] in equation ( 7) over n and m restricts the sum over G in equation ( 7) only to those {G ′ } ⊂ {G}, which are reciprocal vectors simultaneously for direct vectors generated by a1 and a2 and by b1 and b2.For there to be any such reciprocal vectors, vectors b1 and b2 must be linear combinations of a1 and a2 with rational coefficients tij: bi = j tij aj, i, j = 1, 2. Suppose, for brevity, that vectors a1 and a2 are those given in the end of section 2.2.Then the distance between Z1-ions separated by the vector b1 is a l t 2 11 + t 2 12 .On the other hand, the distance between any two ions in a bcc lattice with the cube edge b l is b l l/4 with l integer.Equating and squaring, we observe that for there to be any {G ′ } ⊂ {G}, (b l /a l ) 2 = (Z1/Z2) 2/3 must be rational, which is satisfied for very rare Z1/Z2.For other charge ratios, there are no appropriate G ′ , and the interaction energy between Z1and Z2-half-spaces is zero.Thus, the system is insensitive to their mutual orientation and lateral position (vertical separation is set by ne), so that they are free to move or rotate with respect to each other without any static friction (dynamic friction will be non-zero e.g.due to electron viscosity).We shall refer to such lattices as incommensurate.
In the same fashion, one can easily study a situation, in which one of the crystals (or both) is stressed in such a way that it is no longer cubic, but is commensurate with the lateral lattice spacing of the other crystal.Consequently, the set of G ′ will not be empty, and the interaction energy will not be zero.Moreover, summation over k can be performed then in a closed form with the result: where Uint is the total interaction energy of the crystals per each Z1-ion of the bottom plane.This expression does depend on d .The model of crystal interaction developed in this section is however too simplistic for two reasons.First is the fact, that, in the bottom plane of the Z1-lattice as well as in the top plane of the Z2-lattice, there is a complex set of potential wells, produced by both half-spaces (cf.section 3.1).Hence, it would be more realistic to assume the presence of at least one intermediate plane with a mixed composition between the two half-lattices.Using equation ( 7) twice, one could generate a set of potential wells in this layer due to both half-spaces, which would have an aperiodic quasicrystal arrangement and different depths.It would be possible to maintain overall charge neutrality by counting these wells, specifying the relative fractions of Z1-and Z2-ions, and adjusting the intermediate electron slab thickness accordingly.
Our preliminary study indicates that the static friction in this model depends on ion placement in the wells, but, in any case, is orders of magnitude lower than that for a movement of two halves of a perfect crystal with respect to each other.
Second and more important is the perception that the mutual motion along the interface would be impossible, because the interface would not be strictly planar (Fig. 1).The interface deviation from planarity by just 10 atomic planes is likely enough to ensure that it is not weakened in comparison with the weaker of the two touching crystals.

EPITAXIAL GROWTH OF COULOMB CRYSTALS AND THEIR ELASTIC PROPERTIES
3.1 Epitaxial freezing and crystal stretching in neutron star crusts and white dwarf cores We shall now turn to the problem of solidification in a compact star.Consider gradual freezing of a one-component crystal (100% of ions have charge Ze) and suppose that the freezing direction is vertical, while the surfaces of simultaneous freezing are horizontal.Previous, already frozen layers create a set of potential maxima and minima characteristic of a particular crystallographic plane.By contrast, liquid, on average, is uniform, neutral, and approximately the same near any crystal plane.
Examples of the potential relief for various crystallographic planes are shown in Figs.2a-d under the assumption that the frozen part of the crystal has a well-defined top plane and that crystal ions are located precisely at the nodes of the bcc lattice.Crystallographic notation is explained in detail in section 3.2 and Table 1.The contours in Figs.2a-d are calculated using equation ( 7) and plotted on the plane of the 2D vector d (in units of a l /2).Potential minima are at or near the centres, whereas potential maxima are at the corners and/or at the upper edges of the plots.The difference between maximum and minimum potential (in units of UC) is shown by the solid red curve in Fig. 2e as a function of the respective 3D lattice interplane spacing |a 3⊥ | = η in logarithmic scale.
The periodic potential in Figs.2a-d is a monotonous function of the distance |d ⊥ | from the last frozen plane (e.g.dashed curves in Fig. 2f).However, if one takes into account the potential of the electron slab associated with the newly forming ion layer (this contribution is independent of the lateral coordinate), the combined potential acquires a minimum as a function of the vertical position very close to the middle of the slab (cf.solid curves in Fig. 2f).More precisely, ion deviation from the mid-slab position by a length δ results in an addition of U slab (δ) to the ion potential energy, where For that reason, in Figs.2a-d, the height above the last frozen plane is taken equal to the respective interplane distance (or the electron slab thickness, U slab at δ = η/2 is shown by the dashed blue curve in Fig. 2e.For comparison, let us note that kinetic energy of a typical ion at freezing is ∼ (3/2)UC/175 (e.g.Haensel et al. 2007), which is much less than the typical energy scale of potentials shown in Fig. 2. Thus, lateral ion positions in a newly freezing layer are determined by the ions of the previous layers in analogy with liquid-phase epitaxial crystal growth or solidification from melt, the processes well-known in solid-state physics and semiconductor industry (e.g.Burton, Cabrera & Frank 1951;Chalmers 1954;Billig 1955;Jackson & Chalmers 1956;Pashley 1956;Jackson, Uhlmann & Hunt 1967;Stringfellow 1982;Bauser & Strunk 1984;Kuphal 1991;Scheel 2000, and references therein).Even though crystallization kinetics is a very complicated process (e.g.Dash 1959), we believe, that to intrinsically less hindered nature of 2D nucleation as opposed to 3D nucleation; due to lack of any ion-orientation dependence (typical of covalent bonds); due to long-range almost pure Coulomb forces, extending over a major portion of the elementary cell and capturing unbound ions (cf.especially Figs.2a-c); due to a strong pull on ions to settle at a correct height above the already crystallized surface (related to charge neutrality); and due to extremely slow, with plenty of time for anneal, near-equilibrium nature of freezing in compact stars, it should be even easier to accomplish epitaxial growth of Coulomb crystals in dense matter than, say, silicon in Earth laboratories.
For growth perpendicular to crystallographic surfaces with low Miller indices (Figs.2a-c), the lateral positions of potential minima produced by already frozen half-space are exactly the same as the lateral positions of the ions of the new layer in the infinite 3D lattice (crosses in Figs.2a-d).However, crystals can grow perpendicular to planes of lower symmetry (cf.Elbaum & Chalmers 1955;Weymouth & Soepono 1962;Bauser & Strunk 1984) and, in such cases, it is not necessarily so (e.g.Fig. 2d).It is likely, that the absence of the upper crystallized half-space then results in a built-in deformation near the interface nonuniform in the vertical direction (a deformation of this kind was discussed e.g. by Cabrera 1959).Horizontal deviations (the same for all ions) from the 3D lattice positions decrease with depth increase.As new layers are added to the top, the deformation moves upwards preserving self-similarity.Note, that the red curve in Fig. 2e displays the difference between maxima and minima as calculated for a perfect half-lattice (blue contours) and does not take this deformation, likely affecting several frozen planes at the top, into account.
As we have just seen, lateral positions of ions, being added to a growing crystal, are fixed by previous layers.Vertical positions of the ions, on the other hand, are determined by charge neutrality.Thus, in a freezing star, interplane distances gradually increase, tracing ne decrease, associated with pressure decrease as one moves away from the stellar center.The typical length-scale of ne variation is of the order of the pressure scale height, hP = P/∇P (more on these scales in section 3.2.1).This results in a formation of stretched (more precisely, elongated) ion crystals in the terminology of Baiko & Kozhberov (2017, hereafter Paper I).For self-consistency, we note that, qualitatively, the general picture of Fig. 2 is preserved for stretched matter.
Infinite 3D elongated Coulomb crystals develop unstable phonon modes, if the elongation exceeds a critical value (Paper I).Thus, overstretched crystal layers lose stability and get destructed as soon as their volume properties begin dominating.Since these layers have already cooled below melting, they can then freeze anew, somewhat above the destruction edge, into a stretch-free cubic structure.Above this new cubic seed, the process of freezing with gradual stretching will repeat itself.Below, the cubic solid will form an interface with stretched (but not overstretched) layers from the previous step, which have survived the destruction.
It is difficult to predict the specific structure of the interface.Presumably, the cubic and stretched crystals will be incommensurate, which means zero interaction energy.However, as explained in the end of section 2.3, their mutual motion along the interface will be likely impeded by interface roughness.The shear strength of the interface will be about the same as that of the stretched layers below the interface (cf.section 3.2).
A similar process takes place when the freezing front crosses a shell boundary of the "onion" structure, from Z2to Z1-enriched region with Z1 < Z2.Let us suppose that the composition changes from 100% Z2 to 100% Z1, and the ratio Z2/Z1 is such that a crystal-like structure with a regular arrangement of lattice nodes can exist for any intermediate composition.It is natural to assume that the mixed layer is thin, so that, to zeroth order, there is no change of ne, which is determined by the pressure profile of the star.Then again, already frozen layers anchor the lateral positions of ions of the newly forming layers.There is, however, a drop in the avergage ion charge number Z from Z2 to Z1.Since ne is fixed, this necessitates a reduction of vertical lattice spacing by a factor Z /Z2 and thus, a formation of a contracted crystal (Paper I) with progressively increasing degree of contraction.
Infinite 3D contracted Coulomb crystals are unstable for contraction factors as high as 0.9-0.95(strictly speaking, this result is obtained in Paper I only for one-component crystals).Thus, one can expect breaking of an overcontracted crystal, its refreezing into a stress-free cubic configuration, and the process repeating itself several times before the mixed region is fully crossed.One can expect a formation of several interfaces of incommensurate crystals, whose strength will be determined by the weaker of them, i.e., presumably, by a contracted crystal, with the contraction factor slightly above the critical value and the lowest Z available.
In order to expose as vividly as possible the idea of crystal stretching in solidifying dense stellar matter in response to electron density or Z decrease, we have assumed the simplest possible layer-by-layer mode of the epitaxial growth.We recognize that other growth modes are possible, for instance, edgewise growth of certain special, in particular, close-packed planes, having a non-zero angle with the average crystal-liquid boundary (e.g.Billig 1955;Rosenberg & Tiller 1957;Tiller 1958;Bauser & Strunk 1984).However, the net macroscopic result of such growth is still a single 3D crystal, growing perpendicular to some general plane, constituting the average interface, which, in the case of NS and WD, must be accompanied by stretching to respect charge neutrality.Solidification of solutions adds yet another layer of complexity to the epitaxial process (e.g.Tiller et al. 1953;Burton, Prim & Slichter 1953;Trainor & Bartlett 1961;Mullins & Sekerka 1964;Tiller 1968).
To summarize, we propose, that, in contrast to the standard picture, which is based on comparison of energies of various crystal structures but neglects their growth kinetics, NS crusts and WD cores, for the most part, are made of stretched (elongated and contracted) crystals rather than of the cubic ones.This has several astrophysical implications, the most obvious being for elastic properties of matter.

Breaking strain of dense stretched matter
Let us analyze in some detail the strength of the stretched crystals with respect to shear deformations in planes orthogonal to the stretch.This seems to be relevant for NS and WD physics, where the stretch is aligned with the gravity, whereas the shear, caused, for instance, by magnetic field evolution, is horizontal and does not perturb hydrostatic equilibrium.
In order to find breaking strain, we looked for unstable phonon modes near4 the Brillouin zone center (i.e.modes with imaginary frequencies, see Paper I for details) for crystals, that were stretched by a factor ξ and then sheared in the perpendicular plane.In practice, we have selected two basis lattice vectors 1 and 2 in the plane perpendicular to the stretch direction.The third basis vector goes out of the plane.For stretch, its vertical component is multiplied by ξ, after which, for shear, its horizontal cartesian components are augmented by θ1 = θ cos χ (along the basis vector 1) and θ2 = θ sin χ.Cartesian components of these basis vectors and reciprocal lattice basis vectors as well as some other lattice parameters are summarized in Table 1.
One observes striking differences in breaking strain behaviour with ξ for different stretch directions.A very significant elongation is possible along the cube edge (Paper I).At ξ = √ 2, the lattice, that was originally bcc, turns into face-centred cubic (fcc).For stretches along the cube diagonal, the breaking strain drops abruptly for contractions but gradually for elongations.For stretches along the face diagonal, the situation is exactly opposite.and what is the proper way of deducing crust properties from those of perfect crystallites.A plausible model (polycrystalline model) may be to assume that, in a horizontal layer, there are crystallites stretched vertically by approximately the same factor ξ but oriented more or less randomly.Experiments on bicrystals show that crystals with different liquid-solid interface orientations can grow side by side (e.g.Rosenberg & Tiller 1957, fig. 3).Then, under a shearing deformation, the crystallites will have the same strain to maintain continuity, and the crystallite with the minimum breaking strain will fail first.This may result in a stability loss and failure of neighbouring crystallites8 .
In order to have a more representative breaking strain minimization, we have additionally considered 4 stretch directions of lower symmetry.Lattice properties for these cases, denoted aux 1-4, are also collected in Table 1.We note that the 7 considered cases of crystal growth differ, among other things, by interplane spacing η, which varies (for ξ = 1) from a l / √ 84 for auxiliary case 4 to a l / √ 2 for the face-diagonal stretch (cf.last column of Table 1).
The numerical data for the minimum breaking strain, obtained by minimizing over these 7 stretch directions, is shown in Fig. 3 by dots.To rms accuracy better than 0.1%, the data can be fitted as Breaking strain of 0.12-0.14for shear of an ideal bcc crystal at ξ = 1 has been obtained by Horowitz & Kadau (2009, fig. 1) in MD simulations.This value has been reduced to 0.1 to account for various imperfections and is currently used in applications as the breaking strain for shear deformations (see also discussion in Baiko & Chugunov 2018).Minimum breaking strain for shear of the ideal unstretched bcc crystal obtained in the present work [Fig.3 and equation ( 14)] is ≈ 0.064, which is ∼ 2 times lower than the range quoted above.Reducing it in the same way, we arrive at the updated value of ∼ 0.05 as the breaking shear strain of an imperfect unstretched bcc crystal.

Preferred orientation of growth and macrocrystallites
Liquid-phase epitaxy is a near-equilibrium process, which produces crystal layers of extremely high quality.Also of interest is the Czochralski crystal growth technique (Czochralski 1918;Teal & Little 1950;Buckley 1951;Hurley 1987;Scheel 2000), which nowadays yields ideal crystals of ∼ 2 m size.We note that these methods are subject to severe limitations posed by finite sizes of the apparatus and associated with them nonuniformities, thermal gradients, and stresses (e.g.Dash 1959).It seems then, that in a nearequilibrium crystallization with extremely uniform temperature and composition distributions and plenty of anneal On Earth, natural single crystals as large as ∼ 18 m (and possibly ∼ 50 m) have been found (Rickwood 1981).The abstract of this paper begins with a remarkable statement: "No upper limit on the size of crystals is to be expected . . ." The occurrence of so huge natural crystals indicates, at least, that there are robust self-consistent mechanisms of seeding their growth.
In an illuminating experiment of Rosenberg & Tiller (1957, figs. 6 and 8), purified lead was melted and poured into a mold designed for unidirectional (upward) freezing.Subsequent analysis of the ingot had shown that there were no crystallites other than those, which nucleated on the bottom boundary, and that the ingot bottom surface had 10 times as many crystallites as the top one.This means that 90% of the original crystallites have been crowded out by the surviving ones.A possible theoretical model of such crowding out has been proposed by Tiller (1957, esp. fig. 8).The preferred orientation of crystallite growth was close to perpendicular to {111} planes.With addition of a sufficient quantity of silver impurities, the preferred orientation switched to that perpendicular to {100} planes.However, the conditions of these experiments were clearly far from equilibrium.
According to Bravais's rule, in equilibrium, crystals tend to grow towards a shape bounded by the slowest growing planes (e.g.Trainor & Bartlett 1961).The fastest growing surfaces grow themselves out (e.g.Weymouth & Soepono 1962).The slowest growing planes are the close-packed ones.These are {110} planes for bcc crystals and {111} planes for fcc.The rule is supported by prominence of {111} plane growth in various experiments on Earth oftentimes conducted on fcc materials.Should we then expect that the entire crystallization front in a compact star grows perpendicular to {110} planes of the stretched bcc lattice or close to these directions?This is a distinct possibility9 .However, we note, that in a rare experiment on a bcc material sodium (Weymouth & Soepono 1962, fig. 1), no clear orientational preference has been observed with a slight tendency to directions perpendicular to {111}.
In the case of macroscopic crystallite formation (macrocrystallite model), particular properties of ideal crystallites may have relevance for astrophysical phenomena.For instance, if we assume {111} growth of a bcc lattice (with ξcrit ≈ 1.13), the critical elongation is achieved over the length of ∼ 0.13ne/∇ne, where and we have assumed that the pressure is dominated by ultrarelativistic electron gas.Furthermore, ρ is the mass density in units of g/cc, g is the gravitational acceleration in cm/s 2 , Pe is the electron pressure, and A is the ion mass number.At density 10 9 g/cc, it takes 2 m to grow to a critical height10 .
At ξ not too close to ξcrit, the elastic properties of matter differ quantitatively but not qualitatively from the conventional results.We see however, that the descending portion of the dashed red curve in Fig. 3 is approximately linear from ξ = 1 to ξcrit.This indicates the presence, upon freezing, of layers (with ξ ξcrit) whose breaking strain for shear is much lower than at ξ = 1.Specifically, there is a layer of 20 cm thick, in which the breaking strain is 10 times lower, and a layer of 2 cm thick, in which the breaking strain is 100 times lower, than for the bulk of the crystallite.Let us clarify that, even though these layers are substantially weaker than the unstretched material, they are stable and, by any conventional measure, extremely strong.They would not break spontaneously, but require for that a huge stress, which, however, is 1-2 orders of magnitude lower than at ξ ≈ 1. If, under the action of magnetic or any other such stresses late in a NS evolution, these layers break and then refreeze into a stress-free cubic configuration, they become 10 or 100 times stronger and thus, much harder to break at a later time.
If, on the other hand, the growth is of the face-diagonal type with ξcrit ≈ 1.06, the maximum crystallite height at the same density is ∼ 0.8 m.In the course of stretching, the crystal becomes some 20% stronger and then abruptly loses stability (cf.dot-dashed blue curve in Fig. 3).Due to extremely large derivative |dεcrit/dξ| near ξcrit for elongations, weak layers, if any, are much thinner.

Effective shear modulus
Another quantity of interest for applications is the effective shear modulus µ eff (ξ).In order to calculate it at arbitrary ξ, one could apply an infinitesimal shearing deformation to a crystal stretched in a particular direction, evaluate the energy difference δU , and then derive shear modulus µ for this stretch orientation from the formula: where V is the volume.Given a particular stretch direction, firstly, we have found µ eff averaged over the shearing deformation azimuthal angle χ.This quantity (in units of niUC, assuming single ion type with charge Ze) is shown in Fig. 4 as a function of ξ for 7 different stretch orientations.The variation in some cases is remarkable.For the cube-diagonal stretch (dashed red), µ {111} eff varies from 0.02 at ξ = 0.92 to 0.15 at ξ = 1.13.For cube-edge stretch (solid green), we reproduce classic results for S1212 for bcc (0.1827) and fcc (0.1852) lattices (Fuchs 1936).
Since strain for all crystallites is assumed to be the same, we can further average µ eff over 7 stretch orientations.This results in a slowly varying thin solid black curve in Fig. 4.This kind of averaging neglects the fact that some planes are nominally less numerous than others, for instance, there are only 6 original cubic seed orientations, that produce growth perpendicular to a {100} plane, whereas there are 48 equivalent {421} planes (aux 4).Qualitatively, such simplistic averaging may reflect the truth if the growth of low-index planes is, for whatever reason, preferred, or if they effectively represent larger solid angles.Optionally, we can treat all planes on the same footing and include the number of equivalent planes Neqv into the average shear modulus calculation 11 .Such weighted average results in the slightly different thin solid orange curve in Fig. 4. At ξ = 1, these approaches yield µ eff (1) ≈ 0.116 and 0.114, respectively, which, in both cases, is within 5% of the effective shear modulus µ OI eff = 0.1194 introduced by Ogata & Ichimaru (1990).
11 Neqv is given in the last column of Table 1.

Other possible effects
In this subsection, we shall list a few other possible effects of crystal stretching in NS crusts and WD cores, which have not been studied in any detail but may be relevant for astrophysics.
It is clear, that stretched crystals have higher groundstate energy, than the cubic ones.The difference between the two has been analyzed in some detail by Baiko & Chugunov (2018, figs. 5 and 6), where it was shown to be less than the latent heat associated with freezing into the stress-free cubic lattice.We thus expect, that freezing into a stretched configuration produces less latent heat at crystallization than the standard value.The rest of the energy, in principle, may be released at an arbitrarily late cooling stage.In the context of WD, this may happen at an age, when the stellar luminosity is much lower than at crystallization, and thus noticeably delay WD cooling.
Phonon thermodynamics of stretched crystals in general, and their heat capacity in particular may be subject to appreciable modifications.This is related to the expected appearance of a low-frequency mode, which eventually, at ξcrit, becomes unstable.This also may be important for cooling modelling.Phonon dispersion curves and polarisation vectors affect electron-phonon scattering rates, which, in turn, determine various electron kinetic coefficients such as thermal and electric conductivities.If there is a preferred crystal growth direction, this quantities will become anisotropic.This may be especially pronounced for crystals stretched perpendicular to {100} planes (where ξ can significantly exceed 1) or stretched close to their breaking limits.
Stretching affects nearest neighbour distance in a crystal.At fixed density, for the cube-diagonal stretch, this distance along the stretch direction is different from that in the unstretched lattice by a factor of ξ 2/3 .For the face-diagonal stretch, the nearest neighbour distance in the transverse direction gets multiplied by ξ −1/3 .This may be important for pycnonuclear reaction rates (e.g.Yakovlev et al. 2006), which are very sensitive to subtle details such as Coulomb barrier height and tunneling length, and may easily vary by many orders of magnitude.
Based on thermodynamic and evolutionary arguments, it is well established that NS crusts and WD cores contain adjacent regions with different ion compositions.At low tempertures, the regions crystallize.Boundary layer between these crystals is expected to be relatively thin so that transition between them occurs at essentially constant electron density ne.We have derived a general expression for the interaction energy of such Coulomb crystals in an idealized model.Under broad assumptions regarding the composition and lattice structure, the crystals are incommensurate, and in this case, the interaction energy reduces identically to zero.We hypothesize that, in reality, for both commensurate and incommensurate crystals, their interaction is actually controlled by interface roughness.Consequently, in spite of the interface existence, the combined system behaves as a single crystal, and its breaking properties are determined by those of the weaker part.
Using formulae obtained on the previous step, we have revisited near-equilibrium bottom-up crystallization in NS crusts and WD cores.We argue that this process is in many ways similar to liquid-phase epitaxial (i.e.preserving order of previous layers) crystal growth or crystal pulling from melt in terrestrial laboratories.Moreover, we believe that due to lack of electron bonds and, more generally, due to independence of angular orientation of newly arriving nuclei, due to long-range nature of Coulomb forces, and due to nearly ideally uniform environment of compact star interior, the epitaxial growth of Coulomb crystals in these astrophysical objects may be realized more easily than on Earth.
If this is true, lateral ion positions are locked by already frozen layers, whereas vertical ion positions are governed by charge neutrality.Therefore, interplane spacings in growing crystals either gradually increase, tracing ne decrease, as the crystallization front moves away from the stellar center, or decrease, tracing decrease of the average ion charge number, when the crystallization front crosses boundary between the regions with different compositions.This results in a formation of stretched (elongated and contracted) Coulomb lattices, as opposed to the standard assumption of cubic lattice formation, which is based on energetics arguments but does not take into account crystal growth kinetics.Overelongated and overcontracted crystals develop unstable phonon modes, lose stability, and get destructed.This limits the vertical sizes of growing crystallites.
We study breaking strain of stretched matter for shear deformations perpendicular to the crystal growth direction and find striking differences in its behaviour vs. stretch factor ξ for different stretch orientations.For unstretched material, we find a shear deformation with breaking strain ∼ 2 times lower than the currently known estimate (0.05 instead of 0.1 with account of crystal imperfections).Assuming the presence at a given depth of arbitrarily oriented crystallites stretched by the same factor (polycrystalline model), we minimize breaking shear strain over stretch directions at each ξ and fit the resulting ε min crit (ξ) dependence by a simple analytic formula.
Furthermore, we present some arguments in favor of macrocrystallite formation and preferred crystal growth orientation in dense matter of NS crusts and WD cores (macrocrystallite model).If this is the case, specific properties of crystallites may be directly connected with observed astrophysical phenomena.For instance, for growth parpendicular to {111} planes, we expect the appearance of crustal layers whose shear strength is 10-100 times lower than for unstretched matter of the same density.If, under the action of magnetic or any other stresses late in an NS evolution, these layers break and then refreeze into a stress-free cubic configuration, they become 10-100 times stronger and do not break that easily at a later time.On the contrary, for growth parpendicular to {110} planes, which are the closepacked planes of the bcc lattice, no such weak layers are expected.The diversity of layer sizes and strengths as well as the possibility of significant strengthening of matter after breaking may be related to rich magnetar burst and outburst phenomenology (e.g.Kaspi & Beloborodov 2017) and its extension to lower-B objects.
We have also calculated the dependence of the crystal shear modulus on ξ for growth perpendicular to 3 low-index planes and 4 planes of lower symmetry as well as the shear modulus averaged over growth orientations and fitted main results by analytic expressions.Several other effects of crystal stretching on thermodynamics, kinetics, and pycnonuclear reaction rates in NS crusts and WD cores, which may be important for applications, have been proposed.

Figure 1 .
Figure1.A schematic of the interface of two incommensurate crystals.Shaded blue region is the quasi-crystal zone.Lattice constants for both crystals are greatly exaggerated in comparison with the typical "wavelength" of the interface roughness.

Figure 2 .
Figure 2. Panels a-d: potential relief produced by a semi-infinite bcc Coulomb crystal split along {110}, {100}, {111}, or {221} crystallographic plane, respectively, vs. lateral coordinate d (in units of a l /2).Vertical distance from the top plane |d ⊥ | = |a 3⊥ | = η.Crosses indicate ion position on the next (newly forming) plane in a perfect 3D lattice.Potential minima are at or near the centres of the graphs, maxima are at the corners and/or at the graphs' upper edges.Panel e: the difference between the maxima and minima (solid) and U slab (η/2) (dashed) vs. log η.Panel f: potential of the {100} half-lattice (dashed) and combined potential of the half-lattice and the electron slab (solid) vs. |d ⊥ | with lateral coordinate, corresponding to a minimum (thin) or to a maximum (thick) in panel b.Vertical lines are the electron slab boundaries.Potential zero level is taken mid-slab, i.e. at |d ⊥ | = a l /2, at the lateral position of the minimum.

Figure 3 .
Figure 3. Breaking strain for transverse shear of stretched bcc Coulomb crystal vs. the stretch factor ξ.

5
Stretch direction (A) from Paper I 6 Stretch direction (B) from Paper I 7 in units of a l It is not known what is the actual structure of NS crust

Figure 4 .
Figure 4. Effective shear modulus of stretched bcc Coulomb crystal in units of n i U C .

Table 1 .
Horowitz & Kadau 2009)e, NS crusts and especially WD cores (which literally take eons to freeze) have a much better chance of forming large-scale near perfect crystallites (this possibility has been mentioned earlier byHorowitz & Kadau 2009).