Anisotropic crystal structure of magnetized neutron star crust

TL;DR

This paper investigates how strong magnetic fields cause anisotropic deformation in neutron star crust crystals, potentially triggering magnetar activity and leading to a transition to liquid crystal states under extreme conditions.

Contribution

It reveals that magnetic fields induce anisotropic pressure and deformation in neutron star crust crystals, identifying instability limits and proposing a transition to liquid crystal phases.

Findings

Magnetic fields cause anisotropic stretching or shrinking of crust crystals.

Phonon mode instability limits magnetic field variation without destroying the crystal.

Extreme pressure anisotropy may lead to a liquid crystal phase.

Abstract

Although crystallized neutron star crust is responsible for many fascinating observational phenomena, its actual microscopic structure in tremendous gravitational and magnetic fields is not understood. Here we show that in a non-uniform magnetic field, three-dimensional ionic Coulomb crystals comprising the crust may stretch or shrink while their electrostatic pressure becomes anisotropic. The pressure depends non-linearly on the magnitude of the stretch, so that a continuous magnetic field evolution may result in an abrupt crystal elongation or contraction. This may provide a trigger for magnetar activity. A phonon mode instability is revealed, which sets the limits of magnetic field variation beyond which the crystal is destroyed. These limits sometimes correspond to surprisingly large deformations. It is not known what happens to crust matter subject to a pressure anisotropy exceeding these limits. We hypothesize that the ion system then possesses a long-range order only in one or two dimensions, that is becomes a liquid crystal.

Equations

$$U(\xi ,u_{\alpha \beta })=U(\xi )+V{\sigma }_{\alpha \beta }(\xi )u_{\alpha \beta },$$

$${\sigma }_{\alpha \beta }=-P^{\mathrm{e}}{\delta }_{\alpha \beta }+S_{\alpha \beta }^{\mathrm{st}},$$

$$\frac{\partial }{\partial r_{\beta }}({\sigma }_{\alpha \beta }+{\tau }_{\alpha \beta })-\rho \frac{\partial \psi }{\partial r_{\alpha }}=0,$$

$$\frac{\partial P}{\partial z}$$

$$\frac{\partial P}{\partial x}$$

$$\frac{\partial }{\partial z}(P^{\mathrm{e}}+P_{\parallel })$$

$$\frac{\partial }{\partial x}(P^{\mathrm{e}}+P_{\perp })$$

$$\frac{\partial [P^{\mathrm{e}}(z,0)+P_b^{\mathrm{st}}(z,0)]}{\partial z}=-g\rho (z,0).$$

$$\delta P_{\rho }(z,x)+P_b^{\mathrm{st}}(z,x){\Delta }_{\perp }(\xi (z,x))=\frac{B^2(0)-B^2(x)}{8\pi },$$

$$\delta P_{\rho }(z,x)=P^{\mathrm{e}}(z,x)+P_b^{\mathrm{st}}(z,x)-P^{\mathrm{e}}(z,0)-P_b^{\mathrm{st}}(z,0)$$

$$\frac{\partial \delta P_{\rho }}{\partial z}=-\frac{g}{1+C}\delta \rho ,$$

$${\Delta }_{\perp }(\xi (z,x))=\frac{B^2(0)-B^2(x)}{8\pi P_b^{\mathrm{st}}(z)}.$$

$$B(\delta B{)}_{\mathrm{crit}}\approx 2.2\times 10^{21}{\Delta }_{0.1}{\rho }_6^{4/3}{Z}_{26}^2{A}_{56}^{-4/3}[{\mathrm{G}}^2],$$

Full text

††pagerange: LABEL:firstpage–LABEL:lastpage††pubyear: 2013

1 Introduction

Neutron stars are magnificent astrophysical objects for a wide range of stunning observational phenomena associated with them as well as for extraordinary challenges and opportunities they present to a theorist trying to explain their behavior (e.g., Kaspi, 2010). Among the most striking manifestations of neutron stars, one can mention pulsars, amazing machines converting the star rotation energy into multi-wavelength highly periodic radiation and accelerated particles; millisecond pulsars, the most accurate clocks in the Universe; soft-gamma repeaters, famous for their restless bursting activity and mind-boggling giant flares.

From a theorist’s point of view, the most exciting is probably the core of a neutron star, which is compressed to super-nuclear densities by immense gravity. Essentially, this is a unique laboratory of strong interaction physics (e.g., Haensel, Potekhin & Yakovlev, 2007). The core is hidden from us by an $\sim 1$ km thick crust which may be threaded by magnetic fields of up to at least $\sim 10^{15}$ G. This underscores the need for a microscopic model of magnetized neutron star crust to enable extraction of information coming from the core. Of great theoretical interest is also the origin of such colossal fields, whereas their interaction and joint evolution with the crust are believed to be chiefly responsible for the observed phenomena.

The density $\rho$ of the crust matter spans many orders of magnitude. One usually distinguishes the outer crust at $\rho<4.3\times 10^{11}$ g cm*-3* and the inner crust at densities above that and up to $\sim 1.5\times 10^{14}$ g cm*-3* where the core starts. The bulk of the outer crust consists of fully ionized atomic nuclei of various sorts and strongly degenerate nearly incompressible electron gas. In the inner crust in addition to ions and electrons there are degenerate neutrons dripped from the nuclei. In view of a rapid gravitational separation, it is typically assumed that at any given density in the crust there is only one ion species. If any possible but not firmly established effects associated with dripped neutrons are neglected, such electron-ion system is governed by pure Coulomb forces and is known to crystallize into a Coulomb crystal at a sufficiently low temperature (, Nagara, Nagata & Nakamura1987). This occurs at very early stages of the neutron star life across the entire crust with the exception of the outermost low-density layers.

It is well-known that the body-centered cubic (bcc) lattice has the lowest electrostatic energy among all one-component Coulomb crystals with uniform electron background, however, the energy difference between several stable structures is tiny. Previous studies of Coulomb solids have thus focused on a handful of lattices characterized by near minimum energies and with a possible presence of a uniform magnetic field (e.g., Cohen & Keffer, 1955; Carr, 1961; Nagai & Fukuyama, 1982, 1983; , Baldereschi, Senatore & Oriani1992, Baiko, Potekhin & Yakovlev2001; Kozhberov & Baiko, 2015; Chamel & Fantina, 2016). Formation of the bcc crystal in neutron star envelopes has been confirmed in recent molecular dynamics simulations (, Engstrom, Yoder & Crespi2016). Linear elasticity theory based on the bcc lattice has been used to study neutron star crust oscillations (e.g., Piro, 2005; Sotani, 2016). Molecular dynamics simulations of bcc crystal breaking have been performed (e.g., Horowitz & Kadau, 2009). However, as we demonstrate below, these approaches miss a significant part of the story.

2 Anisotropic crystal pressure

Consider a bcc Coulomb crystal of ions with rigid and uniform electron background in the outer neutron star crust. Suppose the crystal is stretched or shrunken in some direction by a stretch factor $\xi$, which means that projections of all lattice vectors on this direction are multiplied by $\xi$ while the overall scale factor is adjusted to maintain constant ion number density. The system then remains a perfect crystal. If one neglects the ion motion about the lattice nodes and the electron polarisation which are typically valid zero-order assumptions, the system energy $U(\xi)$ is the sum of the energy of strongly degenerate electron gas and the electrostatic (Madelung) energy of ions.

Let us denote lattice vectors of the stretched crystal as $\bm{R}$ and apply a uniform infinitesimal deformation $R_{\alpha}\to R_{\alpha}+u_{\alpha\beta}R_{\beta}$. Then, to first order in $u_{\alpha\beta}$, the energy of the deformed crystal is

[TABLE]

where $V$ is the volume and $\sigma_{\alpha\beta}(\xi)$ is the stress tensor of the crystal stretched by $\xi$. The stress tensor consists of two terms

[TABLE]

where $P^{\rm e}$ is the degenerate electrons’ pressure (e.g., Landau & Lifshitz, 1980; Potekhin & Yakovlev, 2012) and $S^{\rm st}_{\alpha\beta}$ is the electrostatic (ion) contribution. For bcc lattice the latter contribution is isotropic $S^{\rm st}_{\alpha\beta}(\xi=1)=-P^{\rm st}_{b}\delta_{\alpha\beta}$. The electrostatic pressure of the bcc lattice (Fuchs, 1936) is $P^{\rm st}_{b}=-\zeta nZ^{2}e^{2}/(3a)<0$, where $Z|e|$ is the ion charge, $n$ is the ion number density, $\zeta\approx 0.895929255682$ is the Madelung constant of the bcc lattice, and $a=(4\pi n/3)^{-1/3}$ is the ion sphere radius.

For $\xi\neq 1$, one can calculate the ion contribution to the stress tensor exactly using a practical formula (Fuchs, 1936; Baiko, 2011). The pressure is no longer isotropic. In Fig. 1, panel (A), we stretch the crystal along the direction towards the nearest neighbour i.e. along the diagonal of the main lattice cube. In panel (B) the stretch is along the direction towards the second nearest neighbour i.e. along the edge of the main lattice cube. In both cases the pressure is isotropic in the plane perpendicular to the stretch direction. We denote it $P_{\perp}$ and show by dot-dashed (blue) lines in Fig. 1. The pressure along the stretch direction, $P_{\parallel}$, is shown by solid (red) curves and is clearly not the same as $P_{\perp}$. The dependence of pressure components on $\xi$ is seen to be non-linear, especially in panel (B), which limits the accuracy of a standard elastic expansion starting at $\xi=1$.

The stretch in panel (B) at $\xi=\sqrt{2}$ is known (e.g., Nagai & Fukuyama, 1982) to produce a face-centered cubic (fcc) lattice, which is isotropic. The stretch in panel (A) also produces isotropic structures: a simple cubic lattice at $\xi=2$ and the fcc lattice at $\xi=4$.

It is worth noting that the two stretch directions considered coincide with the highest symmetry axes of the bcc crystal. Stretching the crystal along the direction towards the third nearest neighbour i.e. along the diagonal of the face of the main lattice cube yields the pressure tensor with three different eigenvalues. Even though further work is desirable to render a more complete understanding of Coulomb crystals’ behavior under stretches, cases (A) and (B) will be sufficient for our discussion. Also, other types of deformations (e.g., shear) can be studied in a similar fashion.

3 Equilibrium in non-uniform magnetic field

External forces of various nature can cause crystal stretching just discussed. One possibility, which seems to be especially relevant for neutron star crusts, is the action of a non-uniform magnetic field. For illustration of expected effects, let us construct a microscopic model of local hydrostatic equilibrium of crystallized neutron star crust matter in a stationary magnetic field $\bm{B}(\bm{r})$. The momentum conservation equation (e.g., Carroll et al., 1986) can be written as

[TABLE]

where $\tau_{\alpha\beta}=(2B_{\alpha}B_{\beta}-B^{2}\delta_{\alpha\beta})/(8\pi)$, $-\nabla\psi=\bm{g}$ is the gravitational acceleration, and $\rho$ is the mass density.

Of primary interest for us are magnetic field configurations which cannot be in hydrostatic equilibrium with matter characterized by isotropic pressure, i.e. which cannot be stationary in a liquid star. Such configurations may be produced upon freezing of a non-stationary liquid system or be a result of magnetic field evolution already in the solid phase.

Consider the simplest geometry of magnetic field and gravity: $\bm{B}=B(x)\bm{e}_{z}$ and $\bm{g}=-g\bm{e}_{z}$ ($g$ is constant), for which there is no hydrostatic equilibrium in a liquid. To prove the latter statement we observe that in this case $-\tau_{xx}=-\tau_{yy}=\tau_{zz}=B^{2}/(8\pi)$, and if the total pressure $P$ is strictly isotropic, $\sigma_{xx}=\sigma_{yy}=\sigma_{zz}=-P(\rho)$. Then Eq. (3) reduces to

[TABLE]

Equation (5) implies that $P(\rho)=f(z)-B^{2}(x)/(8\pi)$, where $f(z)$ is some function, and thus $\rho$ depends on both, $x$ and $z$. However, Eq. (4) can be written as $\rho=-f^{\prime}(z)/g$ which is independent of $x$. The contradiction means that the stress tensor of matter must be anisotropic and the easiest way to realize this is to impose a crystal stretch.

The stretch direction and magnitude depend on specific boundary conditions. Here we assume, somewhat arbitrarily, that the stretch is along the $z$-axis and $\xi(z,x=0)=1$. Then $\sigma_{zz}=-P_{\parallel}-P^{\rm e}$ and $\sigma_{xx}=-P_{\perp}-P^{\rm e}$, where $P_{\parallel,\perp}$ are determined by both mass density and the stretch factor, while $P^{\rm e}$ is a function of $\rho$ only. Consequently, Eq. (3) can be rewritten as

[TABLE]

Let us introduce functions $\Delta_{\parallel,\perp}(\xi)$ as $P_{\parallel,\perp}=[1+\Delta_{\parallel,\perp}(\xi)]P^{\rm st}_{b}(\rho)$ so that $\Delta_{\parallel,\perp}(1)=0$ (cf. Fig. 1). At $x=0$, Eq. (6) reduces to the equation of hydrostatic equilibrium in the absence of the magnetic field:

[TABLE]

Integrating Eq. (7) we obtain

[TABLE]

where

[TABLE]

is the deviation of the isotropic contribution to the pressure from its field-free value at given $z$. The two terms on the left-hand side of Eq. (9) thus express the fact that the magnetic pressure is compensated by a combined action of a mass density change and a crystal stretch.

To take this consideration a step further, we assume that $\Delta_{\parallel}(\xi)=-C\Delta_{\perp}(\xi)$, where $C\sim 1$ is a positive constant. In particular, $\Delta_{\parallel}$ and $\Delta_{\perp}$ are proportional to each other for a linear dependence of $P_{\parallel,\perp}-P^{\rm st}_{b}$ on $\xi-1$ (cf. Fig. 1). Then we can differentiate Eq. (9) with respect to $z$ and use the derivative of the second term on the left-hand side to express the derivative of $P_{\parallel}$ in Eq. (6). This results in

[TABLE]

where $\delta\rho(z,x)=\rho(z,x)-\rho(z,0)$. Suppose, for certainty, that $B^{2}(x)<B^{2}(0)$. Then $\delta\rho,\delta P_{\rho}\geq 0$ and since $P^{\rm st}_{b}<0$, $\Delta_{\perp}\leq 0$, i.e. $\xi\geq 1$ (cf. Fig. 1). According to Eqs. (11) and (9), $\delta P_{\rho}$ decreases with increase of $z$, whereas $\xi$ increases, so that the relative importance of the first term on the left-hand side of Eq. (9) diminishes with $z$. The actual value of $\delta P_{\rho}$ depends on the assumed boundary condition for $\delta\rho$ at a fixed $z$.

If $\delta\rho=0$ at some $z_{0}$ it will remain zero at $z>z_{0}$ and so will $\delta P_{\rho}$. In this situation, $\rho$, $P^{\rm e}$, and $P^{\rm st}_{b}$ are independent of $x$ and the magnetic pressure is balanced solely by the crystal stretch:

[TABLE]

The stretch factor $\xi(z,x)$ can be found by applying a well-defined function $\Delta_{\perp}^{-1}$ to this formula. Equation (12) represents the maximum effect of the magnetic field on the crystal shape and will be adopted for the sake of the qualitative discussion in the following Section.

4 Limits of crystal anisotropy

Given a Coulomb crystal with an arbitrary lattice, one can calculate its dynamic matrix using the standard formulas (Cohen & Keffer, 1955; Baiko, 2002) and analyse the crystal phonon modes. A bcc Coulomb crystal ($\xi=1$) has three phonon modes at each wavevector $\bm{k}$ in the first Brillouin zone, two of which are transverse acoustic, while the third one is longitudinal optic.

If $\xi\neq 1$ but is close to 1, the same picture persists. However, we have found that as the stretch or contraction increased, one of the acoustic modes developed an instability manifested by the appearance of imaginary frequencies at certain isolated directions of $\bm{k}$ in the vicinity of $k=0$ (i.e., at large wavelengths). In general, such an instability signifies a crystal destruction. The critical $\xi$ (denoted $\xi_{\rm crit}$) are marked by thin solid vertical lines in Fig. 1.

In Fig. 1(A) we observe a weakly non-linear dependence of the pressure tensor components on $\xi$ and the instability occurs at $\xi_{\rm crit}\approx 0.92$ or 1.13 when the anisotropy reaches $\sim 10\%$. Equation (12) applied at $\xi_{\rm crit}$ enables one to estimate magnetic field variation $(\delta B)_{\rm crit}$ which can be supported by the crystal as

[TABLE]

where we have used the explicit expression for $P^{\rm st}_{b}$, $\Delta_{0.1}\equiv\Delta_{\perp}(\xi_{\rm crit})/0.1$, $\rho_{6}$ is the mass density in units of $10^{6}$ g cm*-3*, $Z_{26}\equiv Z/26$, $A_{56}\equiv A/56$, and $A$ is the ion mass number. Hence, for iron at $10^{6}$ g cm*-3* and a typical pulsar field of $10^{12}$ G, $(\delta B)_{\rm crit}/B\sim 2\times 10^{-3}$. For a typical magnetar field of $4\times 10^{14}$ G, $(\delta B)_{\rm crit}/B\sim 10^{-8}$. Assuming a specific magnetic field geometry, one can deduce the maximum size of a perfect crystal, which may be present in such matter, as the distance over which magnetic field varies by $(\delta B)_{\rm crit}$. For instance, in the geometry where $\delta B\sim B$ over the star’s radius $\sim 10^{6}$ cm, there may be $\sim 20$ m crystals in pulsars but only $\sim 0.1$ mm crystals in magnetars at $\rho=10^{6}$ g cm*-3*. If there is a “spot” with $\delta B\sim B$ over $10^{4}$ cm, then the crystals available will be 100 times smaller.

At higher densities a perfect crystal can withstand much stronger field variations. If the dripped neutrons in the inner crust simply add another isotropic term to the total stress, setting $\rho=10^{14}$ g cm*-3*, $Z=40$, and $A=1000$, we obtain $B(\delta B)_{\rm crit}\approx 5.3\times 10^{30}$ G2. It follows that $(\delta B)_{\rm crit}\sim B$ can be accommodated by the perfect crust at $B\lesssim 2\times 10^{15}$ G. Incidentally, this is the same number as the maximum magnetic field presently reported for magnetars (Olausen & Kaspi, 2014) and it is much smaller than the field required for complete inner crust magnetization (Baiko & Yakovlev, 2013).

A surprising result is illustrated in Fig. 1(B). In this case the instability occurs very quickly upon a crystal contraction ($\xi_{\rm crit}\approx 0.94$), but if the crystal is elongated, the pressure stays more or less isotropic (anisotropy does not exceed $2\%$) and the structure remains stable up to $\xi_{\rm crit}\approx 1.5$. The range of stable structures extends all the way to the fcc lattice and beyond. This is in contrast to Fig. 1(A), where the crystal would be destroyed before reaching fcc at $\xi=4$.

From this picture it appears that the instability is not caused by the stretch itself but by the pressure anisotropy associated with it. The value of the available strain is remarkable being several times greater than the typical value of $0.1$, obtained by molecular dynamics simulations of breaking strain in a shearing crystal (Horowitz & Kadau, 2009; Chugunov & Horowitz, 2010, 2012; Johnson-McDaniel & Owen, 2013), even though the respective pressure anisotropy is very weak. It would require $B(\delta B)_{\rm crit}$ approximately 10 times smaller than the estimate (13) in view of different $y$-axis scales in Figs. 1(A) and (B). Comparing critical stresses and strains in Figs. 1(A) and (B) one notices their pronounced dependence on the stretch direction. This calls into doubt the applicability of phenomenological von Mises and Tresca criteria for crystal breaking often used in neutron star crust contexts (e.g., , Ushomirsky, Cutler & Bildsten2000).

Suppose that the crystal is oriented in such a way that a magnetic field variation results in an increase of $\xi$ in the direction analysed in panel (B). This can continue till the minimum of $-P_{\perp}$ is reached but any further variation of $B$ in the same direction (e.g., as a result of field drift or decay) would require a jump of the stretch factor shown schematically in panel (B) by an arrow. Such a strong jump from $\xi\approx 1.08$ to $\approx 1.46$ coupled with a structural transition from a stretched bcc to a stretched fcc lattice might provoke a dramatic crust rearrangement, which could affect the stability of neighbouring crust regions and serve as a trigger mechanism for magnetar activity.

A natural question arises: what happens with the crust if the pressure anisotropy exceeds the values corresponding to $\xi_{\rm crit}$? It is not excluded that another three-dimensional crystal structure is stronger than the bcc lattice but a really substantial improvement seems unlikely. In principle, at low temperatures and strongly anisotropic pressures, the Coulomb system can be in a new, unknown state. To investigate this possibility we have performed a numerical experiment. If we artificially take $\xi$ in the crystal dynamic matrix formula beyond its critical value, the phonon instability will spread over finite ranges of propagation angles and to shorter wavelengths. This is illustrated in Fig. 2 where panels (A) and (B) correspond to super-critical crystal elongation and contraction, respectively, in the direction of the nearest neighbour [as in Fig. 1(A)]. The propagation angles are defined in a spherical reference frame with the $z$-axis in the direction of the stretch. Stars represent the unstable directions at $\xi=1.14$ (A) or $\xi=0.918$ (B), which is very close to the instability onset [cf. Fig. 1(A)]. It turns out, that in panel (A) these are directions towards the third nearest neighbours in the plane perpendicular to the stretch, but in panel (B) they do not coincide with any obvious axes of symmetry. Sequentially expanding contours encircle unstable regions at $\xi=1.15,1.2,1.25,1.3,1.4,1.5$ in panel (A) and at $\xi=0.91,0.9,0.85,0.8,0.7,0.5$ in panel (B). Shading emphasizes maximum unstable configurations attained in Fig. 2.

The white regions correspond to propagation angles at which all three phonon modes remain stable all the way down to $k=0$. In general, these domains shrink with increase of $|\xi-\xi_{\rm crit}|$. At first, they can be described as bundles of planes in $\bm{k}$-space containing the origin as in panel (A) for $\xi\lesssim 1.25$ and in panel (B) for all $\xi$. At higher stretches they may reduce to bundles of lines through the origin as in panel (A).

The unstable branches of the phonon spectrum look similar to collective modes of a Coulomb liquid, which disappear at finite $k$ and do not reach $k=0$ (e.g., Schmidt et al., 1997). This suggests a loss of the long-range order in the respective directions in the crystal. The persistence of the long-wavelength transverse modes in other directions even for relatively large deformations implies, by extension, that the long-range order in one or two dimensions may be preserved. Together these observations can be taken as an indication that at strongly anisotropic pressure the three-dimensional Coulomb system behaves as a liquid crystal.

5 Discussion

In neutron star crust, a bcc Coulomb solid will stretch or contract in response to a wide class of non-uniform magnetic fields. The ensuing crystal will have anisotropic electrostatic pressure, which may depend on the stretch factor in a very non-linear fashion. In some geometries, this may result in a strongly discontinuous response of the crystal structure to a continuous variation of the applied (magnetic) stress and serve as a trigger mechanism for magnetar activity. The three-dimensional Coulomb crystal develops a phonon mode instability (and fails) if the degree of the pressure anisotropy reaches $\lesssim 10\%$ which may correspond to a stretch factor as high as $\sim 150\%$. A respective estimate of the maximum magnetic field sustainable by a perfectly crystallized crust is consistent with the maximum soft gamma repeater field currently inferred from observations. It is unknown what state the crust is in at higher pressure anisotropies. One possibility alluded to here is that it resembles a liquid crystal having a long-range order in one or two dimensions only. Molecular dynamics simulations of Coulomb systems under anisotropic pressure may prove instrumental in discovering these structures.

Anisotropic Coulomb crystals and liquid Coulomb crystals (if exist) will have thermodynamic, kinetic, and elastic properties very different from those of the standard bcc. Such diverse microscopic quantities as neutrino emissivities, diffusion rates, plasma screening factors for nuclear reactions will be seriously affected. The strongest effect on the neutron star structure is expected near the star’s surface, while global characteristics, e.g., the mass–radius relation, most likely will not change due to the relative smallness of the electrostatic pressure.

Acknowledgments

We are grateful to A.I. Chugunov and D.G. Yakovlev for discussions. A.A.K. thanks Leading Science School 9297.2016.2 for support.