text:
# Classical generalized constant coupling model for geometrically frustrated antiferromagnets

## Abstract

A generalized constant coupling approximation for classical geometrically frustrated antiferromagnets is presented. Starting from a frustrated unit we introduce the interactions with the surrounding units in terms of an internal effective field which is fixed by a self consistency condition. Results for the magnetic susceptibility and specific heat are compared with Monte Carlo data for the classical Heisenberg model for the pyrochlore and *kagomé* lattices. The predictions for the susceptibility are found to be essentially exact, and the corresponding predictions for the specific heat are found to be in very good agreement with the Monte Carlo results.

*Introduction.-* In the last several years, geometrically frustrated antiferromagnets (GFAF) have emerged as a new class of magnetic materials with uncommon physical properties and have received a vast amount of attention (see hfm.2000  and references therein). In these materials, the elementary unit of the magnetic structure is the triangle, which makes it impossible to satisfy all the antiferromagnetic bonds at the same time, with the result of a macroscopically degenerated ground state. Examples of GFAF are the pyrochlore and the *kagomé* lattices. In the former, the magnetic ions occupy the corners of a 3D arrangement of corner sharing tetrahedra; in the later, the magnetic ions occupy the corners of a 2D arrangement of corner sharing triangles (see Fig. 1). In the case of materials which crystallize in the pyrochlore structure, the magnetic susceptibility follows the Curie–Weiss law down to temperatures well below the Curie–Weiss temperature. At this point, usually of two orders of magnitude smaller than the Neél point predicted by the standard mean field (MF) theory, some systems exhibit some kind of long range order (LRO), whereas others experience a transition to a spin glass state (SG). This is a striking feature for a system with only a marginal amount of disorder. Finally, there some pyrochlores which do not exhibit any form of order whatsoever, and are usually regarded as spin liquids. In the case of the *kagomé* lattice, even though there are very few real systems where this structure is realized, the magnetic properties fall in two great categories: the vast majority of the compounds studied show a transition to a LRO state with non collinear configuration of spins, and a few systems exhibit no LRO, but a SG like transition.

For these reasons, it is easy to understand the large amount of attention, both from the experimental and the theoretical point of view, these systems have attracted. From the theoretical point of view, a number of techniques have been used to try to understand the origin of the puzzling properties mentioned above. All the theoretical results seem to indicate that for these geometries the classical Heisenberg model with only nearest neighbor interactions does not order at any finite temperature, in agreement with Monte Carlo (MC) results moessner.1999 ; moessner.1998 ; Reimers.1992 ; Chalker.1992 ; Reimers.1993 ; Garanin.1999 . There are also a relatively few works which have dealt with the quantum effects in these systems harris.1991 . However, the main interest during last years has been in the classical GFAF.

In any case, the current theories have to deal with two main problems: the first one is that in real systems there are additional effects, which can be dipole–dipole interactions, small anisotropies, or further neighbor interactions, which are not easy to include in the theory, and so make it very difficult to compare with experimental results. However, this fact can be circumvented by using MC data calculated for the particular assumptions of the model. The second problem is that these models only provide qualitative agreement even when compared with MC data. In this sense, the single unit approximation, is the model which probably gives the best quantitative agreement when compared with susceptibility MC data for the pyrochlore lattice moessner.1999 . However, in order to reach such a good agreement, it is necessary to introduce an *ad hoc* rescaling procedure with an unclear physical justification. In the other hand, the infinite component spin vector model provides qualitative agreement for the *kagomé* case Garanin.1999 .

Obviously, the difficulty to find reliable models is rooted in the complexity of the geometry of these systems. Therefore, we think that, before embarking on the development of complicated models (which receive small support even from MC data obtained for the particular assumptions of the model), we should try to find simpler reliable approximation schemes which, as a first step, are in good qualitative and quantitative agreement with MC data. At the same time, these models should be as simple as possible (mathematically speaking), in order to be easily generalized to include additional interactions that play an important role in real systems. Only in this way, can we compare the models with data obtained from real systems, without this comparison being obscured by the complexity of the mathematical formalism.

In this work, we deal with the first of the two steps mentioned in the previous paragraph. We introduce a generalization of the well known constant coupling (CC) method kastel.1956 , which has been successfully applied to the study of 3D standard ferro and antiferromagnets. Our method, which we will call generalized constant coupling (GCC) method, takes into account the geometrical peculiarities of the frustrated structures mentioned above. In spite of its mathematical simplicity, the susceptibilities and specific heat calculated for the pyrochlore and *kagomé* lattices, in the paramagnetic region, are in excellent agreement with MC data for the classical Heisenberg model in these systems.

It is important to note that, even though we have focused the discussion of the present work in the context of geometrical frustration in Heisenberg systems, this approach represents, in fact, a general technique for the investigation of the thermodynamic properties of spin Hamiltonians in frustrated geometries and, thus, its applicability is not limited to the GFAF, but could be useful in the investigation of other physical systems where geometrical frustration is relevant (see ref. hfm.2000,  for a discussion of other problems where geometrical frustration is relevant).

*The model.-* The Heisenberg model with only nearest neighbor (NN) interactions in the presence of a magnetic applied field $`H_0`$ is described by the Hamiltonian smart.1966

$$H=J\underset{i,j}{}\stackrel{}{s}_i\stackrel{}{s}_jH_0\underset{i}{}s_{z_i},$$
(1)
where $`J`$ is the positive antiferromagnetic coupling, $`\stackrel{}{s}_i`$ and $`\stackrel{}{s}_j`$ represent classical spins of modulus $`s_0`$ located in a pyrochlore or *kagomé* lattice and $`s_{z_i}`$ the corresponding component along the applied field, and the sum is done over pairs of NN.

The idea of our approximate method is based on the experimental fact that the spin-spin correlations of the classical Heisenberg in the pyrochlore lattice are always short ranged moessner.1999 . Therefore, it is a reasonable approximation to start by considering isolated units (tetrahedra or triangles, for the pyrochlore and *kagomé* lattices, respectively), and later add the interactions with the surrounding units by in an approximate way (this is in contrast with the standard CC method, in which the elementary magnetic unit is taken as a pair of magnetic ions). Thus, it is important to first study the properties of the individual units. This task has been carried out by Moessner and Berlinsky moessner.1999 , and the result obtained for the partition function of an isolated unit with $`p`$ spins in the absence of applied field is given by

$$Z_p=\frac{(4\pi )^{p+1}}{(4\pi \beta \stackrel{~}{J})^{3/2}}_0^{\mathrm{}}x^2\left(\frac{\mathrm{sin}x}{x}\right)^p\mathrm{exp}\left(\frac{x^2}{2\beta \stackrel{~}{J}}\right)𝑑x,$$
(2)
where $`\beta =1/T`$ (we will use units of the Boltzmann constant $`k_B`$ throughout this work, so the energies are expressed in absolute temperature units), and we have introduced the effective coupling $`\stackrel{~}{J}=Js_0^2`$. Expressions for units with 2, 3, and 4 ions can be found in the original work by those authors.

The susceptibility per spin, $`\chi _p`$, can be easily calculated from the partition function by using the fluctuation–dissipation theorem, and is given by

$$\chi _p(\stackrel{~}{T})=\frac{S_p^2}{3pT}=\frac{2\stackrel{~}{T}}{3pJ}\frac{}{\stackrel{~}{T}}\mathrm{ln}Z_p,$$
(3)
where $`S_p^2`$ represents the thermal average of the total spin of the unit with $`p`$ ions, and we have introduced the adimensional temperature $`\stackrel{~}{T}=\frac{T}{\stackrel{~}{J}}`$.

Next, we introduce the interaction with neighboring units as an unknown internal effective field, $`H_1`$, created by the $`(p1)`$ NN ions outside the unit. The CC approximation consists of taking this internal field as $`H_1=(p1)H^{}`$, where $`H^{}`$ is the average internal field acting on a spin due to each of its NN. For example, in the case of the standard CW model, the internal effective field is given by $`H_1=2(p1)H^{}`$, as the ions are considered separately, and each has $`2(p1)`$ NN in the corner sharing structures considered in this work (see fig. 1). The internal field is evaluated by imposing the self consistent condition of equating the magnetization per spin in the field with that of a unit in the field, which can be mathematically stated as

$$s_0L(s_0(H_0+2(p1)H^{})/T)=\frac{m_p(H_0+(p1)H^{})}{p},$$
(4)
where the left side of the equation corresponds to the value of the magnetization per spin in the classical limit of the Curie–Weiss model, with $`L(x)=\mathrm{coth}x\frac{1}{x}`$ the so called Langevin function smart.1966 . The right side corresponds to the magnetization per spin of the isolated unit under the influence of the internal field. In the general case, equation (4) can only be solved numerically. However, in the paramagnetic regime observation , we can put, for small fields, $`\left(H_0+H_1\right)/T1`$,

$$\frac{m_p(H_0+(p1)H^{})}{p}\chi _p(T)(H_0+(p1)H^{}),$$
(5)
for the magnetization per spin in the cluster, where the susceptibility per spin is given by expression (3), and

$$s_0L(s_0(H_0+2(p1)H^{})/T)\frac{s_0^2}{3T}(H_0+2(p1)H^{}).$$
(6)
Taking into account eqs. (5) and (6), we can solve equation (4) very easily in terms of the function

$$\epsilon _p(\stackrel{~}{T})=\frac{S_q^2}{ps_0^2}1=\frac{2\stackrel{~}{T}^2}{p}\frac{}{\stackrel{~}{T}}\mathrm{ln}Z_q(\stackrel{~}{T})1,$$
(7)
to give $`H^{}=\frac{\epsilon _p(\stackrel{~}{T})}{(p1)\left[1\epsilon _p(\stackrel{~}{T})\right]}H_0`$, and the corresponding susceptibility per ion is then given by

$$\chi _p^{gcc}(\stackrel{~}{T})=\frac{1}{3\stackrel{~}{T}}\frac{1+\epsilon _p(\stackrel{~}{T})}{1\epsilon _p(\stackrel{~}{T})}.$$
(8)
Another interesting quantity we can readily evaluate in this model is the specific heat in zero applied field. In this limit, the internal field will be identically zero, and the specific heat will be identical to the specific heat calculated for non interacting units. Taking into account the relation $`\stackrel{}{s}_i\stackrel{}{s}_j=\frac{s_0^2}{p1}\epsilon _p(\stackrel{~}{T})`$ it can be easily shown that the specific heat per spin is given by

$$C_p=\frac{}{\stackrel{~}{T}}\epsilon _p(\stackrel{~}{T}).$$
(9)
*Comparison with MC data.-* Let us first consider the susceptibility obtained from the classical GCC model for the pyrochlore lattice. As can be seen from the observation of Fig. 2(a), the agreement between the theoretical curve and MC data moessner.1999  can be considered as exact in whole the temperature range. Moreover, at $`T=0`$ K, the value predicted by our model is $`1/8\stackrel{~}{J}`$, which is the value obtained from MC calculations. In the inset of that figure, we can see how the theory even predicts the kink in the susceptibility as one approaches very small temperatures.

In oder to corroborate that this surprising agreement is not a coincidence, we also compared the susceptibility calculated from the GCC model with MC data for the *kagomé* lattice Reimers.1993 , which can be seen in Fig. 2(b). Again, the curve predicted by our model can be considered as exact in whole the temperature range. Moreover, the value predicted by this model at $`T=0`$ K is $`1/6\stackrel{~}{J}`$, again the exact value.

It is important to stress that the curves depicted in Figs. 2(a) and 2(b) *have not been rescaled in any sense, nor contain any fit parameter*. As far as we know, no other model is able to predict a such accurate susceptibility in the whole temperature range, considering both the pyrochlore and the *kagomé* lattice on an equal footing.

In order to further check the accuracy of the model, we decided to compare its predictions for the specific heat. The results are depicted in Fig. 3, together with MC data Reimers.1992 ; Chalker.1992  for these systems. Again, the agreement is excellent. For the pyrochlore lattice, the calculated specific heat is essentially exact down to 0 K, where a value of 3/4 is obtained moessner.1998 . However, for the *kagomé* lattice, the exact value is $`11/12`$, whereas our model predicts a value of 1. In any case, this deviation is expected to occur, as our model does not include the effect of zero modes, which are known to be especially important in the *kagomé* case, where they give rise to the phenomenon known as order by disorder Chalker.1992 ; harris.1991 . In contrast, MC data suggest that this mechanism is absent in the Heisenberg pyrochlore.

*Conclusions.-* In this work we have presented a generalization of the so called constant coupling approach for geometrically frustrated antiferromagnets with a pyrochlore or *kagomé* structure. As a preliminary approach, we have considered the most simple description of these systems, which consists of a classical Heisenberg model with only nearest neighbors interactions. It seems clear from Monte Carlo results that, in this approximation, these systems remain paramagnetic in the whole temperature range. Analytical expressions for the static magnetic susceptibility and specific heat can be obtained in the framework of this model. These have been compared with Monte Carlo data for these systems. The predicted susceptibility is essentially exact in whole the temperature range for both types of lattice. The specific heat for the pyrochlore lattice is again found to be essentially exact down to zero temperature, whereas it deviates from the real value in the *kagomé* case, as our model does not include the so called order by disorder phenomenon, which is especially important in this system at very low temperatures. It is specially remarkable that the curves that are compared with Monte Carlo data do not contain any fit parameter, or any rescaling factor.

Obviously, this is only a first step in the understanding of the uncommon features present in frustrated geometries. There are still a lot of puzzling questions regarding the nature of the spin glass and spin liquids observed in these systems, which we do not think can be understood in the framework of the constant coupling theory.

In any case, we think that the present model provides an excellent starting point to understand open questions that remain even in the paramagnetic regime, as are, to cite some examples, the effects of small anisotropies, dipole–dipole interactions, or the effects of a small amount of dilution with non–magnetic impurities in the magnetic lattice. Furthermore, the corresponding quantum generalized constant coupling method provides some insight on how the quantum effects manifest themselves in the magnetic properties of these systems Garcia.2000 . It is important to note again, that the applicability of the present technique is not limited to the study of Heisenberg systems, but is useful for the investigation of the physical properties of other systems where geometrical frustration is relevant.

###### Acknowledgements.

Angel García Adeva wants to thank M.A. Cazalilla for very stimulating discussions. This work has been partially supported by the Spanish MEC, under the Subprograma General de Formación de Personal Investigador en el Extranjero.


format: text

text:
# On the Causality and Stability of the Relativistic Diffusion Equation

## I Introduction

The relativistic generalization of the hydrodynamic theory, as developed by Landau and Lifshitz , leads to differential equations of the parabolic type. Consider for instance the diffusion equation for the viscous shear flow. In a frame in which the background equilibrium state is at rest (“comoving frame”), it is

$$[c^2(e+p)_t\eta _x^2]\delta u^i(x,t)=0,$$
(1)
where $`\delta u^i`$ with $`i=y,z`$ is the perturbation in the transverse velocity. Equation (1) is the archetype of a parabolic equation. In non-relativistic physics, it provides an excellent description of a wide range of physical phenomena, as countless experiments have shown. Within the relativistic framework, however, it seems to fail, as two deficiencies become evident: The first, acausality, refers to the fact that Eq. (1) allows for propagation of signals with arbitrarily large velocities. (This should be, and indeed has already been, a worry in the Galilean hydrodynamics, since one has very definite ideas of the velocities of the constituent microscopic particles, which represent an upper limit of the signal velocity in a dilute system of hard-core interaction.) The second defect is an instability found by Hiscock and Lindblom , who showed that Eq. (1) develops a solution that grows exponentially with time in any non-comoving Lorentz frames. It is noteworthy that the growth time scale was found to be microscopically short.

To overcome these deficiencies, extended fluid theories were put forward which start from the hydrodynamic theory but include additional dynamic variables; see e.g. . The resultant larger set of phenomenological coefficients can be chosen such that all the equations are hyperbolic, ensuring causality and stability. The price for this is twofold: A rather more complicated theory and the difficulty of finding a universally accepted set of additional variables (except perhaps in dilute systems). In fact, recently it has been shown by Geroch  and Lindblom  that the complicated dynamical structure which ensures causality is unobservable. The evolution of any physical fluid state according to any causal theory results in energy-momentum tensors and particle currents that are experimentally indistinguishable from the respective hydrodynamic expressions.

In this paper, we take a step back and again focus on the simpler and more universal parabolic equation (1). While acausality is an expected feature of parabolic differential equations, already discussed in the non-relativistic context , instability is not. Being absent in a comoving frame, one is astonished at its appearance in non-comoving frames. Our main purpose thus is to examine the origin and physical relevance of the instability.

The paper is organized as follows: In Sec. II some general mathematical aspects of partial differential equations are reviewed. We especially recall the intuitive meaning of the characteristics. Section III is devoted to the one-dimensional diffusion equation in a comoving frame. We discuss the problem of causality and examine the two types of Cauchy problems that can be formulated with respect to a parabolic equation. In Sec. IV the discussion is generalized to non-comoving Lorentz frames. We especially scrutinize the solution that Hiscock and Lindblom have identified as an instability, and show that it is a result of an ill-posed initial value problem. In Sec. V we briefly examine the general case, in which the spatial dimension of the diffusion equation is greater than one.

## II Preliminaries

Let us first review some general aspects of partial differential equations that can be found in standard textbooks on mathematical physics (eg. ). Consider a linear partial differential equation of second order for the unknown function $`\vartheta (x,t)`$. It can be written, most generally, in the form

$$\left[A_x^2+2B_x_t+C_t^2\right]\vartheta +F(x,t,\vartheta ,_x\vartheta ,_t\vartheta )=0,$$
(2)
where $`A`$, $`B`$ and $`C`$ are given functions of the two independent variables $`x`$ and $`t`$. Depending on the value of $`DB^2AC`$ at a given point, Eq. (2) is referred to be of the elliptic ($`D<0`$), parabolic ($`D=0`$), or hyperbolic ($`D>0`$) type in this point. In the following we shall restrict ourselves to the cases in which $`A`$, $`B`$ and $`C`$ are constants. The type of Eq. (2) then remains unchanged throughout the entire region.

The equation

$$A(_x\phi )^2+2B(_x\phi )(_t\phi )+C(_t\phi )^2=0$$
(3)
is called the equation of characteristics of the partial differential equation (2). Correspondingly, the family of curves,

$$\phi (x,t)=const,$$
(4)
with $`\phi (x,t)`$ satisfying Eq. (3), is called the family of characteristics. We collect the following facts:

(i) An equation of the hyperbolic type has two distinct families of real characteristics, an equation of the parabolic type has only one; an elliptic equation does not have real characteristics.

(ii) The equation of characteristics is invariant with respect to arbitrary transformations of the independent variables, $`\stackrel{~}{x}=\stackrel{~}{x}(x,t)`$, $`\stackrel{~}{t}=\stackrel{~}{t}(x,t)`$. This implies that, if $`\phi (x,t)`$ is a solution of Eq. (3), and if $`\phi (x,t)`$ transforms into $`\stackrel{~}{\phi }(\stackrel{~}{x},\stackrel{~}{t})`$, then $`\stackrel{~}{\phi }(\stackrel{~}{x},\stackrel{~}{t})`$ is a solution of the equation of characteristics accompanying the transformed differential equation.

(iii) The outer real characteristics that pass through a given point $`(x_0,t_0)`$ bound the domain of influence $`\mathrm{\Omega }_0`$ of this point. If we consider the variable $`t`$ as the time and think of the solution $`\vartheta (x,t)`$ as a quantity that varies in $`x`$-space with time $`t`$, then this means that for $`t>t_0`$ the solution in the region outside $`\mathrm{\Omega }_0`$ is not influenced by the initial data given at $`(x_0,t_0)`$. As an example, consider the telegraph equation

$$_t\vartheta \alpha \left(_x^2\frac{1}{\upsilon ^2}_t^2\right)\vartheta =0,$$
(5)
with given constants $`\alpha ,\upsilon >0`$. Since $`D=\alpha ^2/\upsilon ^2>0`$, it belongs to the hyperbolic type. Its equation of characteristics is $`(_x\phi )^2(_t\phi )^2/\upsilon ^2=0`$, which gives two families of characteristics,

$$\phi (x,t)=x\pm \upsilon t=const.$$
(6)
Taking some point $`(x_0,\mathrm{\hspace{0.17em}0})`$, the respective domain of influence is thus bounded by $`x\pm \upsilon t=x_0`$, implying that the effects of the initial data propagate at finite velocity $`\upsilon `$.

For equation (2), the most general Cauchy problem is formulated in the following way: Let $`S`$ be some smooth curve given in the space of the variables $`x`$, $`t`$. With each point $`(x,t)S`$ there is associated some direction $`n`$ not tangent to $`S`$. The problem now consists of finding, in some neighborhood of $`S`$ (either on one or both sides of the curve), a solution $`\vartheta (x,t)`$ which satisfies the prescribed Cauchy conditions

$$\vartheta |_S=\mathrm{\Theta }_0(x,t),\frac{\vartheta }{n}|_S=\mathrm{\Theta }_1(x,t).$$
(7)
It is important to note that the domain in which the unknown solution has to be determined is not specified beforehand. So generally, the initial manifold $`S`$ lies within the domain of definition of the solution.

A problem is said to be well-posed if it has the following properties: The solution (i) exists, (ii) is uniquely determined, and (iii) depends continuously on the assigned data. The last requirement is imposed in connection with the fact that the initial data of physical problems are determined experimentally and so small errors occur. It is thus necessary to be sure that the solution does not depend essentially on the measurement errors of these data.

## III The Diffusion Equation in a Comoving Frame

A well-known example of Eq. (2) is given by the parabolic diffusion equation

$$_t\vartheta \alpha _x^2\vartheta =0,$$
(8)
where $`\alpha >0`$. With $`A=\alpha `$ and $`B=C=0`$ its equation of characteristics (3) takes the form $`(_x\phi )^2=0`$. Hence, the characteristics are given by the one-parameter family of lines

$$\phi (x,t)=t=const.$$
(9)
It is obvious from this that the diffusion equation allows for the propagation of disturbances with infinite velocity. Indeed, the initial value of $`\vartheta `$ at the point $`(x_0,t_0)`$ has influence on the solution $`\vartheta (x,t)`$ in the whole semi-infinite region $`tt_0`$.

This fact, however, does not in practice cause any complications: The superluminal propagation speeds are associated only with variations on microscopically small time and length scales (and with amplitudes of the order of thermodynamic fluctuations). On these scales a macroscopic description loses validity, and the diffusion equation and its solutions break down accordingly. The acausal consequences are therefore precluded by restricting the solutions of Eq. (8) to the hydrodynamic range of validity,

$$\left|\frac{_x\vartheta }{\vartheta \vartheta _0}\right|\frac{1}{\xi },\left|\frac{_t\vartheta }{\vartheta \vartheta _0}\right|\frac{1}{\tau },$$
(10)
with $`\vartheta _0`$ denoting the constant part of $`\vartheta `$; and $`\xi `$, $`\tau `$ the characteristic distance and time between collisions of particles, or elementary excitations. For gases this has been first demonstrated by Weymann .

As an example, consider the solution of Eq. (8) that satisfies the initial condition $`\vartheta |_{t=0}=\vartheta _0+(A/d)e^{x^2/2d^2}`$. It is

$$\vartheta (x,t)=\vartheta _0+(A/\sqrt{d^2+2\alpha t})e^{x^2/2(d^2+2\alpha t)},$$
(11)
which represents a Gaussian distribution with width $`\sqrt{d^2+2\alpha t}`$. As can be seen from the first of Eq. (10), this solution is of physical significance only in the interval

$$|x|d^2/\xi +(2\alpha /\xi )t.$$
(12)
The maximum speed at which measurable information is transmitted is thus of the order of $`\upsilon _{\mathrm{max}}\alpha /\xi `$, which is far less than the speed of light for typical values of the kinetic coefficient $`\alpha `$ and the microscopic length $`\xi `$. (Taking $`\alpha `$ as the heat conductivity, one has for iron at room temperature: $`\alpha 10^5\mathrm{m}^2/\mathrm{s}`$ and $`\xi 10^8`$m. So $`\upsilon _{\mathrm{max}}10^3`$m/s, which is of the order of the speed of sound.)

The most common initial-value problem for Eq. (8) is to find a solution $`\vartheta (x,t)`$ for which the Cauchy conditions (7) are prescribed at the initial manifold $`t=0`$. However, since the line $`t=0`$ is a characteristic here, the Cauchy data cannot be prescribed independently, but must satisfy a compatibility condition . For instance, if the Cauchy conditions are given by $`\vartheta |_{t=0}=\mathrm{\Theta }_0(x)`$ and $`(_t\vartheta )|_{t=0}=\mathrm{\Theta }_1(x)`$, we get from Eq. (8), for $`t=0`$, $`\mathrm{\Theta }_1(x)=\alpha _x^2\mathrm{\Theta }_0(x)`$. The characteristic Cauchy problem for the diffusion equation is therefore posed in the following way: In the region $`t>0`$ find a solution $`\vartheta (x,t)`$ satisfying the initial condition

$$\vartheta |_{t=0}=\mathrm{\Theta }_0(x).$$
(13)
As is generally known , this problem is well-posed for arbitrary, smooth functions $`\mathrm{\Theta }_0(x)`$ that have a well-defined Fourier transform. In the special case of a pure exponential, $`\mathrm{\Theta }_0=Ae^{ikx}`$ ($`k\mathrm{IR}`$), the solution takes the form

$$\vartheta (x,t)=Ae^{ikx+\mathrm{\Gamma }t},\mathrm{\Gamma }(k)=\alpha k^2,$$
(14)
which is unique and stable for $`t>0`$.

Quite another situation arises when the Cauchy data are prescribed on a non-characteristic curve, say eg. $`x+bt=0`$, with some finite $`b\mathrm{IR}`$. The diffusion equation then gives two modes, one bound and the other divergent with $`|x+bt|\mathrm{}`$. A well-known example is the so-called ”sideways problem”: In the half-space $`x>0`$ ($`\mathrm{}<t<\mathrm{}`$) we seek a function $`\vartheta (x,t)`$, which satisfies Eq. (8), and which attains on the non-characteristic line $`x=0`$ the Cauchy conditions

$$\vartheta |_{x=0}=Ae^{i\omega t},(_x\vartheta )|_{x=0}=Be^{i\omega t},$$
(15)
with $`\omega \mathrm{IR}`$. As can be easily seen, this problem is generally solved by

$$\vartheta (x,t)=\left(\delta \vartheta _1e^{(i+1)\mathrm{\Lambda }_1x}+\delta \vartheta _2e^{(i+1)\mathrm{\Lambda }_2x}\right)e^{i\omega t},$$
(16)
$$\mathrm{\Lambda }_{1,2}(\omega )=\pm \sqrt{\omega /2\alpha },$$
(17)
which represents a superposition of two modes

$$\vartheta _{1,2}(x,t)=\delta \vartheta _{1,2}(\omega )e^{i(\pm x\sqrt{\omega /2\alpha }+\omega t)}e^{\pm x\sqrt{\omega /2\alpha }},$$
(18)
with the amplitudes $`\delta \vartheta _{1,2}(\omega )`$ being determined by the two Cauchy conditions (15). One has

$$\delta \vartheta _{1,2}(\omega )=\frac{1}{2}A\pm \frac{1}{4}B(i1)\sqrt{2\alpha /\omega }.$$
(19)
So, without imposing further restrictions, the general solution (16) explodes exponentially as $`x`$ increases. Moreover, it shows discontinuous dependence on the initial data : As $`\omega \mathrm{}`$, Eq. (16) is bounded on the initial line $`x=0`$ but grows like $`\mathrm{exp}(x\sqrt{\omega /2\alpha })`$ for any $`x>0`$. Consequently, the non-characteristic Cauchy problem (15) is not well-posed.

Recalling the intuitive concept of the characteristics (cf. Sec. II), the reason for this becomes obvious. The initial data given at some point $`(0,t_0)`$ affect the value of the solution exactly in those points which lie in the domain of influence $`\mathrm{\Omega }_0=\{(x,t)|x\mathrm{IR},tt_0\}`$. Solving Eqs. (8), (15) thus yields the solution for $`x>0`$ as well as for $`x<0`$. Confer the two modes of Eq. (18). They describe damped waves which carry the initial data given at $`x=0`$ to the left and to the right. Physical intuition suggests, however, that $`\vartheta _1`$ and $`\vartheta _2`$ do exist only for $`x<0`$ and $`x>0`$, respectively. A superposition in the form of Eq. (16) does not make any physical sense.

Now, the sideways problem will not be well-posed unless the behavior at infinity is prescribed. In fact, the physically realistic assumption that $`\vartheta `$ be bounded as $`x\mathrm{}`$  leads to a solution that exists, is unique, and depends continuously on the initial data given at $`x=0`$. The correct formulation of the sideways problem thus is the following: Find the bounded solution of Eq. (8) in the region $`x>0`$ ($`\mathrm{}<t<\mathrm{}`$), satisfying

$$\vartheta |_{x=0}=Ae^{i\omega t},\omega \mathrm{IR}.$$
(20)
Note that the boundedness condition ”replaces” the second Cauchy condition.

Equation (15) is an example of a non-characteristic Cauchy problem with the initial data given on a line that is timelike. For the following, it is important to consider also the case where the Cauchy data are given on a spacelike line. As we shall see, this problem is not well-posed either.

In the half-space $`\{(x,t)|x\mathrm{IR},ct\beta x>0\}`$ we seek the solutions of Eq. (8) satisfying periodic initial data on the non-characteristic spacelike line $`ct\beta x=0`$ ($`\beta \mathrm{IR}`$, $`0<\beta <1`$). Here $`c`$ is the speed of light. With $`\gamma (1\beta ^2)^{1/2}`$, $`x^\mu =(x,ct)`$, $`\eta ^{\mu \nu }=\mathrm{diag}(1,1)`$, and $`n^\mu =\gamma (\beta ,\mathrm{\hspace{0.17em}1})`$, $`e^\mu =\gamma (1,\beta )`$ respectively denoting the timelike and spacelike unit vector normal and parallel to the initial line $`ct\beta x=n^\mu x_\mu =0`$, the general Ansatz

$$\vartheta (x,t)e^{ike^\mu x_\mu +\mathrm{\Gamma }n^\mu x_\mu /c},k\mathrm{IR},\mathrm{\Gamma }\text{ }\text{ }\mathrm{C}$$
(21)
takes the form

$$\vartheta (x,t)e^{ik\gamma (x\beta ct)+\mathrm{\Gamma }\gamma (t\beta x/c)},$$
(22)
with $`\vartheta |_{ct\beta x=0}\mathrm{exp}(ik\gamma ^1x)`$. Inserting Eq. (22) into the diffusion equation (8) yields

$$\gamma \alpha \frac{\beta ^2}{c^2}\mathrm{\Gamma }^2\left(1+2i\gamma \alpha \frac{\beta }{c}k\right)\mathrm{\Gamma }\gamma \alpha k^2+i\beta ck=0.$$
(23)
For $`k0`$ the two roots $`\mathrm{\Gamma }_{1,2}(k)`$ are complex. The real parts are given by

$$\mathrm{\Gamma }_{\mathrm{R1},\mathrm{R2}}=\frac{1}{2}C\pm \sqrt{\frac{1}{8}C^2+\sqrt{\frac{1}{64}C^4+\frac{1}{4}C^3\gamma ^3\alpha k^2}},$$
(24)
where $`C(\gamma \alpha \beta ^2/c^2)^1`$. From this one finds the inequalities

$$\mathrm{\Gamma }_{\mathrm{R1}}+\mathrm{\Gamma }_{\mathrm{R2}}=C>0,$$
(25)
$$\mathrm{\Gamma }_{\mathrm{R1}}\mathrm{\Gamma }_{\mathrm{R2}}=\frac{1}{8}C^2\sqrt{\frac{1}{64}C^4+\frac{1}{4}C^3\gamma ^3\alpha k^2}0,$$
(26)
which imply that exactly one mode grows exponentially with time. Now, as $`k\mathrm{}`$, the general solution diverges like $`\mathrm{exp}(\gamma \sqrt{c^3/2\alpha \beta }[t\beta x/c]\sqrt{k})`$ for any $`t>\beta x/c`$, while the initial values are bounded on $`t\beta x/c=0`$. Hence small changes in the initial data would cause considerable changes in the solution, indicating that the problem is not well-posed.

## IV The Diffusion Equation in a General Lorentz Frame

Next we turn our attention to the diffusion equation in an inertial frame $`\stackrel{~}{K}`$, in which the medium moves with constant velocity $`v`$ in the negative $`x`$-direction. Employing the Lorentz transformation rules $`_t=\gamma _{\stackrel{~}{t}}\gamma v_{\stackrel{~}{x}}`$ and $`_x=\gamma _{\stackrel{~}{x}}\gamma (v/c^2)_{\stackrel{~}{t}}`$, we get from Eq. (8) the boosted diffusion equation

$$\gamma (_{\stackrel{~}{t}}v_{\stackrel{~}{x}})\vartheta \alpha \gamma ^2\left(_{\stackrel{~}{x}}^22\frac{v}{c^2}_{\stackrel{~}{x}}_{\stackrel{~}{t}}+\frac{v^2}{c^4}_{\stackrel{~}{t}}^2\right)\vartheta =0.$$
(27)
It is straightforward to verify that this equation still belongs to the parabolic type. In a covariant language, it is rewritten as

$$u^\mu _\mu \vartheta \alpha \mathrm{\Delta }^{\mu \nu }_\mu _\nu \vartheta =0,$$
(28)
where $`\mathrm{\Delta }^{\mu \nu }=\eta ^{\mu \nu }+c^2u^\mu u^\nu `$, $`u^\mu =\gamma (v,c)`$, $`_\mu =(_{\stackrel{~}{x}},_{c\stackrel{~}{t}})`$, and $`\eta ^{\mu \nu }=\mathrm{diag}(1,1)`$. According to Eq. (3), the equation of characteristics becomes

$$\gamma ^2\left(_{\stackrel{~}{x}}\stackrel{~}{\phi }\frac{v}{c^2}_{\stackrel{~}{t}}\stackrel{~}{\phi }\right)^2=0,$$
(29)
so that the family of characteristics comes out as the general solution of the ordinary differential equation $`\gamma vd\stackrel{~}{x}+\gamma c^2d\stackrel{~}{t}=0`$, or equivalently,

$$u_\mu dx^\mu =0.$$
(30)
One finds

$$\stackrel{~}{\phi }(\stackrel{~}{x},\stackrel{~}{t})=\gamma \stackrel{~}{t}+\gamma \frac{v}{c^2}\stackrel{~}{x}=t=const,$$
(31)
with $`t`$ being the proper time measured in a comoving frame $`K`$. Recalling Eq. (9), we explicitly see that the characteristics are invariant under Lorentz boosts; cf. Sec. II.

With Eq. (31), the relativistic equivalent to the characteristic Cauchy problem (13) is the following: Find a solution to Eq. (27) which satisfies prescribed values at the initial manifold $`\stackrel{~}{t}+(v/c^2)\stackrel{~}{x}=0`$. Clearly, this problem is well-posed only for $`\stackrel{~}{t}+(v/c^2)\stackrel{~}{x}>0`$. Its exponential solution can easily be found from Eq. (14) by making use of the Lorentz transformation, $`t=\gamma \stackrel{~}{t}+\gamma (v/c^2)\stackrel{~}{x}`$ and $`x=\gamma \stackrel{~}{x}+\gamma v\stackrel{~}{t}`$.

Now, the solutions of Eq. (27) that Hiscock and Lindblom  examine satisfy periodic initial data $`e^{i\stackrel{~}{k}\stackrel{~}{x}}`$ ($`\stackrel{~}{k}\mathrm{IR}`$) on the non-characteristic line $`\stackrel{~}{t}=0`$. Taking

$$\vartheta (\stackrel{~}{x},\stackrel{~}{t})e^{i\stackrel{~}{k}\stackrel{~}{x}+\stackrel{~}{\mathrm{\Gamma }}\stackrel{~}{t}},$$
(32)
they find the dispersion relation

$$\gamma \alpha \frac{v^2}{c^4}\stackrel{~}{\mathrm{\Gamma }}^2\left(1+2i\gamma \alpha \frac{v}{c^2}\stackrel{~}{k}\right)\stackrel{~}{\mathrm{\Gamma }}\gamma \alpha \stackrel{~}{k}^2+iv\stackrel{~}{k}=0,$$
(33)
see Eq. (67) of Ref. . For $`\stackrel{~}{k}0`$ the two roots $`\stackrel{~}{\mathrm{\Gamma }}_{1,2}(\stackrel{~}{k})`$ are complex. The real parts satisfy the conditions

$$\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{R1}}+\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{R2}}=\frac{c^4}{\gamma \alpha v^2}>0,$$
(34)
$$\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{R1}}\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{R2}}=\left(\mathrm{\Gamma }_{\mathrm{I}\mathrm{\hspace{0.17em}1}}\frac{c^2\stackrel{~}{k}}{v}\right)^20,$$
(35)
which imply that one of the two modes diverges as $`\stackrel{~}{t}\mathrm{}`$. ¿From this the authors of Ref.  conclude that the Landau-Lifshitz theory of relativistic hydrodynamics is unstable (in the sense that small spatially bounded departures from equilibrium will diverge with time) and hence must be abandoned as a physically realistic theory.

We disagree with this conclusion. The reason is that the above non-characteristic Cauchy problem is not well-posed, in complete analogy to its non-relativistic equivalent, Eqs. (8), (22). \[Note that Eq. (23) equals Eq. (33) if $`\beta =v/c`$.\] Its solution does not depend continuously on the initial data and therefore does not meet one of the three mathematical requirements to be posed with regard to realistic physical problems (cf. Sec. II). So it is the type of the initial-value problem here which is physically unacceptable, and not the instability of the resulting solution.

Nevertheless, it is illuminating to interpret the Hiscock-Lindblom solution in terms of wave propagation (while ignoring the fact that the domain of definition of the above Cauchy problem is restricted to the half space $`\{(\stackrel{~}{x},\stackrel{~}{t})|\stackrel{~}{x}\mathrm{IR},\stackrel{~}{t}0\}`$). The general solution consists of two damped waves traveling to the regions $`\stackrel{~}{t}>0`$ and $`\stackrel{~}{t}<0`$, respectively. Referring again to the discussion of the sideways problem, Sec. III, this is a consequence to be expected. Since the initial manifold crosses the characteristics (31), the information (given at $`\stackrel{~}{t}=0`$) propagates (along the characteristic lines $`\stackrel{~}{t}+[v/c^2]\stackrel{~}{x}=const`$) in both directions. One thus gets two modes, each transporting the effects of the initial data in the respective region.

It is important to note that the appearance of a mode running backwards in time in non-comoving frames is directly related to the infinite signal speeds accompanying the parabolic equation. In fact, it is well-known  that “superluminal acausality” in a Lorentz frame $`K`$ causes “chronological acausality” in any other Lorentz frames $`\stackrel{~}{K}`$. For example, consider a process in a frame $`K`$ whereby an event $`𝒫_1=(x,t)`$ causes another event $`𝒫_2=(x+\mathrm{\Delta }x,t+\mathrm{\Delta }t)`$ at superluminal velocity $`U>c`$. Let the time difference be $`\mathrm{\Delta }t>0`$ so that $`\mathrm{\Delta }x=U\mathrm{\Delta }t>0`$. Then, in a frame $`\stackrel{~}{K}`$, moving with velocity $`v`$ relative to $`K`$, we have $`\mathrm{\Delta }\stackrel{~}{t}=\gamma \mathrm{\Delta }t\gamma (v^2/c)\mathrm{\Delta }x=\gamma \mathrm{\Delta }t(1vU/c^2)`$. If $`c^2/U<v<c`$, this yields $`\mathrm{\Delta }\stackrel{~}{t}<0`$. This means that in $`\stackrel{~}{K}`$ the signal goes backward in time, or equivalently, that the response $`𝒫_2`$ precedes the stimulus $`𝒫_1`$.

Now, as we have seen in Sec. III, the superluminal acausality in $`K`$ does not lead to any physical consequences; it is automatically precluded by a restriction to the hydrodynamic range of validity. Accordingly, we see from Eqs. (34) and (35) that the damping rate of the mode which propagates in the past is bounded below by $`\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{R}+}c^4/\gamma \alpha v^2`$. With $`\tau `$ denoting the characteristic microscopic time in $`K`$ (referred to as the collision time in dilute systems), and $`\stackrel{~}{\tau }=\gamma \tau `$ the respective time in $`\stackrel{~}{K}`$, we thus have

$$\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{R}+}\stackrel{~}{\tau }=\frac{c^4\tau }{\alpha v^2}1.$$
(36)
(For instance, taking $`\alpha `$ as the heat conductivity of a piece of iron, one has $`\alpha 10^5\mathrm{m}^2/\mathrm{s}`$ and $`\tau 10^{12}`$s, so that $`\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{R}+}\stackrel{~}{\tau }10^{10}`$.) This shows that the mode running backwards in time $`\stackrel{~}{t}`$ decays on a time-scale that is much shorter than the microscopic one, and hence is far outside the hydrodynamic regime.

Let us now come back to the non-characteristic Cauchy problem. Since the domain of definition of the solution is $`\{(\stackrel{~}{x},\stackrel{~}{t})|\stackrel{~}{x}\mathrm{IR},\stackrel{~}{t}0\}`$, the mode reflecting the chronological acausality is defined only for positive times $`\stackrel{~}{t}`$. As a consequence, the superposed solution does not depend continuously on the initial data, and the problem is ill-posed. With other words, prescribing Cauchy conditions on a non-characteristic hypersurface in a 2-dimensional flat spacetime does not lead to a realistic physical problem. Such a problem might be obtained only if one prescribes the value of the solution at infinity (cf. the sideways problem in Sec. III), or if the Cauchy data are prescribed on one of the characteristics.

Now, what are the physical conclusions after the mathematical ones have been drawn? Given the fact that the standard Cauchy problem is ill-posed in the Landau-Lifshitz theory, that one cannot arbitrarily prescribe initial data of physical quantities on a boosted time slice, are we to conclude that the theory is defunct – and to be replaced by one in which the Cauchy problem is well-posed? Is the capability to accommodate arbitrary initial data on boosted time slices a sine qua none for a healthy, physical theory? We believe the answer is no, for two different reasons, although we concede that an affirmative answer may also be upheld. Our reasons are

* All differential equations belonging the the hydrodynamic theory and accounting for coarse-grained, irreversible physics – including the diffusion equation – possess a preferred inertial system, in which the material is at rest. This is not true with respect to microscopic theories for which the Cauchy problem is known to be well-posed in any inertial frame, consider for instance the vacuum Maxwell equations.
* More importantly, in our simple example considered above, we actually know what the initial data on a boosted time slice $`\stackrel{~}{t}=0`$ correspond to in the rest frame, namely to an ill-posed, sideway problem arising from data on a non-characteristic spacelike line $`ct\gamma (v/c)x=0`$. This connection is completely general, and especially independent of the Landau-Lifshitz theory. Therefore, an insistence on the well-posedness of Cauchy problems for boosted systems amounts to requiring a new type of differential equation for the rest frame, ones for which the sideway problem is well-posed. Now, there are few equations in physics which are better confirmed than the rest frame diffusion equations as we know them: for temperature, velocity and concentration. And it seems highly unlikely that one can change something as basic as the ill-posedness of the sideway problem without destroying the agreement with the experiments.

Resigning to the fact that the standard Cauchy problem on a boosted time slice may be ill-posed, then clearly, there is no reason whatever to abandon the Landau-Lifshitz theory.

## V Generalization to cases of more than one spatial dimension

Let us finally examine the general case in which the spatial dimension of the diffusion equation is greater than one. Consider the covariant equation (28), with $`\mu `$ running now from 1 to 4. As we have seen in Sec. IV, its characteristics are determined by the ordinary differential equation

$$u_\mu dx^\mu =0,$$
(37)
which assigns to each point $`x^\mu `$ an infinitesimal characteristic 3-flat normal to the direction of the time-like velocity $`u^\mu `$ at that point. The crucial point now is that generally these infinitesimal 3-flats do not integrate to 3-surfaces $`\phi (x^\mu )=const`$.

To see this, let us suppose that a general integral $`\phi (x^\mu )=const`$ exists. Then, we have

$$d\phi =(_\mu \phi )dx^\mu =0,$$
(38)
and hence

$$_\mu \phi =\lambda u_\mu ,$$
(39)
with some $`\lambda =\lambda (x^\mu )`$. Since $`_\mu _\nu \phi =_\nu _\mu \phi `$, this yields

$$_\mu u_\nu _\nu u_\mu =\lambda ^1(u_\mu _\nu \lambda u_\nu _\mu \lambda ),$$
(40)
from which $`\epsilon ^{\mu \nu \alpha \beta }u_\nu _\alpha u_\beta =0`$ or

$`u_\alpha (_\nu u_\mu _\mu u_\nu )`$ (41)
$`+u_\mu (_\alpha u_\nu _\nu u_\alpha )+u_\nu (_\mu u_\alpha _\alpha u_\mu )=0`$ (42)
results. Equation (42) is a necessary condition for integrability. It can be proved (see eg. ) that it is also sufficient, ie., if it is satisfied, a general integral exists. \[Note that Eq. (42) is identically satisfied for $`\alpha ,\mu ,\nu =1,\mathrm{\hspace{0.17em}2}`$. Thus, in the case of one spatial dimension, characteristics can always be found.\]

Multiplying the condition of integrability with $`u^\alpha `$, one finds

$$\mathrm{\Delta }^{\mu \alpha }\mathrm{\Delta }^{\nu \beta }(_\alpha u_\beta _\beta u_\alpha )=0.$$
(43)
So, characteristic hypersurfaces exist only if the background equilibrium state is non-rotating. This fact has led various authors to believe that generally the parabolic diffusion equation lacks an initial-value formulation, and hence turns out not to be viable. This, however, would be true only if the non-characteristic Cauchy problem is ill-posed, just like it is in the case of one spatial dimension. Unfortunately, very little seems to be known about the non-characteristic Cauchy problem in spacetimes with more than one spatial dimension. Further work should give clarity here.


text:
# Neutron Star Structure and the Neutron Radius of 208Pb

\[

## Abstract

We study relationships between the neutron-rich skin of a heavy nucleus and the properties of neutron-star crusts. Relativistic effective field theories with a thicker neutron skin in 208Pb have a larger electron fraction and a lower liquid-to-solid transition density for neutron-rich matter. These properties are determined by the density dependence of the symmetry energy which we vary by adding nonlinear couplings between isoscalar and isovector mesons. An accurate measurement of the neutron radius in 208Pb—via parity violating electron scattering—may have important implications for the structure of the crust of neutron stars.

\]

It is an extrapolation of 18 orders of magnitude from the neutron radius of a heavy nucleus—such as 208Pb with $`R_n5.5`$ fm—to the approximately 10 km radius of a neutron star. Yet both radii depend on our incomplete knowledge of the equation of state of neutron-rich matter. Therefore, an accurate measurement of the neutron radius in 208Pb may have important implications for some neutron star properties.

Heavy nuclei are expected to have a neutron-rich skin. This important feature of nuclear structure arises because of the large neutron excess and because the Coulomb barrier reduces the proton density at the surface. The thickness of the neutron skin depends on the pressure of neutron-rich matter: the greater the pressure, the thicker the skin as neutrons are pushed out against surface tension. The same pressure supports a neutron star against gravity . Thus models with thicker neutron skins often produce neutron stars with larger radii.

Neutron stars are expected to have a solid crust of nonuniform neutron-rich matter above a liquid mantle. The phase transition from solid to liquid depends on the properties of neutron-rich matter. Indeed, a high pressure implies a rapid rise of the energy with density making it energetically unfavorable to separate uniform matter into regions of high and low densities. Thus a high pressure typically implies a low transition density from a solid crust to a liquid mantle. This suggests an inverse relationship: the thicker the neutron-rich skin of a heavy nucleus, the thinner the solid crust of a neutron star.

In this letter we study possible “data-to-data” relations between the neutron-rich skin of a heavy nucleus and the crust of a neutron star. These relations may impact neutron star observables. Indeed, properties of the crust are important for models of glitches in the rotational period of pulsars , for the shape and gravitational radiation of non-spherical rotating stars  and for neutron-star cooling . Note that the skin of a heavy nucleus and the crust of a neutron star are composed of the same material: neutron-rich matter at similar densities.

The Parity Radius Experiment (PREX) at the Jefferson Laboratory aims to measure the neutron radius in 208Pb via parity violating electron scattering . Parity violation is sensitive to the neutron density because the $`Z^0`$ boson couples primarily to neutrons. The result of this purely electroweak experiment could be both accurate and model independent. In contrast, all previous measurements of bulk neutron densities used hadronic probes that suffer from controversial uncertainties in the reaction mechanism (see for example Ref. ). PREX should provide a unique observational constraint on the thickness of the neutron skin in a heavy nucleus. In this letter we explore some of the implications of this measurement on the structure of neutron stars.

Microscopic calculations of the energy of neutron matter constrain both the neutron skin in 208Pb and the crust of a neutron star; see for example Ref. . However, these calculations of infinite neutron matter are not directly tested by observable properties of finite nuclei such as their charge densities or binding energies. Moreover, nonrelativistic calculations of symmetric nuclear matter have not succeeded in predicting the saturation density. It thus becomes necessary to fit some properties of a three-body force in order to reproduce nuclear saturation. Indeed, the properties of $`A=8`$ pure neutron drops calculated in Ref.  may depend on the three-nucleon force used. Thus we feel that it is important to distinguish direct finite-nucleus measurements—such as PREX—from theoretical neutron-matter “observables” based solely on calculations. Indeed, PREX may provide an important test of these calculations .

We start with a relativistic effective field theory  that provides a simple description of finite nuclei and a Lorentz covariant extrapolation for the equation of state of dense neutron-rich matter. The theory has an isoscalar-scalar $`\varphi `$ (sigma) meson field and three vector fields: an isoscalar $`V`$ (omega), an isovector $`b`$ (rho), and the photon $`A`$. We now supplement the Lagrangian with new nonlinear sigma-rho and omega-rho couplings. These couplings allow us to change the density dependence of the symmetry energy which changes both the thickness of the neutron skin in 208Pb and the neutron-star crust.

The interacting Lagrangian density is given by ,

$`_{\mathrm{int}}=\overline{\psi }\left[g_\mathrm{s}\varphi \left(g_\mathrm{v}V_\mu +{\displaystyle \frac{g_\rho }{2}}𝝉𝐛_\mu +{\displaystyle \frac{e}{2}}(1+\tau _3)A_\mu \right)\gamma ^\mu \right]\psi `$ (1)
$`{\displaystyle \frac{\kappa }{3!}}(g_\mathrm{s}\varphi )^3{\displaystyle \frac{\lambda }{4!}}(g_\mathrm{s}\varphi )^4+{\displaystyle \frac{\zeta }{4!}}g_\mathrm{v}^4(V_\mu V^\mu )^2+{\displaystyle \frac{\xi }{4!}}g_\rho ^4(𝐛_\mu 𝐛^\mu )^2`$ (2)
$`+g_\rho ^2𝐛_\mu 𝐛^\mu \left[\mathrm{\Lambda }_4g_\mathrm{s}^2\varphi ^2+\mathrm{\Lambda }_\mathrm{v}g_\mathrm{v}^2V_\mu V^\mu \right].`$ (3)
We consider a number of different parameter sets. First, we note that the nonlinear rho coupling $`\xi `$ will modify the density dependence of the rho mean field and this could change the neutron-skin thickness. However, unless $`\xi `$ is made very large, this term was found to have a small effect . Thus for simplicity, we set $`\xi 0`$ in all our parameter sets. Second, note that we could have added a cubic sigma-rho interaction of the form: $`_3=M\mathrm{\Lambda }_3(g_\mathrm{s}\varphi )(g_\rho ^2𝐛_\mu 𝐛^\mu )`$. A nonzero $`\mathrm{\Lambda }_3`$ does change the thickness of the neutron skin in 208Pb—but at the expense of a change in the proton density. Therefore, we set $`\mathrm{\Lambda }_30`$ and focus exclusively on $`\mathrm{\Lambda }_4`$ and $`\mathrm{\Lambda }_\mathrm{v}`$.

We start with the original NL3 parameter set of Lalazissis, König, and Ring . (Note that a small adjustment of the $`NN\rho `$ coupling constant was needed to fit the symmetry energy of nuclear matter at a Fermi momentum of $`k_F=1.15\mathrm{fm}^1`$; see text below). The NL3 set has $`\zeta =\mathrm{\Lambda }_4=\mathrm{\Lambda }_\mathrm{v}=0`$ and provides a good fit to the ground-state properties of many nuclei. In this model symmetric nuclear matter saturates at $`k_F=1.30`$ fm-1 with a binding energy per nucleon of $`E/A=16.25`$ MeV and an incompressibility of $`K=271`$ MeV. All other parameter sets considered here have been fixed to the same saturation properties.

We now add the new nonlinear couplings $`\mathrm{\Lambda }_4`$ and/or $`\mathrm{\Lambda }_\mathrm{v}`$ between the isoscalar and the isovector mesons. At the same time we adjust the strength of the $`NN\rho `$ coupling constant ($`g_\rho `$) from its NL3 value to maintain the symmetry energy of nuclear matter unchanged (see text below). Note that neither $`\mathrm{\Lambda }_4`$ nor $`\mathrm{\Lambda }_\mathrm{v}`$ affect the properties of symmetric nuclear matter since $`𝐛_\mu 0`$. Hence, the saturation properties remain unchanged. Our goal is to change the neutron density and the neutron-skin thickness in 208Pb while making very small changes to the proton density which is well constrained by the measured charge density . The two new couplings ($`\mathrm{\Lambda }_4`$ and $`\mathrm{\Lambda }_\mathrm{v}`$) change the skin thickness in 208Pb by similar amounts. Yet they have different high-density limits. For $`\mathrm{\Lambda }_\mathrm{v}=0`$ the symmetry energy is proportional to the baryon density $`\rho `$ in the limit of very high density, while it only grows like $`\rho ^{1/3}`$ for nonzero $`\mathrm{\Lambda }_\mathrm{v}`$. This can produce noticeable differences in neutron star radii, as we show below.

Let us start with $`\mathrm{\Lambda }_4=0`$. For a given omega-rho coupling $`\mathrm{\Lambda }_\mathrm{v}`$ we readjust only the $`NN\rho `$ coupling constant $`g_\rho `$ in order to keep an average symmetry energy fixed. The symmetry energy at saturation density is not well constrained by the binding energy of nuclei. However, some average of the symmetry energy at full density and the surface energy is constrained by binding energies. As a simple approximation we evaluate the symmetry energy not at full density $`k_F1.30`$ fm-1 but rather at an average density corresponding to $`k_F=1.15`$ fm-1. Thus, all our parameter sets have a symmetry energy of $`23.50`$ MeV at $`k_F=1.15`$ fm-1. Note that the original NL3 parameter set predicts a symmetry energy of 37.4 MeV at full saturation density and close to 23.50 MeV at $`k_F=1.15`$ fm-1 . This simple procedure produces a nearly constant binding energy per nucleon for 208Pb as $`\mathrm{\Lambda }_\mathrm{v}`$ is changed, as can be seen in Table I. Moreover, Table I shows that increasing $`\mathrm{\Lambda }_\mathrm{v}`$ reduces the neutron-skin thickness significantly—while maintaining the proton radius nearly constant. In the following we plot our results for a range of $`\mathrm{\Lambda }_\mathrm{v}`$ values for which the proton radius is within 0.01 fm of its $`\mathrm{\Lambda }_\mathrm{v}=0`$ value.

To study the solid crust of a neutron star we make a simple random-phase-approximation (RPA) calculation of the transition density below which uniform neutron-rich matter becomes unstable against small amplitude density fluctuations. This provides a lower bound to the true transition density . We start with the longitudinal dielectric function $`ϵ_L`$, as defined in Eq. (68) of Ref. , evaluated at an energy transfer $`q_0=0`$ and at an arbitrary momentum transfer $`q`$. That is,

$$ϵ_L(q_0=0,q)=det\left(1D_L\mathrm{\Pi }_L\right).$$
(4)
Here $`\mathrm{\Pi }_L`$ is a longitudinal polarization matrix describing particle-hole excitations of a uniform system of protons, neutrons, and electrons in beta equilibrium, as given in Eq. (56) of Ref. . The matrix $`D_L`$, describing meson and photon propagation, follows from Eq. (57) of Ref. —but includes additional terms to account for the nonlinear nature of the meson self-couplings . We estimate the transition density $`\rho _c`$ by computing the largest density at which $`ϵ_L(0,q)<0`$ for any given $`q`$.

FIG. 1. Estimate of the transition density from nonuniform to uniform neutron-rich matter versus neutron-minus-proton radius in 208Pb. The curves are for the four parameter sets described in the text.

In Fig. 1 we display the transition density for various parameter sets (see Table II) as a function of the predicted difference in the root-mean-square neutron and proton radii $`R_nR_p`$ in 208Pb. The curves are parameterized by different values of $`\mathrm{\Lambda }_\mathrm{v}`$, as shown in Table I. The NL3 parameter set saturates nuclear matter with a relatively small value of the nucleon effective mass: $`M^{}Mg_s\varphi =0.59M`$. The parameter set S271 saturates nuclear matter as NL3 but with $`M^{}=0.70M`$. This set also has $`\zeta =0`$. The two remaining curves in the figure are for parameter sets having $`\zeta =0.06`$ and both saturate nuclear matter with $`M^{}=0.80M`$. (Set Z271v has a nonzero $`\mathrm{\Lambda }_\mathrm{v}`$ while set Z2714 uses a nonzero $`\mathrm{\Lambda }_4`$). Note that the scalar mass $`m_\mathrm{s}`$ for parameter sets S271, Z271v, and Z2714 is adjusted to reproduce the proton radius in 208Pb as computed with NL3. Figure 1 displays a clear inverse correlation between the transition density and the neutron-skin thickness $`R_nR_p`$. The transition density expressed in fm-3 is about,

$$\rho _c0.160.39(R_nR_p),$$
(5)
with the skin thickness expressed in fm. Moreover, this correlation seems to be insensitive to $`M^{}`$ or to using $`\mathrm{\Lambda }_4`$ or $`\mathrm{\Lambda }_\mathrm{v}`$ to change $`R_nR_p`$. These results suggest that a measurement of the neutron radius in 208Pb will provide considerable information on the transition density.

Note that Fig. 1 only shows our results. Yet all other calculations that we are aware of also give consistent results. For example, the nonrelativistic microscopic equation of state of Friedman and Pandharipande has a transition density of $`\rho _c=0.096`$ fm-3 according to Lorenz et al. . For this equation of state Brown finds $`R_nR_p=0.16\pm 0.02`$ fm . These numbers are in excellent agreement with Eq. (3).

Brown also finds a linear relation between $`R_nR_p`$ and the derivative of the energy of neutron matter versus density $`dE/d\rho `$ evaluated at $`\rho =0.1`$ fm-3 . He considers a large variety of nonrelativistic Skyrme interactions. Our results for $`dE/d\rho `$ versus $`R_nR_p`$ are completely consistent. We expect these common $`dE/d\rho `$ values to give similar $`\rho _c`$ values consistent with Eq. (3) for these Skyrme interactions. This is because $`dE/d\rho `$ is related to the pressure while $`\rho _c`$ depends on the density dependence of the pressure.

Finally, for the relativistic interaction TM1 of Sugahara and Toki , we calculate from Eq. (2) $`\rho _c0.059`$ fm-3 and $`R_nR_p=0.27`$ fm. Again, the numbers are again in good agreement with Eq. (3).

In Fig. 2 we show the electron fraction per baryon $`Y_e`$ versus density for uniform neutron-rich matter in beta equilibrium. We include results only for the S271 parameter set as all other sets yield similar results. The different curves are for different values of $`\mathrm{\Lambda }_\mathrm{v}`$ which predict the indicated $`R_nR_p`$ values. The curves start near the transition densities displayed in Fig. 1. The electron fraction $`Y_e`$ is determined by the symmetry energy while $`R_nR_p`$ is sensitive to the density dependence of the symmetry energy. Therefore a measurement of $`R_nR_p`$ constrains the growth of $`Y_e`$ with density. If $`R_nR_p`$ is greater than about 0.24 fm, $`Y_e`$ becomes large enough to allow the direct URCA process  to cool down a 1.4 solar mass neutron star.

FIG. 2. Electron fraction $`Y_e`$ versus baryon density for uniform neutron-rich matter in beta equilibrium using the S271 parameter set (other sets yield similar results). The curves are for different values of $`\mathrm{\Lambda }_\mathrm{v}`$ that predict the indicated values of $`R_nR_p`$ for 208Pb, in fm.

We now consider the radius $`R`$ of a 1.4 solar mass neutron star. For a given parameter set, $`R`$ increases with $`R_nR_p`$. However, as one changes the parameter set to increase $`M^{}`$ or $`\zeta `$ the equation of state becomes softer at high density. As a result $`R`$ decreases for fixed $`R_nR_p`$. Also, the high density equation of state is softer with $`\mathrm{\Lambda }_\mathrm{v}`$ than with $`\mathrm{\Lambda }_4`$ so Z271v gives slightly smaller stars than parameter set Z2714. We conclude that $`R`$ is not uniquely constrained by a measurement of the neutron-skin thickness because $`R_nR_p`$ only depends on the equation of state at normal and lower densities while $`R`$ is also sensitive to the equation of state at higher densities. Yet one may be able to combine separate measurements of $`R_nR_p`$ and $`R`$ to obtain considerable information about the equation of state at low and high densities. For example, if $`R_nR_p`$ is relatively large while $`R`$ is small this could suggest a phase transition. A large $`R_nR_p`$ implies that the low density equation of state is stiff while a small $`R`$ suggests that the high density equation of state is soft. The transition from stiff to soft could be accompanied by a phase transition.

FIG. 3. Radius of a 1.4 solar mass neutron star versus neutron-minus-proton radius in 208Pb for the four parameter sets described in the text.

In conclusion: 1) It is possible to fit nuclear observables—such as charge densities, binding energies, and single particle spectra—with effective field theories that predict a range of neutron-skin thicknesses. This can be done by adding nonlinear couplings between isoscalar and isovector meson fields or, in general, by adding interactions that modify the density dependence of the symmetry energy. We conclude that the neutron-skin thickness is not tightly constrained by these observables. Yet a measurement of the skin thickness will constrain the density dependence of the symmetry energy.

2) The density dependence of the symmetry energy is adjustable in our relativistic effective field models while still reproducing nuclear-matter properties and other ground-state observables. Indeed, our models can provide a Lorentz-covariant extrapolation for the high density equation of state with a symmetry energy that rises slower with density relative to earlier relativistic mean-field models.

3) The electron fraction $`Y_e`$ of neutron-rich matter in beta equilibrium is correlated with the neutron-skin thickness in 208Pb. The thicker the neutron skin the faster $`Y_e`$ rises with density. In our models a neutron-skin thickness of the order of 0.24 fm or larger suggests that $`Y_e`$ will become large enough to allow a direct URCA process to cool down a $`1.4M_{}`$ neutron star.

4) We have found an inverse correlation between the neutron-skin thickness and the density of a phase transition from nonuniform to uniform neutron-rich matter. In our models the transition density obeys the approximate relation: $`\rho _c0.160.39(R_nR_p)`$, with $`R_nR_p`$ in fm and $`\rho _c`$ in fm-3. This suggests that a neutron skin measurement in 208Pb can provide important information on the thickness and other properties of the crust of a neutron star.

5) Microscopic calculations of the energy of neutron matter constrain the density dependence of the symmetry energy and hence the neutron-skin thickness in 208Pb. Therefore a neutron skin measurement may provide an important observational check on such calculations that is not provided by other observables. Moreover, a neutron-skin measurement may constrain three-body forces in neutron-rich matter.

We acknowledge useful discussions with A. Brown, R. Furnstahl, J. Lattimer, S. Reddy and B. Serot. We thank A. Brown for sharing results prior to publication. This work was supported in part by DOE grants DE-FG02-87ER40365 and DE-FG05-92ER40750.