text:
# Mechanical Properties of a Model of Attractive Colloidal Solutions

## 1 Introduction

In recent years the idea that colloidal particles can form glassy structures has been established in a number of very interesting experimental and theoretical works. Most attention has been focused on colloidal particle systems that are dominated by repulsive interactions, for which, at high packing fraction values, the glass represents an alternative packing to the usually more favourable crystal structure. These systems are well represented by a simple hard sphere model . Where repulsive interactions dominate, the loss of ergodicity is lost due to blocking of the movements of particles by the quite dense surrounding cages formed by their nearest neighbours. In colloidal systems, Mode Coupling Theory  has played a leading role, interpreting and rationalizing some of the observations , and achieving quite acceptable numerical agreement in comparison to experiments .

Recently a new type of glass has arisen as the focus of interest . This has been called the attractive glass. By attractive glass we imply that the loss of ergodicity is driven largely by strong short-ranged attractive interactions, in other words the ‘stickiness’ of the particles to each other eventually dominates the thermal motions, and the system freezes. Thus, close packing is no longer necessary for a glass to be stable, and it transpires that such glasses can form at densities much lower than close packing.

Once having established the distinct nature of the attractive glass phenomenon, it is unsurprising that new phenomena should be associated with the system. For instance, MCT predicts that, for well defined values of the temperature or packing fraction, the repulsive and attractive glasses differ only by their dynamical properties (as opposed to structural differences). Many questions remained unanswered in this arena, and more research, using different approaches must be developed before one has great confidence in the conclusions drawn so far. Nevertheless, the basic results emerging from MCT seem reasonable, and are reproducible, using a variety of different input static structure factors .

The situation in experimental studies is even less clear. It may be that colloidal glasses driven by attractive interactions have been observed and studied for some time, without clear recognition of their distinctive nature. Certainly some of the reported kinetically arrested states in colloidal systems, usually identified as ‘gels’, appear to arise at densities considerably lower than close packing . Given that such observations appear to be associated with quite strong short-ranged attractions driven by depletion forces, or grafted chains onto larger particles, it would seem likely that indeed such systems are examples of those we currently discuss. On the other hand, as yet, there appear to be no clear reports of colloidal glass-glass transition in the experimental literature, but the typical logarithmic decay of density correlations, that MCT predicts to happen close to an $`A_3`$ singularity  has been observed in the past, for the glass transition of some polymeric systems . More recently a logarithmic decay has been reported in a micellar system with short-ranged attractive interactions, and this may be related to a proximate glass-glass transitions as these authors point out . The clear recognition of two different glasses would be one clear and unambiguous signal that a distinctive attractive glassy state has been observed. Since this is expected to be an issue of considerable interest to experimentalists, it will be important to discuss those measurements that could differentiate the glasses, and this is one aim of the present work.

In this paper it is our intention to discuss the mechanical properties of the glasses formed when repulsive and attractive interactions are present, and where, under some conditions, the latter may become dominant. We shall work within the framework of the ideal MCT, and present results for the linear elastic shear modulus at zero frequency , $`G^{}(\omega =0)`$, and for the longitudinal stress modulus in the hydrodynamic limit, corresponding to zero frequency and long wavelengths, $`m_0`$, to illustrate our point of view. The study of the elastic shear modulus has previously been addressed in , but there a different potential is used.

Our reasoning for presenting these results has two underlying strands of rationale. Firstly, we have outlined the interesting and novel features we believe to be present for colloidal systems where there are strong short-ranged potentials. However, we note that experiments in dense colloidal systems are far from simple, and the number of techniques that can be reliably applied to examine these questions are much smaller than in dilute solutions due to multiple scattering and other limitations implied at high concentrations. Furthermore, when one works on the boundary between liquids and amorphous solids there are further restrictions in the options for experimental techniques that may be applied. In fact, the mechanical responses of an amorphous materials are amongst the most simple and reliable methods of characterisation. Even the bulk shear and longitudinal moduli are quite helpful in forming an initial assessment of the material itself, and of course they are themselves amongst the simplest parameters characterising the transition from a liquid to a glass, and as we shall see, even between different glasses. Combined with this, it is clear that most practical applications of amorphous materials composed from colloidal particles will involve a strong interest in mechanical properties of the glass. Having said all this, the reader should be aware that many attractive glasses may be quite fragile, and some thought will have to be applied to finding systems and methods where measurements can be made. In any case, these comments comprise the first strand of reasoning for our interest in the mechanical moduli of attractive glasses.

The second reason for our interest is quite different, and involves a deeper analysis of the basis of the theory that we use to describe glassy systems. That is, rarely is MCT applied to a system simple enough that many issues can be worked out in detail, but at the same time that there be a complete knowledge of the ‘phase diagram’ and its dynamics. In our particular case, we shall study the square-well potential model, with very short-ranged attractions. Some of the basic predictions for structure and dynamics of this model have been worked out in detail within the MCT formalism . We have observed that the attractive glass introduces a new richness into the study, and it naturally becomes of interest to understand what presumptions MCT is making about such phenomena. We shall be interested in such questions as, what part of the physics drives the attractive glass transition within MCT, and whether this understanding can shed some light on the nature of the theory itself. Finally, as an aside we shall also seek to understand which physical characteristics of the colloidal particle and of the whole system sets the scale for the mechanical properties of a colloidal glass. Again, if MCT is correct, this will be of considerable practical interest.

In the light of these comments it is now possible to proceed to introduce the fundamental equations that define MCT, and the mechanical responses that we calculate from it.

## 2 Theory

Mode Coupling Theory provides a description of the structural relaxation of super-cooled liquids. The variables of interest are the normalized density correlators defined as,

$$\varphi _q(t)=\rho _𝐪^{}(0)\rho _𝐪(t)/|\rho _𝐪|^2$$
(1)
where $`\rho _𝐪(t)=_{j=1}^Ne^{i𝐪𝐫_j(t)}`$, with $`N`$ being the number of particles in the system. By using the Mori-Zwanzig formalism, it can be shown  that the equations of motion for the variables $`\varphi _q(t)`$ are,

$$\ddot{\varphi _q}(t)+\mathrm{\Omega }_q^2\varphi _q(t)+\nu _q\dot{\varphi _q}(t)+\mathrm{\Omega }_q^2_0^tm_q(tt^{})\dot{\varphi _q}(t^{})𝑑t^{}=0$$
(2)
for Newtonian dynamics. The two quantities $`\mathrm{\Omega }_q`$ and $`\nu _q`$ are respectively the characteristic frequency of the phonon-type motions of the fluid, and a term that describes instantaneous damping, the latter arising from the ‘fast’ contribution to the memory function. They are defined as,

$$\mathrm{\Omega }_q=\frac{q^2k_BT}{mS(q)}$$
(3)
$$\nu _q=\nu _1q^2$$
(4)
and typically one chooses $`\nu _1=1`$ in the calculations. Equations (2) are formally exact for a set of N particles.
To describe the dynamics of colloidal suspensions, they have been modified neglecting the inertia term and including the solvent contributions . Thus, we have

$$\dot{\varphi _q}(t)+q^2D_q^s\left\{\varphi _q(t)+_0^tm_q(tt^{})\dot{\varphi _q}(t^{})𝑑t^{}\right\}=0$$
(5)
where $`D_q^s`$ is the Brownian short-time diffusion.
The crucial approximation of MCT consists of giving an explicit factorized form for the memory kernel in eqs. (2) and (5) as,

$$m_q=\frac{1}{2}\frac{d^3k}{(2\pi )^3}V(𝐪,𝐤)\varphi _k(t)\varphi _{|𝐪𝐤|}(t)$$
(6)
and the vertex functions are the coupling constants of the theory,

$$V(𝐪,𝐤)=\frac{\rho }{q^4}\left(𝐪(𝐪𝐤)c_{|𝐪𝐤|}+𝐪𝐤c_k\right)^2S_qS_kS_{|𝐪𝐤|}$$
(7)
In the static limit $`t\mathrm{}`$, independently on the type of microscopic dynamics, the density correlators $`\varphi _q(t)`$ tend to a finite value $`f_q=\rho _q^{}(0)\rho _q(\mathrm{})/|\rho _q(0)|^2`$, known as the non-ergodicity factor, if the system is kinetically arrested. This loss of ergodicity for the density correlators is seen as the transition to a kinetic glassy state within MCT. Thus, the equations (2) become in the static limit,

$$\frac{f_q}{1f_q}=\frac{1}{2}\frac{d^3k}{(2\pi )^3}V(𝐪,𝐤)f_kf_{|𝐪𝐤|}$$
(8)
It is clear that $`f_q=0`$ always corresponds to a solution of equations (8). This corresponds to an ergodic state of the system in which the correlations decay for long time. For some critical values of the thermodynamic parameters(control parameters), such as temperature and density, there appear bifurcations of the solutions of equations (8), that produce non-zero solutions. These correspond to the non-ergodic states of the system, and given that there is no positional order in the system we identify these solutions as glasses. The bifurcations can be multiple, up to the number of control parameters of the system. Thus, when a bifurcation gives rise to more than two solutions of equations (8), there will exist multiple solutions with finite non-ergodicity factors. In this case, MCT predicts that only the state corresponding to the largest value of $`f_q`$ is the long-time limit solution of the equations .

With an input structure factor, we can now solve the MCT equations for the non-ergodicity parameter, and it is possible to calculate the ‘phase diagram’ of the system, localizing the regions in the thermodynamic parameter space where the system is in the fluid ($`f_q=0`$) or in the glassy state ($`f_q>0`$), and also some mechanical properties of the glass itself. Note that by ‘phase diagram’ we mean here that the fluid and glassy states of the system are identified, the latter being non-equilibrium states of matter.

A particularly interesting physical quantity is the elastic shear modulus, $`G^{}(\omega )`$. From the shear viscosity for a colloidal system , it is possible to evaluate the elastic shear modulus, within the MCT approximation, in the limit $`\omega 0`$, to give ,

$$G^{}(\varphi ,T)=\frac{d^3}{60\pi ^2}_0^{\mathrm{}}𝑑kk^4\left(\frac{dlnS_k}{dk}f_k\right)^2$$
(9)
where $`G^{}`$ is in units of $`(k_BT)/d^3`$, where $`T`$ is the temperature of the system and $`d`$ is the diameter of a particle.

Another property that can be examined by experimentalists is the longitudinal stress modulus $`m_0`$. We will discuss only the hydrodynamic approximation of this quantity , that is easily obtained taking the limits $`q0`$ and $`t\mathrm{}`$ in the equations (6) and (7), giving that,

$$m_{q=0}(\omega =0)=m_0=_0^{\mathrm{}}V_kf_k^2𝑑k$$
(10)
where the long wavelength limit of the vertex function is given by,

$$V_k=\rho S_{q=0}\left(\frac{S_kk}{2\pi }\right)^2\left[c_k^2+\frac{2}{3}\left(k\frac{c_k}{k}\right)c_k+\frac{1}{5}\left(k\frac{c_k}{k}\right)^2\right]$$
(11)
It is possible to relate $`m_0`$ to the velocity of sound in the glass (solid) compared to that in the liquid, at the transition line. Thus, we observe that the speed of sound in a liquid is given in terms of the compressibility of the liquid whilst the formation of a solid leads to a finite memory kernel at long times, and consequently an increment, $`m_0`$:

$$\frac{v_{\mathrm{}}}{v_0}=\sqrt{1+m_0}$$
(12)
where $`v_0`$ and $`v_{\mathrm{}}`$ are respectively the limiting low and high-frequency sound speeds.

Since the longitudinal modulus is a close equivalent of the shear modulus, but for extensional distortions, it is natural that we should seek to calculate and compare the two quantities. We would like to point out to experimentalists who may be interested in this field that the theory, as applied to ratios of quantities at the transition may be quite well modeled by the theory, and less dependent on model details. We therefore recommend some thought on how experiments of this type might be carried out.

We have solved  the MCT equations for the ideal glass-transitions (8) for a system of mono-disperse colloidal particles, interacting via a square well potential, with a very narrow-range attractive well, defined as

$$V(r)=\{\begin{array}{cc}\mathrm{}r<d& \\ \beta u_0d<r<d+\mathrm{\Delta }& \\ 0d+\mathrm{\Delta }<r& \end{array}$$
(13)
where $`\beta =\left(k_BT\right)^1`$, with $`k_B`$ the Boltzmann’s constant, and $`\mathrm{\Delta }`$ is the well-width, which can be related to the small parameter $`ϵ=\mathrm{\Delta }/(d+\mathrm{\Delta })`$.

To calculate numerically the equilibrium structure factors $`S(q)`$ of such a system, one can use the Mean-Spherical Approximation (MSA), Percus-Yevick, or other closure relations for the Ornstein-Zernike integral equation . Ref. compares many of the important features of the statics and dynamics of these systems using MSA and Percus-Yevick (PY) approximations. In view of the good agreements found between these closures, we shall here present only calculations based on the PY closure. The phase diagram for different values of the well-width of the potential has been presented in Ref. . In particular, for the mechanical results, we will focus our attention on the case where $`ϵ=0.03`$. Here we report in Figure 1 more details and present both the glass transitions as well as the liquid-gas spinodal, which may be regarded as an approximation to the equilibrium phase diagram. This provides the reader with an overview of where we expect to find all of the phenomena discussed later. Thus, in Fig. 1 we can see the curves labelled respectively $`B_1`$ representing the fluid-repulsive glass, and $`B_2`$ representing the fluid-attractive glass transition. In the inset, we have shown in more detail the attractive-repulsive glass transition curve with its end-point, labelled as $`A_3`$. Beyond this point, the two types of glasses become indistinguishable. Also, the dashed curve in the low volume fraction region represents the gas-liquid spinodal, calculated via an expansion in $`ϵ`$ around the Baxter potential  and in good agreement with numerical calculations for the square-well potential, at least to the right-hand side of the critical point (larger volume fractions). We comment here that the left-hand side branch of the spinodal is singular for a Baxter model, because it corresponds to complex values of the characteristic parameters $`\lambda `$, and also for such low packing fractions and temperatures the numerical PY solution of the square-well potential was not reliable. Thus, the only true meaningful branch of the spinodal in Fig. 1 is the ‘liquid’ branch, that corresponds approximately to $`\varphi 0.14`$, where we estimated that the critical point is located, corresponding to a critical temperature of approximately $`0.3`$.

It is interesting to note that, whilst the curve $`B_1`$ is almost vertical, $`B_2`$ is largely horizontal and the two meet forming a non-zero angle. The fact that the repulsive glass-liquid transition curve is vertical is ensured by the fact that glassification is driven only by the hard core, which lacks any energy scale. On the contrary, the attractive glass-liquid one, being fairly horizontal, implies that there is a single well-characterized energy scale that drives the glassification.

## 3 Results

To discuss our results we shall frequently refer to the phase diagram in Fig. 1.

Let us begin by looking at what is at first sight one of the more striking predictions of the calculations, the transition between two types of glasses and their merging at an end-point that has been labelled as the $`A_3`$ point , since it represents a third order singularity of the MCT equations (8).

The reader should be aware that, within MCT, on the transition lines, whether they be liquid-glass or glass-glass, the two states of matter do not have any difference in density or in structure and are differentiated in terms of the non-ergodicity factors $`f_q`$. Since this is a non-equilibrium property of the system, we note that there is essentially no equilibrium quantity that establishes the difference in phases, and the order parameter must therefore be composed of $`f_q`$.

Since the repulsive glass and the attractive glass illustrated in Fig. 1, and inset to that figure, are differentiated only by the changes of their $`f_q`$ it is therefore natural to ask what differences are implied in the shear modulus and other mechanical properties of these two glasses, and how these difference vanish as we approach the end-point where the two glasses become identical. We have examined this question for the example of the shear modulus.

Thus, in Figure 2 we plot $`G^{}`$ as a function of increasing temperature for two fixed volume fractions $`\varphi =0.539672`$ and $`\varphi =0.544052`$, both of them involving a crossing of the glass-glass transition curve. For the smaller of these two volume fractions, it is evident from the figure that, upon crossing the transition, there appears a sharp discontinuity in the shear modulus. For the larger one, that is very close to the end-point value packing fraction, $`\varphi _{A_3}`$, crossing the curve there is negligible difference between the shear moduli of the two glasses. Evidently, it is of interest to define $`\mathrm{\Delta }G^{}`$, the difference in $`G^{}`$ found in the two glasses at the transition, and to examine this as a function of $`\mathrm{\Delta }\varphi `$, and $`\mathrm{\Delta }T`$, respectively the differences in volume fraction and temperature from their end-point values, which we evaluated as $`\varphi _{A_3}0.5441`$ and $`T_{A_3}1.09975`$. We numerically find, over the whole range of the glass-glass transition, that the laws connecting these two quantities are,

$`\mathrm{\Delta }G^{}(\mathrm{\Delta }\varphi )^p;`$ $`\mathrm{\Delta }\varphi =\varphi \varphi _{A_3}`$ (14)
$`\mathrm{\Delta }G^{}(\mathrm{\Delta }T)^q;`$ $`\mathrm{\Delta }T=TT_{A_3}`$ (15)
where $`p`$ is approximately $`(0.33\pm 0.03)`$ and $`q`$ approximately $`(0.32\pm 0.08)`$. The errors in these estimates could be reduced with larger computational effort, but there do indeed appear to be power laws to high precision. These exponents can be explained by a simple argument. We know that near an $`A_2`$ singularity $`f_q`$ is a solution of a second order polynomial, while near an $`A_3`$ point it is solution of a cubic one. This implies,

$`f_qf_q^{A_2}(ϵ_{A_2})^{1/2}`$ (16)
$`f_qf_q^{A_3}(ϵ_{A_3})^{1/3}`$ (17)
where $`ϵ_{A_i}`$ can be both $`TT_{A_i}`$ or $`\varphi \varphi _{A_i}`$. Thus, evaluating the expression for $`G^{}`$ (9) in the proximity of the singularity $`A_3`$, this gives, in the leading order, an exponent $`1/3`$.

It will be of interest for experimentalists to seek the proposed phenomenon of glass-glass transition, and to deduce estimates of any exponents that arise. Amongst the most accessible of these is that of the shear modulus discussed above. However, we do note the reservation that it is likely that there is some associated structural and density relaxation at the glass-glass transition, and this may disturb the simple power law outlined above.

Evidently, a fairly practical comment that emerges from these results is that the repulsive glass stiffness with respect to shear is much smaller than the attractive glass one, as indicated by the large vertical discontinuity shown in the inset of Fig. 1 . This reflects the fact that particles in the attractive glass are bonded by the stickiness of the potential, whilst in the repulsive one there is no real bonding between them, thus implying that they are more easily broken apart under shear. Also, as expected, the attractive glass shear modulus increases considerably with the decrease of temperature, the attractions between particles becoming more relevant, whilst for the repulsive glass there are no significant changes with temperature, there being no energy scale involved in its formation.

Having outlined our results for the shear modulus differences between the two glasses, it is perhaps worth while to revisit the physical meaning and implications of these results. Firstly, all of the differences in mechanical properties here come from the non-ergodicity factor alone, rather than static structure since the two states have the same structure factor, that is essentially the structure factor of the liquid to which the MCT glass is referred. It is therefore interesting to note that the differences in $`f_q`$’s between the two types of glass, reflecting as they do the nature and efficiency of the relaxation processes at different length scales, lead to such large differences in moduli.

Let us now turn to another part of the phase diagram of these systems. We have earlier alluded to the work by Verduin and Dhont  for a very short-ranged attractive colloidal system, where they found some non-ergodic states at low temperature. These they interpreted as gel states, but in line with the ideas laid out in this paper they might also be viewed as attractive glasses, as it was already discussed in . In fact, the horizontal portion of the glass-liquid curve in their phase-diagram already is indicative that the energy scale is playing a leading role, rather than packing forces. Therefore, we believe that it will transpire that it is quite feasible to prepare and study attractive glasses in some detail.

In particular it will be possible to study many properties along the liquid-glass transition. We comment first that any attempts to design colloidal particles with the appropriate well shape will most likely involve some uncertainties. Thus, the obvious methods of coating the spheres with some attractive layer, or using depletion forces in polymer solutions, cannot hope to exactly reproduce the parameters and shape of the square well. Therefore, the curve of liquid-glass transitions corresponding to the square well, may be fit to experimental data. Though it may not be correct in all details, it should be reasonably accurate in its behaviour as a function of packing fraction of spheres, and in reproducing many of the attendant phenomena. Thus, irrespective of the extent of agreement between the square well potential and the effective potential that is ultimately tested in experiments, we may be fairly certain that some predictions of our theory will be more robust than others. Amongst these we might again include the typical evolution of the mechanical properties as a function of packing fraction.

Therefore it becomes of some interest to construct the evolution of $`G^{}`$ on the glass side of the liquid-glass transition across an extensive range of packing fractions encompassing all of what would be viewed as the attractive glass, to be able to compare these results with experiments. Thus, in Figure 3 we present the curve of $`G^{}`$ as a function of $`\varphi `$, along the attractive glass line, labelled $`B_2`$ in Fig. 1. We note that this curve extends to quite low volume fractions, and that it is linear in a large range of volume fractions. In fact this is true for all those values of $`\varphi `$ where the attractive interactions as considered to be completely dominant. Ultimately the curve turns downwards towards typical repulsive glass values of $`G^{}`$ when we go to higher volume fraction. Even so, we might have expected that the downward bending of the curve would have happened very close to the critical hard sphere packing fraction, $`\varphi 0.52`$, whereas it commences around $`\varphi 0.485`$. Very similar characteristics for the shear modulus, along the attractive glass transition line, have also been found for the attractive Yukawa potential .

Finally as a matter of curiosity let us draw attention to one aspect of the shear modulus plot in Figure 3. Consider the two circled points on the curve in Fig. 3. They correspond to states having approximately the same shear modulus but quite different packing fractions, i.e. $`\varphi _10.34`$ and $`\varphi _20.535`$. Thus, we represent in Fig. 4 the respective non-ergodicity parameters of the two states for comparison. It is interesting to note that, whilst the range of the two is almost the same as it should be since they both represent states of attractive glass, the one correspondent to the lower packing fraction exhibits a more pronounced non-ergodicity of the system at every length scale. The fact that we have the same shear modulus is a reflection that both the non-ergodicity and static structure factors are relevant for the modulus, and in this case we see that they compensate each other so that we can ‘build’, from the same system, two glasses with the same stiffness with respect to shear, but having completely different packings, and different structure. This would be an interesting phenomenon, if confirmed by experiments.

Now it is worth considering the MCT prediction that glasses can exist at volume fractions less than $`10`$%, as indicated in the phase diagram, Fig. 1. Certainly attractive glasses can and should form at much lower fractions than repulsive glasses, but we must not accept the results of MCT blindly. Note carefully the limitation that MCT does not permit the self-consistent relaxation of structure in the glass phase, and consequently the density cannot adapt to the new situation as the ergodicity is lost. The consequences of this have been outlined for the glass-glass transition, but one should be aware also for the possibility for qualitatively mistaken predictions in the region of a liquid-glass transition where there is a meta-stable liquid-gas phase separation, as in Fig. 1. Thus we must acknowledge the tendency for separation into more and less dense phases that could be superimposed on the liquid-glass transition curve $`B_1`$. That this does occur is not proven, but if it did, our MCT calculation would not accommodate it.

Given this concern, we have calculated the mean number of bonded spheres $`n_b`$ to a central sphere at the formation of the attractive glass, along the curve $`B_2`$. We consider as sufficient condition for bonding that the distance between two particles is up to the attractive well width. Thus, it is easy to estimate $`n_b`$ by carrying out the integration over the pair correlation function such that only those spheres falling within the attractive well are included in the integration. From the definition of the potential (13), we then have

$$n_b=_0^{(d+\mathrm{\Delta })}\rho g(r)𝑑𝐫$$
(18)
For convenience we also calculated one estimate of the mean coordination number $`n_c`$ around a sphere, irrespective of whether the spheres are bonded or not.

Nevertheless we have,

$$n_c=_0^r^{}\rho g(r)𝑑𝐫$$
(19)
where $`r^{}`$ corresponds to our estimation of the first peak of the pair correlation function $`g(r)`$.

The mean bonding number is well defined quantity, and is calculated exactly here, whereas the mean co-ordination number is not so well defined. The former may be relied on, the latter is used only as a comparison to the mean bonding number.

Both results are plotted in Figure 5 as a function of the volume fraction, along the curve $`B_2`$. It is striking that both curves are nearly linear. However, the interesting thing to note is that the mean bonding number at $`\varphi 0.10`$ is around $`1.5`$. Also, we can roughly estimate that only at $`\varphi 0.17`$ we find $`n_b2`$. Thus, it seems unlikely that extended mechanically stable structures can exist for such small bonding numbers. Indeed, a mean bonding number of $`2`$ would imply polymeric structures, and for bonding numbers a bit larger than $`2`$ we may have enough cross-links to establish a network, and finite shear moduli. Just how large the mean bonding number has to be before a finite shear modulus is obtained, we cannot say. Very simple arguments based on counting degrees of rotational and vibrational freedom, and requiring that there be no zero modes of the energy, indicate that the minimum number of bonds to form an extended structure for a chemical glass, should be about $`2.4`$ . However these arguments have not been shown to be relevant to the present case, so there is little more progress that can be made now.

The discussion should, however, alert the reader to the possibility that somewhere along the $`B_1`$ line, possibly when the bonding number is a little less than $`2`$, the reduction of any further bonds in the system might cause it to decompose, presumably into a less dense and a more dense phase. These remain open questions for the moment, but as MCT becomes applied more to glasses driven by attractions, it will be important to consider them in future.

Of course, another fundamental question that should be addressed in future is the applicability of idealized MCT to colloidal systems with attractive interactions. Indeed, whilst for hard-sphere type systems the idealized theory has been found in good agreement with experimental results, it is not yet clear whether for attractive systems, at lower densities, activated processes may become important, so that hopping should be included in the theory for a better description of this type of systems. Our present opinion is that the idealised theory is likely to be useful for moderate densities, where the bonding number is somewhat larger than 2. We do not preclude the possibility that, for much lower densities, the theory may still be suggestive, and indeed there have been already attempts to interpret these regimes .

Thus, all of this discussion leads us to seek a somewhat deeper understanding of exactly what the basis of the mode coupling predictions are. By this we mean, we seek to understand the essential features of the theory as it relates to attractive glasses. One of the points that we would like to probe a little more is to understand just how local MCT is in its understanding of the loss of ergodicity, and in its estimates of the main features of mechanical moduli. The simplest possible proposition would be that, essentially with knowledge of the number of bonds and strength of associations between nearest neighbour particles, we would estimate a correct MCT transition, and, perhaps, even estimate a correct mechanical response. We now seek to test this hypothesis. We shall examine this question in two steps, beginning with the calculation of the non-ergodicity factor that determines the glass transition itself, then progressing to the shear modulus.

For simplicity, we call $`q^{}`$ the wave-vector corresponding to the first peak of the structure factor. We then proceeded recalculating the liquid-glass transition using two approximations. Firstly, we used the correct short wavelengths behaviour of the structure factor, for $`(q>q^{})`$, in solving the MCT equations (8), but using, for the long wavelength part of the structure factor (i.e. $`q<q^{}`$), the correspondent part of it for a packing fraction that we chose arbitrarily. We both considered $`\varphi 0.31`$ and $`\varphi 0.39`$, obtaining no significant changes in the results. We found that the essential features of the trends and quantitative evolution of the transition are set only by the short-ranged part of the structure factor, then irrespectively of what we chose for the long wavelength part of it. Thus, length scales involving only the hard core and attractive well of a spherical particle are sufficient to determine the glass-liquid transition within MCT. Indeed, we investigated a few points along the curve $`B_2`$ and we also found that the transition temperatures are reproduced within a few percents error.

To confirm these results, we also tried to evaluate the curve $`B_2`$ using an opposite strategy. We then used the correct long wavelength part of the structure factor ($`q<q^{}`$), and as short-ranged part the correspondent for the same values of packing fractions as before. This lead not only to much larger errors in the glass-liquid curve $`B_2`$, but, what is more important, to the wrong evolution of that curve with volume fraction, this meaning that the curve was not only shifted but also changed in its shape. We may therefore conclude that, whatever the full content of MCT, one can estimate the location of the transition curve of the attractive glass to high precision by knowing only local information about the particles around a central particle. Long-ranged effects therefore appear not to play a central role in this aspect of the MCT description. The same observation had been made earlier for the Yukawa model and for the square-well potential, both solved within MSA . Thus, this is a general feature of the theory itself, independently on the model or the approximation chosen to solve it.

We now turn to another natural question about the predictions of MCT. After having analyzed what the dominant contributions are for the prediction of the location of the glass transition, we want to investigate what is important in the determination of the mechanical properties of the glass within the theory. To do so, we still refer to the results of the shear modulus $`G^{}`$ in Fig. 3.

Thus, we note that, in equation (9), there are evidently two important contributions to the shear modulus, the non-ergodicity factor, and the logarithmic derivative of the structure factor. We may ask if only one of these features provides the main contribution to the shear modulus. To check this we may proceed as before, first retaining the correct structure factor, and using the non-ergodicity factor of a chosen value of packing fraction ($`\varphi =0.39`$), and then vice versa. The question may be further refined by asking whether it is possible to further locate the main driving force of the shear modulus as being short ($`q>q^{}`$) or long ($`q<q^{}`$) length scales.

The answer is striking. With almost quantitative accuracy we find that the linear behaviour of the shear modulus with volume fraction, in the purely attractive region, which we can define being in the range of volume fractions between $`0.17`$ and $`0.4`$, originates solely in the short length-scale ($`q>q^{}`$) behaviour of the structure factor, as shown if Figure 6. We represent here the curve, already reported in Fig. 3 , with superimposed the data obtained by using the chosen $`f_q`$ and the true structure factor (crosses), and those obtained by using the true non-ergodicity factor and only the exact short length scale contribution for $`S(q)`$ (circles). Thus, providing the non-ergodicity factor is finite, and on the correct scale, the shear modulus is relatively insensitive across the whole attractive glass region to its details, being only shifted of a small amount. We note that in the figure all the sets of data are coincident at $`\varphi =0.39`$, since this is the chosen value of reference for the different cases. Since we find, as before, that only the short length scale picture of the system is really involved in the determination of $`G^{}`$ (circles), we can explain the linear behaviour of the shear modulus with packing fraction in the attractive region, being directly related to the number of bondings. Thus, we represent $`G^{}`$ in function of $`n_b`$ in the inset of Fig. 6, and we show that there also exists a linear relation between them.

Finally, we turn to analyze what happens for the phenomenology of the shear modulus at higher packing fractions, beyond what we believe is the pure attractive region, where $`G^{}`$ undergoes a rather rapid decline. We have pointed out already that, along the glass-liquid curve in Fig. 3 this rapid decrease of the shear modulus happens as we approach $`\varphi 0.48`$, that is a value quite below the close-packed structure appropriate for the hard core. It may be argued that this phenomenon is hardly surprising since we know that, finally, the shear modulus must decrease to values characteristic of a repulsive glass. These are much smaller since, as we have earlier pointed out, there are no attractive interactions, and no effectively bonded particles to the central one. However, it is nevertheless interesting to understand why the shear modulus softens quite dramatically at that particular packing fraction, and by what means the effective attractions are being screened in the system in this region of densities.

To understand this point, we have studied the thermodynamic pressure for a square well potential, given in ref. . In particular, we have examined the pressure and compressibility of the liquid along the glass-liquid curve, in order to exclude the possibility that there be any anomalies in the liquid itself in the relevant region. We found no anomalies within the PY approximation in the region of packings from approximately $`0.45`$ up to beyond the endpoint value $`\varphi _{A_3}`$. It is then interesting to note that, whilst the liquid is perfectly normal, the proximate glass, undergoes this softening at well defined values of $`\varphi `$ and $`T`$, corresponding to the decrease of the shear modulus. In essence we find that the softening of the glass, whilst it certainly originates in cancellations between hard core and attractive parts of the potential as the density increases, occurs mechanistically within the MCT memory kernel.

To be more precise, we have used the same strategy as in the purely attractive region, to examine the different contributions of the structure factor and of the non-ergodicity factor to the shear modulus. So, in Fig. 6 the data represented with crosses, corresponding to the correct structure factor and a fixed arbitrary non-ergodicity factor, are also represented at higher volume fractions. It is evident that their linear behaviour, entirely due to the structure factor, persists with increasing density. Conversely, by a similar analysis, where the structure factor was chosen arbitrary and the non-ergodicity factor was the correct one, we have found that the principal origin of the decrease in the modulus is the change in the non-ergodicity factor, since in this region it is changing from a characteristic attractive to a repulsive form. The smaller integral resulting in eq. (9) reflects the fact that the attractive glass is less mobile due to the formation of attractive bonds. Thus, it can only be the non-ergodicity factor, solution of the MCT equations (8), the responsible for the softening of the glass.

We now turn our attention to the behaviour of the longitudinal stress modulus. We have already shown that it presents a discontinuous behaviour, similar to the elastic shear modulus, upon crossing the glass-glass transition in . Now we show $`m_0`$ along the curve $`B_2`$, i.e. along the attractive glass-liquid transition line in Figure 7, on the glass side of the transition, in the same way as we did for the shear modulus. We conclude that, for increasing packing fraction, though the glass becomes more rigid with respect to shear, reaching a maximum at $`\varphi 0.48`$, the longitudinal modulus decreases continuously along that same curve, coming close to its characteristic liquid value very close to the point where the two curves $`B_1`$ and $`B_2`$ meet, corresponding approximately to $`\varphi 0.536`$. At first sight this seems counter-intuitive since we expect that the extensional rigidity should increase with density, at least up to that density where the elastic shear constantly increases.

If we examine the longitudinal modulus as a function of density, for temperatures that cause the system to be within the repulsive glass region only, the normal expectation of increased longitudinal modulus with increasing density is confirmed, as illustrated in the $`a`$-part of figure 8. In the $`b`$-part of fig. 8, instead, we have the opposite behaviour of the longitudinal modulus with density for a much lower temperature than in the $`a`$-part. Indeed, from eq. (10), it is clear that the modulus depends on the zero-momentum limit of the structure factor, appearing as a prefactor, and the data correspond to a temperature for which the system passes very closely to the critical point, or spinodal curve, of the underlying liquid-gas transition. In fact, MCT implies that the long length scale associated with the proximate liquid, when it is near a critical point, or otherwise the underlying spinodal, is quenched into the solid glass because both liquid and glass have the same structure factor. Within MCT the large modulus in the attractive glass region therefore derives essentially from this large quenched correlation length in the glass, rather that any microscopic interaction. We note that there is no such quenched long length-scale dependence implied by the shear modulus, either in the formula (9), or from the results showed in Figs. 3 and 6. In fact, as we showed earlier, the shear modulus in the attractive region is mainly determined by nearest neighbour adhesions. The striking difference in the behaviours of the two modulii is an interesting prediction of MCT, when applied to attractive glasses. Possibly some of this variation originates because long length scale structures are quenched into the glass. Whilst we do not know if it is true in nature, the strong distinction between longitudinal modulus, and shear modulus along the liquid-glass curve, will be readily testable in experiments on colloidal systems.

## 4 Conclusions

We have studied the kinetic glass transitions of particles with a model square well potential using Mode Coupling Theory. This interaction potential should, we believe, be a reasonable approximation to that found in particles with grafted chains and systems with strong depletion interactions. All of the non-trivial results arise when the attractive piece of the interaction potential is strong, but of very short range, and we have studied a typical example of this type. Based on this, we have been able to describe some mechanical properties of two types of colloidal glass, the one resulting mainly from attractive and the other mainly from repulsive interactions between particles. We and others  have earlier proposed a distinction between the two kinds of glasses based on their dynamical behaviour However, in this paper we have shown that this difference may be probed experimentally using the difference in mechanical properties, and that there are exponents that describe the merging of these glasses into a single glass beyond some critical density.

In general terms, glasses dominated by attractions have a stronger rigidity under shear than those originating simply from packing forces.

Also, we have studied the behaviour of the zero-frequency elastic shear modulus $`G^{}`$ and the longitudinal stress modulus $`m_0`$ along the liquid-attractive glass transition curve. Besides being of intrinsic interest, this part of the study was chosen in the belief that comparison to experiments will be more reliable along this curve, and less dependent on model potential details. The predictions are striking. This is a positive feature since even confirmation of trends will be of some interest in evaluating the MCT approach to these systems.

At lower packing fraction, where we believe attractions to be completely dominant, there is a linear increase of the shear modulus with packing. This reflects the fact that also the mean number of bonds to a particle in the system is fairly linear in $`\varphi `$ and the only relevant physical mechanism to determine the shear modulus in this region of the phase diagram occurs at short length scale. We have argued that here the shear modulus is determined by nearest neighbour adhesions. As an aside, we note also that within MCT the attractive glass-liquid transition curve itself is essentially determined by the the same factors as those determining the shear modulus.

On the contrary, at higher values of packing fraction, the shear modulus decreases quite dramatically towards typical repulsive glass values. This phenomenon originates from the changes in the non-ergodicity factor that compete in equation (9) with those of the structure factor and lead the system to soften with respect to shear. The implication is that there should be a maximum of the shear modulus along the liquid-glass transition.

Now we turn to the longitudinal modulus of the system. Here the result is quite different from that for the shear modulus. Along the glass-liquid transition curve $`m_0`$ decreases with increasing density. We have observed that the transition curve, passes close to the critical point and some portion of the spinodal of the underlying liquid-gas transition, and consequently the input structure factor at long wavelength is very large. The prefactor in formula (11) therefore becomes very large. The physical implication is that large scale density fluctuations quenched into the glass cause a large increase in the longitudinal modulus, and the longitudinal modulus is dominated by these for much of the attractive glass regime.

We note, in passing that MCT does not always reliably predict the stability of the attractive glass. By this we mean that when we check the mean bonding number to a particle in what has been predicted to be a glass, we sometimes find that this number is even less than 2. We conclude that in these cases, MCT is over-emphasizing the stability of the glass, probably for a variety of reasons. In any case, we use this independent calculation as a rough check to exclude regions of MCT glass that are clearly unphysical. It will be important to address these issues in future.

There is no doubt that the model and the means by which we have studied it are simple, and there are numerous limitations implied thereby. Nevertheless, the model is in some sense canonical in that it contains the essential physical input of strong short-ranged interaction, and repulsion. We believe that it encompasses many essential ideas regarding the mechanical properties of colloidal glasses and we have laid these out for consideration by further theoretical, but mainly experimental researches.

We thank W. Götze and M. Fuchs for discussions and suggestions. F.S. and P.T. are supported by INFM-PRA-HOP and MURST PRIN98, E.Z., G.F. and K.D. by INCO-Copernicus Grant IC15CT96-0756 and COST P1.


format: text

text:
# RELATIVISTIC THERMAL BREMSSTRAHLUNG GAUNT FACTOR FOR THE INTRACLUSTER PLASMA. III. FREQUENCY INTEGRATED EMISSIVITY

## 1 INTRODUCTION

Nozawa, Itoh, & Kohyama (1998) have recently carried out accurate calculations on the relativistic thermal bremsstrahlung Gaunt factor for the intracluster plasma. Their calculation is based on the method of Itoh and his collaborators (Itoh, Nakagawa, & Kohyama 1985; Nakagawa, Kohyama, & Itoh 1987; Itoh, Kojo, & Nakagawa 1990; Itoh et al. 1991,1997). In calculating the relativistic thermal bremsstrahlung Gaunt factor for the high-temperature, low-density plasma, Nozawa, Itoh, & Kohyama (1998) have made use of the Bethe-Heitler cross section (Bethe & Heitler 1934) corrected by the Elwert factor (Elwert 1939). They have also calculated the Gaunt factor using the Coulomb-distorted wave functions for nonrelativistic electrons following the method of Karzas & Latter (1961). They have presented the results of the numerical calculations in the form of extensive numerical tables. The grid of their numerical table is made fine enough so that smooth interpolations can be carried out for the applications of the results to the analyses of various astrophysical observational data.

Other references on the calculation of thermal bremsstrahlung Gaunt factor include Culhane (1969), Culhane & Acton (1970), Raymond & Smith (1977), Gronenschild & Mewe (1978), Mewe, Lemen, & van den Oord (1986), and Carson (1988).

Itoh et al. (2000) presented accurate analytic fitting formulae which summarize the results of the extensive numerical tables presented by Nozawa, Itoh, & Kohyama (1998). Analytic fitting formulae can be readily implemented in the computer programs for the analyses of the observational data. Therefore, we recommend these analytic fitting formulae be used widely by the observers for the analyses of their observational data.

In this paper we will present accurate analytic fitting formulae for the frequency integrated emissivity due to the relativistic as well as nonrelativistic thermal bremsstrahlung which is relevant to the intracluster plasma. The fitting is carried out for $`1Z_j28`$, $`6.0\mathrm{log}T8.5`$, where $`Z_j`$ is the charge number of the ion, and $`T`$ is the electron temperature in kelvins. The present analytic fitting formula will be useful for the analysis of the X-ray emission which comes from the intracluster plasma as well as the other X-ray sources.

The present paper is organized as follows. We will give the analytic fitting formulae for the frequency integrated Gaunt factor for the relativistic thermal bremsstrahlung in $`\mathrm{\S }`$2. We will give the analytic formula for the frequency integrated Gaunt factor for the nonrelativistic thermal bremsstrahlung in $`\mathrm{\S }`$3. Concluding remarks will be given in $`\mathrm{\S }`$4.

## 2 FREQUENCY INTEGRATED GAUNT FACTOR FOR RELATIVISTIC THERMAL BREMSSTRAHLUNG

The thermal bremsstrahlung emissivity is expressed in terms of the thermally averaged relativistic Gaunt factor $`g_{Z_j}`$ (Nozawa, Itoh, & Kohyama 1998) by

$`<W(u)>du`$ $`=`$ $`1.426\times 10^{27}g_{Z_j}(T,u)n_en_jZ_j^2T^{1/2}`$ (1)
$`\times `$ $`e^udu\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^3,`$
$`u`$ $``$ $`{\displaystyle \frac{\mathrm{}\omega }{k_BT}},`$ (2)
where $`\omega `$ is the angular frequency of the emitted photon, $`T`$ is the temperature of the electrons (in kelvins), $`n_e`$ is the number density of the electrons (in cm-3), and $`n_j`$ is the number density of the ions with the charge $`Z_j`$ (in cm-3). In Nozawa, Itoh, & Kohyama (1998), the nonrelativistic Gaunt factor $`g_{NR}`$ has been also calculated. This Gaunt factor is exact in the low-temperature limit. At intermediate temperatures, it has been confirmed by Nozawa, Itoh, & Kohyama (1998) that the relativistic Gaunt factor which has been calculated with the use of the Bethe-Heitler cross section corrected by the Elwert factor shows excellent agreement with the nonrelativistic exact Gaunt factor. At higher temperatures, the nonrelativistic Gaunt factor deviates from the relativistic Gaunt factor because of the insufficiency of the nonrelativistic approximation. These two Gaunt factors have been tabulated in Nozawa, Itoh, & Kohyama (1998).

At sufficiently high temperatures, we adopt the relativistic Gaunt factor. At sufficiently low temperatures, we adopt the nonrelativistic exact Gaunt factor. At intermediate temperatures, these two Gaunt factors coincide with each other for small values of $`Z_j`$. For larger values of $`Z_j`$, the two Gaunt factors show small discrepancies even at intermediate temperatures. Therefore, we generally interpolate between the two Gaunt factors smoothly at intermediate temperatures. To be more precise, we find the point at which the discrepancy between the two Gaunt factors (for fixed values of $`Z_j`$ and u) is the smallest as a function of the temperature. Then we interpolate between the two Gaunt factors smoothly using a sine function. The temperature range for the interpolation is $`\mathrm{\Delta }`$ $`\mathrm{log}T`$ = $`\pm `$ 0.1 to $`\pm `$ 0.5 with respect to the central temperature at which the discrepancy is the smallest depending on the minimum value of the discrepancy. The analytic fitting formulae of $`g_{Zj}(T,u)`$ for $`Z_j=128`$, $`6.0\mathrm{log}T8.5`$, $`4.0\mathrm{log}u1.0`$ have been presented in Itoh et al. (2000).

Now we will integrate the emissivity over the whole frequency range. Thus we obtain

$`W`$ $``$ $`{\displaystyle _0^{\mathrm{}}}<W(u)>𝑑u`$ (3)
$`=`$ $`1.426\times 10^{27}g_{Z_j}(T)n_en_jZ_j^2T^{1/2}\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^3,`$
$`g_{Z_j}(T)`$ $``$ $`{\displaystyle _0^{\mathrm{}}}e^ug_{Z_j}(T,u)𝑑u.`$ (4)
In Figure 1 we show the frequency integrated Gaunt factor $`g_{Z_j}(T)`$ as a function of temperature for various values of $`Z_j`$.

We give an analytic fitting formula for $`g_{Z_j}(T)`$. The range of the fitting is $`1Z_j28`$, $`6.0\mathrm{log}T8.5`$. We express the frequency integrated Gaunt factor by

$`g_{Z_j}(T)`$ $`=`$ $`{\displaystyle \underset{i,k=0}{\overset{10}{}}}a_{ik}z^it^k,`$ (5)
$`z`$ $``$ $`{\displaystyle \frac{1}{13.5}}\left(Z_j14.5\right),`$ (6)
$`t`$ $``$ $`{\displaystyle \frac{1}{1.25}}\left(\mathrm{log}T7.25\right).`$ (7)
The coefficients $`a_{ik}`$ are presented in Table 1. The accuracy of the fitting is generally better than 0.1%.

## 3 FREQUENCY INTEGRATED GAUNT FACTOR FOR NONRELATIVISTIC THERMAL BREMSSTRAHLUNG

The thermal bremsstrahlung emissivity in the nonrelativistic limit is expressed in terms of the nonrelativistic exact Gaunt factor $`g_{\mathrm{NR}}`$ (Nozawa, Itoh, & Kohyama 1998) by

$`<W(u)>_{\mathrm{NR}}du`$ $`=`$ $`1.426\times 10^{27}g_{\mathrm{NR}}(\gamma ^2,u)n_en_jZ_j^2T^{1/2}`$ (8)
$`\times `$ $`e^udu\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^3,`$
$`u`$ $``$ $`{\displaystyle \frac{\mathrm{}\omega }{k_BT}},`$ (9)
$`\gamma ^2`$ $``$ $`{\displaystyle \frac{Z_{j}^{}{}_{}{}^{2}\mathrm{Ry}}{k_BT}}=Z_{j}^{}{}_{}{}^{2}{\displaystyle \frac{1.579\times 10^5\mathrm{K}}{T}},`$ (10)
where $`\omega `$ is the angular frequency of the emitted photon, $`T`$ is the temperature of the electrons (in kelvins), $`n_e`$ is the number density of the electrons (in cm-3), and $`n_j`$ is the number density of the ions with the charge $`Z_j`$ (in cm-3). It should be noted that the thermal bremsstrahlung emissivity in the nonrelativistic limit is a function of $`\gamma ^2`$ and $`u`$ only. It does not depend on $`Z_j`$ and $`T`$ separately, but on the ratio $`Z_{j}^{}{}_{}{}^{2}/T`$. This is a remarkable fact for nonrelativistic electrons.

The analytic fitting formula for the nonrelativistic exact Gaunt factor $`g_{NR}(\gamma ^2,u)`$ for the ranges $`3.0\mathrm{log}\gamma ^22.0`$, $`4.0\mathrm{log}u1.0`$ has been presented in Itoh et al. (2000).

Now we will integrate the nonrelativistic emissivity over the whole frequency range. Thus we obtain

$`W_{NR}`$ $``$ $`{\displaystyle _0^{\mathrm{}}}<W(u)>_{NR}𝑑u`$ (11)
$`=`$ $`1.426\times 10^{27}g_{NR}(\gamma ^2)n_en_jZ_j^2T^{1/2}\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^3,`$
$`g_{NR}(\gamma ^2)`$ $``$ $`{\displaystyle _0^{\mathrm{}}}e^ug_{NR}(\gamma ^2,u)𝑑u.`$ (12)
In Figure 2 we show the frequency integrated nonrelativistic Gaunt factor $`g_{NR}(\gamma ^2)`$ as a function of $`\gamma ^2`$. A similar graph has been shown in Karzas & Latter (1961).

We give an analytic fitting formula for $`g_{NR}(\gamma ^2)`$. The range of the fitting is $`3.0\mathrm{log}\gamma ^22.0`$. We express the frequency integrated nonrelativistic Gaunt factor by

$`g_{NR}(\gamma ^2)`$ $`=`$ $`{\displaystyle \underset{i=0}{\overset{10}{}}}b_i\mathrm{\Gamma }^i,`$ (13)
$`\mathrm{\Gamma }`$ $``$ $`{\displaystyle \frac{1}{2.5}}\left(\mathrm{log}\gamma ^2+0.5\right).`$ (14)
The coefficients $`b_i`$ are presented in Table 2. The accuracy of the fitting is generally better than 0.1%.

## 4 CONCLUDING REMARKS

We have presented accurate analytic fitting formulae for the frequency integrated Gaunt factors for relativistic as well as nonrelativistic thermal bremsstrahlung. The analytic fitting formulae have been constructed from the numerical results of the calculation reported in Nozawa, Itoh, & Kohyama (1998). The accuracy of the fitting is generally better than 0.1%. The present fitting formulae can be used widely for the analysis of thermal bremsstrahlung radiation.

We thank Professor Y. Oyanagi for allowing us to use the least square fitting program SALS. This work is financially supported in part by the Grant-in-Aid of Japanese Ministry of Education, Science, Sports, and Culture under the contract #10640289.

Figure Legends

* Fig. 1. Frequency integrated Gaunt factor $`g_{Z_j}(T)`$ for relativistic thermal bremsstrahlung as a function of temperature for various values of $`Z_j`$.
* Fig. 2. Frequency integrated Gaunt factor $`g_{NR}(\gamma ^2)`$ for nonrelativistic thermal bremsstrahlung as a function of $`\gamma ^2`$.


text:
# Optical and infrared interferometry

## 1. Introduction

It is somewhat disappointing that a meeting on “Galaxies and their Constituents at the Highest Angular Resolutions” contains so few results from optical interferometry111In this review, the term “optical” is also meant to include infrared wavelengths.. To some extent, this reflects the interests of those on the Scientific Organising Committee, but it must also be taken as a sign that optical interferometry has yet to deliver on its promises. Why is this? One reason is surely the difficultly of the technology and its vulnerability to Hofstadter’s Law (1980):

> Hofstadter’s Law — It always takes longer than you expect, even when you take into account Hofstadter’s Law.

But progress is being made, and there is a large number of interferometry projects in various stages of development and operation. The following is the current list of operating interferometers, listed in the order in which they obtained first fringes:

* GI2T/REGAIN (Grand Interféromètre à 2 Télescopes).
* SUSI (Sydney University Stellar Interferometer)
* COAST (Cambridge Optical Aperture Synthesis Telescope)
* ISI (Infrared Spatial Interferometer)
* FLUOR (Fiber Linked Unit for Optical Recombination)
* IOTA (Infrared-Optical Telescope Array)
* NPOI (Navy Prototype Optical Interferometer)
* PTI (Palomar Testbed Interferometer)
* CHARA (Center for High Angular Resolution Astronomy)

The following are under construction (in alphabetical order, since I am not game to predict which will be first to get fringes!):

* Keck Interferometer
* LBT (Large Binocular Telescope)
* MIRA-II (Mitaka IR Array)
* VLTI (Very Large Telescope Interferometer)

Finally, there are planned space missions, including SIM, TPF and DARWIN.

More details of all these projects can be found via links from the Web-based “Optical Long Baseline Interferometry News,” currently maintained by Peter Lawson222http://huey.jpl.nasa.gov/olbin/. Here, I simply want to stress the number and range of the projects, and to point out that many of them are producing scientific results. The proceedings of recent conferences (Unwin & Stachnik 1999; Lena & Quirrenbach 2000) testify to the vigorous activity in this field. In Section 3., I highlight some areas in which recent technological developments have been made. First, however, it is appropriate to make some general remarks about optical interferometry.

## 2. Interferometers measure fringes!

Realising that interferometers measure fringes is the key to understanding their capabilities and limitations (see, e.g., Bedding 1997; Quirrenbach 1997). To observe an object, an interferometer must be able to detect fringes and track them. This requires that the object has a sufficiently compact bright component, or that a suitable reference source (e.g., a star) is located nearby on the sky.

The main observable is the fringe contrast – the visibility – which is a dimensionless number between zero and one that indicates the extent to which a source is resolved on the baseline being used. Importantly, the signal-to-noise with which the visibility can be measured is a function of $`NV^2`$, where $`N`$ is the number of photons detected per sub-aperture per integration time, and $`V`$ is the measured visibility (including the inevitable reduction by atmospheric and instrumental effects). This dependence on $`NV^2`$ implies that interferometry becomes increasingly difficult for faint sources, but also for those with complex structures. For example, a drop in visibility from 1 to 0.1 as a source becomes resolved is equivalent, in terms of signal-to-noise, to a drop in brightness by a factor of 100 (5 magnitudes)!

The $`NV^2`$ limit can be addressed in several ways:

* Using larger sub-apertures. However, this requires adaptive optics if the telescope diameter becomes much larger than $`r_0`$, the average scale over which the wavefront is flat.
* Bootstrapping to increase integration times, which can be achieved by tracking fringes over a connected series of short baselines to allow low-visibility fringes to be measured on the longest. This approach requires lots of telescopes, which must be deployed rather wastefully in terms of $`(u,v)`$ coverage.
* Tracking fringes on a point source to increase integration time on the target. This is best done with a dual-feed system.

## 3. Recent developments

### 3.1. Fibres and integrated optics

Single-mode optical fibres have been used very successfully for spatially filtering light beams (Coude Du Foresto et al. 1997). This allows visibility amplitudes to be measured to high precision (about one percent), by reducing the sensitivity to atmospheric fluctuations. So far, results have been demonstrated with IOTA in the $`K`$ band (2.2 micron; Perrin et al. 1999) and now in the $`L`$ band (3.75 micron; Mennesson et al. 1999). It seems clear that fibres have an important role to play in interferometry, especially in the infrared.

A related development is in the field of integrated optics, which potentially allows large tables of bulk optics to be replaced by miniature devices. White-light interference has been achieved in the laboratory with two (Berger et al. 1999) and three beams (Haguenauer et al. 2000). It is too soon to estimate the impact of this technology on optical interferometers, with applications in space perhaps being the most likely.

### 3.2. Mid-infrared interferometry

The only interferometer operating in the mid-infrared is the ISI, which has a heterodyne system. A recent development has been the application of ISI to spectral-line observations of molecules around giant stars (Monnier et al. 2000). Both VLTI and Keck are targeting the mid-infrared, which is very attractive because the 8–10-m apertures are essentially diffraction-limited at 10–20 microns. A lot of excellent science can be expected from those instruments, provided the problems with thermal background can be overcome, but the amount of observing time will obviously be limited.

### 3.3. Imaging

Imaging is often seen as the main goal of interferometry, despite the fact that a lot of scientific questions can be answered without reconstructing an image. We must not expect high-quality VLA-type images from optical interferometers, given the much smaller number of array elements and the much larger atmospheric phase fluctuations. Instead, the correct parallel to draw is with VLBI, in which closure-phase imaging is a well-established technique. The first crude images from optical arrays have been produced for binary stars (Baldwin et al. 1996; Hummel et al. 1998) and barely-resolved single stars (Young et al. 2000).

The key to good imaging is a large number of array elements. This is well known by radio astronomers, and the potential in the optical is shown by the spectacular images produced with aperture-masking of single telescopes such as the Keck (Tuthill et al. 2000). Images from current long-baseline arrays will not be this good, since none has more that six elements. The next generation of imaging interferometers will need at least ten, and preferably fifteen elements.

### 3.4. Differential phase measurements

Apart from fringe visibility, the other main observable for an interferometer is fringe phase. With three or more telescopes, closure phases can be used to reconstruct complex objects, as mentioned above. But even with two telescopes, useful phase measurements can be made. These involve measuring the differential phase of one object relative to another (or one object relative to itself at different wavelengths), which allows very precise astrometry. A nice feature of differential phase measurements is that much interesting science can be done with point sources (i.e., stars), so that the $`NV^2`$ limitation is overcome. Examples include parallaxes and proper motions, such as small motions induced by the gravitational pull of unseen planets and, for dual-wavelength observations, the shift in photo-centre caused by the presence of a “hot Jupiter” planet.

The key technological element is a dual feed, in which two stars are observed simultaneously so that their differential phase can be measured. Such a system also allows dual-wavelength operation, in which a single star is observed at two wavelengths with the two feeds. Finally, the second feed can instead be used for fringe tracking on a reference star, giving longer integration times on a fainter and/or more complex source. All these techniques have been demonstrated for the first time with the PTI (Lane et al. 2000) and, although problems remain, it appears that plans to apply these techniques to the Keck Interferometer and space missions are feasible.

### 3.5. Nulling

Methods for using an interferometric null to allow a faint object (e.g., a planet) to be seen next to a bright star are being developed by many groups (Hinz et al. 1998; Guyon et al. 1999; Böker & Allen 1999; Velusamy et al. 1999; Ollivier et al. 2000; Wallace et al. 2000). Deep and stable nulls have now been demonstrated in the laboratory and with telescopes, and nulling capabilities will be included in the Keck Interferometer. The scientific aims include direct detection of exoplanets, and also of exo-zodiacal dust. Eventually, it is hoped that space interferometers will provide images and even spectra of Earth-like exoplanets. And that somewhat optimistic – but perhaps achievable – scenario seems a good place to end this review!

#### Acknowledgments.

For travel support, I thank the Australian Research Council and the Science Foundation for Physics in the University of Sydney.

## References

Baldwin, J. E., Beckett, M. G., Boysen, R. C., Burns, D., Buscher, D. F., Cox, G. C., Haniff, C. A., Mackay, C. D., Nightingale, N. S., Rogers, J., Scheuer, P. A. G., Scott, T. R., Tuthill, P. G., Warner, P. J., Wilson, D. M. A., & Wilson, R. W., 1996, A&A, 306, L13.

Bedding, T. R., 1997, In Paresce, F., editor, ESO workshop on Science with the VLT Interferometer, page 3. Berlin: Springer-Verlag. astro-ph/9609092.

Berger, J. P., Rousselet-Perraut, K., Kern, P., Malbet, F., Schanen-Duport, I., Reynaud, F., Haguenauer, P., & Benech, P., 1999, A&AS, 139, 173.

Böker, T., & Allen, R. J., 1999, ApJS, 125, 123.

Coude Du Foresto, V., Ridgway, S., & Mariotti, J.-M., 1997, A&AS, 121, 379.

Guyon, O., Roddier, C., Graves, J. E., Roddier, F., Cuevas, S., Espejo, C., Gonzalez, S., Martinez, A., Bisiacchi, G., & Vuntesmeri, V., 1999, PASP, 111, 1321.

Haguenauer, P., Severi, M., Schanen-Duport, I., Rousselet-Perraut, K., Berger, J. P., Duchêne, Y., Lacolle, M., Kern, P., Malbet, F., & Benech, P., 2000, In Lena, P. J., & Quirrenbach, A., editors, Interferometry in Optical Astronomy, Proc. SPIE vol. 4406, page 1107.

Hinz, P. M., Angel, J. R. P., Hoffmann, W. F., McCarthy, D. W., McGuire, P. C., Cheselka, M., Hora, J. L., & Woolf, N. J., 1998, Nat, 395, 251.

Hofstadter, D. R. Gödel, Escher, Bach: an Eternal Golden Braid. Harmondsworth: Penguin, 1980.

Hummel, C. A., Mozurkewich, D., Armstrong, J. T., Hajian, A. R., Elias, N. M., & Hutter, D. J., 1998, AJ, 116, 2536.

Lane, B. F., Colavita, M. M., Boden, A. F., & Lawson, P. R., 2000, In Lena, P. J., & Quirrenbach, A., editors, Interferometry in Optical Astronomy, Proc. SPIE vol. 4406, page 452.

Lena, P. J., & Quirrenbach, A., editors. Interferometry in Optical Astronomy, Proc. SPIE vol. 4406, 2000.

Mennesson, B., Mariotti, J.-M., Coudé Du Foresto, V., Perrin, G., Ridgway, S., Ruilier, C., Traub, W. A., Carleton, N. P., Lacasse, M. G., & Mazé, G., 1999, A&A, 346, 181.

Monnier, J. D., Fitelson, W., Danchi, W. C., & Townes, C. H., 2000, ApJS, 129, 421.

Ollivier, M., Mariotti, J. M., Sekulic, P., Michel, G., Léger, A., Bouchareine, P., Brunaud, J., Coudé Du Foresto, V., Mennesson, B., Borde, P., Amy-Klein, A., Vanlerberghe, A., Lagage, P.-O., Artzner, G. E., & Malbet, F., 2000, In Lena, P. J., & Quirrenbach, A., editors, Interferometry in Optical Astronomy, Proc. SPIE vol. 4406, page 354.

Perrin, G., Coudé Du Foresto, V., Ridgway, S. T., Mennesson, B., Ruilier, C., Mariotti, J. M., Traub, W. A., & Lacasse, M. G., 1999, A&A, 345, 221.

Quirrenbach, A., 1997, In Paresce, F., editor, ESO workshop on Science with the VLT Interferometer, page 35. Berlin: Springer-Verlag.

Tuthill, P. G., Monnier, J. D., Danchi, W. C., & Wishnow, 2000, pasp, 116, 2536.

Unwin, S., & Stachnik, R., editors. Working on the Fringe: Optical and IR Interferometry from Ground and Space, volume 194, 1999. ASP Conf. Ser.

Velusamy, T., Beichman, C. A., & Shao, M., 1999, In Unwin, S., & Stachnik, R., editors, Working on the Fringe: Optical and IR Interferometry from Ground and Space, volume 194, page 430. ASP Conf. Ser.

Wallace, K., Hardy, G., & Serabyn, E., 2000, Nat, 406, 700.

Young, J. S., Baldwin, J. E., Boysen, R. C., Haniff, C. A., Lawson, P. R., Mackay, C. D., Pearson, D., Rogers, J., St.-Jacques, D., Warner, P. J., Wilson, D. M. A., & Wilson, R. W., 2000, mnras, 315, 635.


text:
# Maximally incompressible neutron star matter

## I introduction

The problem of the maximum possible mass of neutron stars has a long history. Baade and Zwicky in 1934 proposed a star composed of densely packed neutrons as the final state resulting from the supernova process . Oppenheimer and Volkoff in 1939 demonstrated a star composed of noninteracting neutrons is supported against gravitational collapse to a black hole by the Fermi degeneracy pressure only for stellar masses up to 0.72 solar masses ($`M_{}`$) . The discovery of the first pulsar in 1967 and the realization pulsars are highly magnetized rotating neutron stars triggered substantial theoretical work towards understanding the structure of neutron stars. Pulsar masses can be measured from observations of Doppler shifts of periodic emissions from neutron stars in binary systems. Some 20 radio pulsar masses have been accurately determined so far, with most masses clustering around 1.4 $`M_{}`$ . This is well beyond the value of 0.72 $`M_{}`$ computed for noninteracting neutrons, so the short-range repulsion of the strong nuclear force must make a substantial contribution to the pressure supporting neutron stars against gravitational collapse. There is also evidence for neutron stars with substantially larger masses than the radio pulsars due to accretion of matter from a binary companion. The mass of the x-ray pulsar in the Vela X-1 binary has been estimated to be 1.9 $`M_{}`$ , and the x-ray burster Cygnus X-2 has been determined to be 1.8 $`M_{}`$ . Quasiperiodic oscillations in the x-ray emissions from neutron stars accreting matter from low mass companions have been argued to imply a mass of up to 2.3 $`M_{}`$ for the neutron stars in these systems .

An accurate theoretical determination of the maximum possible mass of neutron stars is of practical interest for identifying as a black hole any compact object with a larger mass. The maximum mass value remains uncertain because the equation of state of neutron star matter is not well understood at the high density values found in neutron star interiors. Modern methods develop the equation of state by fitting experimental nucleon-nucleon scattering data and properties of light nuclei to two and three-body interaction potentials . Confidence in the results of these methods is high near normal nuclear density (energy density = 152 MeV/fm3 = 2.7 X 1014 g/cm3, baryon density = 0.16 fm-3) because, at low density, experimental verification is possible with existing facilities. The validity of extrapolating the experimentally verified low density equation of state to high density, or as an alternative theoretically modeling the physical processes thought to occur at high density, is uncertain because of the lack of high density laboratory data. The sources of the uncertainties in the high density equation of state are incomplete knowledge of the three-nucleon interactions, the contributions of mesonic and other baryonic besides nucleonic degrees of freedom, and a possible quark deconfinement transition .

In light of these uncertainties in the exact form of the high density equation of state, several authors have derived an upper limit on the maximum neutron star mass using only general restrictions on the equation of state. Oppenheimer and Volkoff suggested , and Rhoades and Ruffini rigorously proved , the maximum mass of a stable neutron star results when the stiffest equation of state compatible with thermodynamic stability and causality is used. Rhoades and Ruffini derived an upper limit on the maximum neutron star mass of 3.2 $`M_{}`$ by using the experimentally verified Harrison-Wheeler equation of state  up to 1.7 times normal nuclear energy density, and the stiffest equation of state compatible with causality at high density. They defined a casual equation of state to be one having an adiabatic sound speed less than the speed of light, so they took the most incompressible equation of state to have the sound speed equal to light speed. Kalogera and Baym  updated the result of Rhoades and Ruffini, using an equation of state developed by Wiringa, Fiks, and Fabrocini  up to twice normal nuclear energy density. They derived an upper mass limit of 2.9 $`M_{}`$ using the same maximally incompressible equation of state as Rhoades and Ruffini above twice normal nuclear energy density.

The purpose of this paper is to apply insights from relativistic kinetic theory to the problem of the maximum possible mass of neutron stars. Kinetic theory used to model viscosity and heat conduction shows requiring the sound speed to be less than the speed of light is a necessary but not sufficient condition for causality. The equation of state must satisfy a set of constraining conditions for the fluid to have stable equilibrium states with perturbations from equilbrium that propagate causally via a system of hyperbolic equations . If any one of these constraints is violated, there will be at least one fluid mode with a superluminal propagation speed. In particular, it is possible for the fluid to be acausal in a transverse mode or a different longitudinal mode even though the sound speed is below light speed. These constraints generally fix the stiffest possible equation of state to have an adiabatic sound speed significantly less than the speed of light. The main result of the present paper is application of these constraints to neutron star matter at high density softens the maximally incompressible equation of state and reduces the maximum neutron star mass significantly below the value found using a high density equation of state having the adiabatic sound speed equal to the speed of light. In the next section the thermodynamic constraint conditions from relativistic kinetic theory are reviewed. Then the maximally incompressible equation of state of neutron star matter is developed from these constraint conditions. Finally, neutron star models are constructed using this stiffest possible equation of state at high density, and the resulting maximum neutron star mass is determined. Gravitational units having $`G=c=1`$ are used.

## II relativistic kinetic theory constraints on the equation of state

The influences of dissipation in the form of bulk and shear viscosity and thermal conductivity in neutron star matter are modeled here using relativistic kinetic theory based on the Grad method of moments  as developed by Israel and Stewart . Reference  contains a concise overview of the Israel-Stewart theory. The simpler relativistic fluid theories of Eckart  and Landau-Lifshitz  are not appropriate for this problem because they have been demonstrated to exhibit pathological properties. These theories do not possess stable equilibrium states regardless of the equation of state used , and there are fluid modes having propagation speeds greater than the speed of light .

The Israel-Stewart theory has been shown to possess stable equilibrium states, and linear perturbations from equilibrium propagate casually via a hyperbolic system of equations and form a well-posed initial value problem . The thermodynamic variables describing the fluid must satisfy a set of conditions in order for the fluid to be thermodynamically stable, causal, and hyperbolic:

$$\frac{1}{\rho +p}\left(\frac{\rho }{p}\right)_s0,$$
(1)
$$\frac{1}{\rho +p}\left(\frac{\rho }{s}\right)_p\left(\frac{p}{s}\right)_\mathrm{\Theta }0,$$
(2)
$$\left(\rho +p\right)\left[1\left(\frac{p}{\rho }\right)_s\right]\frac{1}{\beta _0}\frac{2}{3\beta _2}\frac{K^2}{\mathrm{\Omega }}0,$$
(3)
$$\rho +p\frac{\beta _1+2\beta _2+2\alpha _1}{2\beta _1\beta _2\alpha _1^2}0,$$
(4)
$$\beta _00,$$
(5)
$$\mathrm{\Omega }0,$$
(6)
$$\beta _1\frac{\alpha _1^2}{2\beta _2}0,$$
(7)
$$\beta _20,$$
(8)
where

$$K=1+\frac{\alpha _0}{\beta _0}+\frac{2\alpha _1}{3\beta _2}\frac{n}{T}\left(\frac{T}{n}\right)_s,$$
(9)
$$\mathrm{\Omega }=\beta _1\frac{\alpha _0^2}{\beta _0}\frac{2\alpha _1^2}{3\beta _2}\frac{1}{nT^2}\left(\frac{T}{s}\right)_n.$$
(10)
Satisfaction of these conditions is necessary and sufficient for stability, causality, and hyperbolicity of linear perturbations. In Eqs. (1)-(10), $`\rho ,p,n,T`$, and $`s`$ are, respectively, the energy density (including mass-energy), pressure, particle density, temperature, and entropy per particle, each as measured in a comoving frame. The thermodynamic potential $`\mathrm{\Theta }`$ is the relativistic chemical potential divided by the temperature:

$$\mathrm{\Theta }=\frac{\rho +p}{nT}s.$$
(11)
$`\beta _1,\beta _2`$, and $`\beta _3`$ are effectively relaxation times for, respectively, bulk viscous, thermal, and shear viscous dissipation. $`\alpha _0`$ is the strength of the bulk stress-heat flow vector coupling, and $`\alpha _1`$ is the shear stress-heat flow vector coupling strength (see Ref.  for details). A complete fluid description requires specification of an equation of state and functional forms for the bulk and shear viscosities, thermal conductivity, relaxation times, and bulk-thermal and shear-thermal couplings. Whether or not a fluid state is thermodynamically stable, causal, and hyperbolic is determined solely by the equation of state, the relaxation times, and the viscous stress-heat flow vector coupling strengths.

For neutron star matter it is convenient to choose the particle (baryon) density $`n`$ and the temperature $`T`$ as the independent thermodynamic variables. The energy density and pressure can then be written as

$$\rho (n,T)=n\left[m+E_0(n)+I(n,T)\right]$$
(12)
and

$$p(n,T)=n^2\left[\frac{dE_0}{dn}+\left(\frac{I}{n}\right)_s\right].$$
(13)
In neutron stars the thermal energy per baryon $`I`$ is small compared to the baryon rest mass energy $`m`$ and the ground state energy per baryon $`E_0`$.

The functional forms for the relaxation times and the viscous-thermal coupling strengths are currently unknown for fully interacting nuclear matter. Their forms have been calculated for noninteracting degenerate Fermi gases :

$$\alpha _0=\frac{9\left(1+\nu \right)\left(\nu ^2+2\nu +2\right)\left(\beta \nu \right)^2}{A_0\pi ^2\nu ^2\left(\nu ^2+2\nu \right)^{3/2}}+O\left(\beta \nu \right)^0,$$
(14)
$$\alpha _1=\frac{6\left(1+\nu \right)}{A_0\left(\nu ^2+2\nu \right)^{5/2}}+O\left(\beta \nu \right)^2,$$
(15)
$$\beta _0=\frac{405\left(1+\nu \right)^5\left(\beta \nu \right)^4}{16A_0\pi ^4\nu ^4\left(\nu ^2+2\nu \right)^{1/2}}+O\left(\beta \nu \right)^2,$$
(16)
$$\beta _1=\frac{9\left(1+\nu \right)\left(\beta \nu \right)^2}{A_0\pi ^2\nu ^2\left(\nu ^2+2\nu \right)^{3/2}}+O\left(\beta \nu \right)^0,$$
(17)
$$\beta _2=\frac{15\left(1+\nu \right)}{2A_0\left(\nu ^2+2\nu \right)^{5/2}}+O\left(\beta \nu \right)^2,$$
(18)
where

$$A_0=\frac{m^4g}{2\pi ^2\mathrm{}^3}$$
(19)
and $`g`$ is the spin weight. Equations (14)-(18) are valid for any inverse dimensionless temperature ($`\beta =m/kT`$) provided the degenerate limit $`\beta \nu >>1`$ is satisfied. The dimensionless potential $`\nu `$ is the nonrelativistic chemical potential per particle divided by the baryon rest energy:

$$\nu =\frac{\rho +p}{nm}\frac{s}{k\beta }1.$$
(20)
The shear stress relaxation time Eq. (18) has recently been used to show the coupling predicted by kinetic theory between vorticity and shear viscous stress may cause colder neutron stars undergoing accertion spin-up to remain stable against $`r`$-mode instabilities at higher spin rates . Here the noninteracting degenerate Fermi gas forms will also be used to approximate the relaxation times and viscous-thermal couplings for neutron star matter. The thermal energy per baryon $`I`$ is poorly understood for neutron star matter, and the noninteracting degenerate Fermi gas form (which is proportional to $`1/\beta ^2`$) will be used.

When Eqs. (14)-(18) are substituted into Eqs. (1)-(8) and the low temperature limit $`\beta \mathrm{}`$ is taken, the result is

$$\rho 0,$$
(21)
$$p0,$$
(22)
$$\left(\rho +p\right)v_s^20,$$
(23)
$$v_s^2\frac{\rho \frac{1}{3}p}{\rho +p},$$
(24)
$$p3\rho .$$
(25)
If any of Eqs. (21)-(25) are violated, kinetic theory predicts the corresponding equilibrium state will be unstable and the fluid will possess at least one mode with a propagation speed exceeeding light speed. Both the energy density and pressure must be nonnegative for stability. Further, the adiabatic sound speed $`v_s\left[v_s^2=\left(p/\rho \right)_s\right]`$ must be real so that pressure does not decrease with increasing density: if Eq. (23) is violated, small elements of matter would be predicted to spontaneously collapse. The equation of state of neutron star matter was assumed to satisfy Eqs. (21)-(23), along with $`v_s^21`$, by Rhoades and Ruffini  and Kalogera and Baym  in deriving the maximum neutron star mass limit. Here the adiabatic sound speed will be required to satisfy Eq. (24) rather than $`v_s^21`$. Equation (24) implies the fluid state is acausal, even though the adiabatic sound speed is less than the speed of light, when

$$\frac{\rho \frac{1}{3}p}{\rho +p}<v_s^2<1.$$
(26)
Equation (25) does not affect the value of the maximum neutron star mass because it is not violated until the central density is well above the value yielding the maximum mass limit.

For neutron stars the thermal energy can be neglected and Eqs. (21)-(25) expressed solely in terms of the ground state energy per baryon and the rest energy:

$$m+E_00,$$
(27)
$$\frac{dE_0}{dn}0,$$
(28)
$$n\frac{d^2E_0}{dn^2}2\frac{dE_0}{dn}0.$$
(29)
$$m+E_0n^2\frac{d^2E_0}{dn^2}\frac{7n}{3}\frac{dE_0}{dn}0,$$
(30)
$$3m+3E_0n\frac{dE_0}{dn}0,$$
(31)
The zero temperature adiabatic sound speed can be written as

$$v_s^2=\frac{n^2\frac{d^2E_0}{dn^2}2n\frac{dE_0}{dn}}{m+E_0+n\frac{dE_0}{dn}}.$$
(32)
## III maximally incompressible equation of state

In a fully relativistic theory of neutron star matter fluid perturbations would propagate causally. One mode for propagation of fluid perturbations is via adiabatic sound waves, and many of the equations of state which are based on fitting experimental nucleon scattering data predict sound waves will travel faster than the speed of light at high density. For example, the beta-stable A18 + $`\delta v`$ \+ UIX equation of state of Akmal, Panharipande, and Ravenhall  has an adiabatic sound speed that is greater than the speed of light above $`n=0.86`$ fm-3. At zero temperature the pressure is a function of only the energy density, and so $`v_s^2=dp/d\rho `$. The pressure must not increase with increasing energy density faster than

$$p(\rho )=\rho +\text{const}$$
(33)
if the adiabatic sound speed is to not exceed light speed. Using the expressions for the energy density and pressure written as functions of the ground state energy per baryon \[Eqs. (12) and (13)\], Eq. (33) at zero temperature can be written as a first-order ordinary differential equation for $`E_0`$:

$$n^2\frac{dE_0}{dn}=n\left(m+E_0\right)+\text{const}.$$
(34)
The general solution of Eq. (34) is

$$E_0(n)=C_1\left(\frac{n}{\stackrel{~}{n}}\right)+C_2\left(\frac{n}{\stackrel{~}{n}}\right)^1m.$$
(35)
Equation (35) can be regarded as the maximally incompressible equation of state consistent with a subluminal sound speed.

Relativistic kinetic theory has shown an adiabatic sound speed less than the speed of light is a necessary but not sufficient condition for the causal propagation of fluid perturbations. The equation of state must satisfy Eqs. (21)-(25) for all densities relevant for neutron star models for the equation of state to be causal. In this paper the maximally incompressible equation of state is taken as the form that satisfies Eq. (24) with an equality:

$$v_s^2=\frac{\rho \frac{1}{3}p}{\rho +p}.$$
(36)
Equations (12), (13), and (32) can be used to write Eq. (36) as a second-order ordinary differential equation for $`E_0`$:

$$n^2\frac{d^2E_0}{dn^2}+\frac{7}{3}n\frac{dE_0}{dn}E_0=m.$$
(37)
The corresponding differential equation for the total energy per baryon $`E=m_0+E_0`$ is homogeneous:

$$n^2\frac{d^2E}{dn^2}+\frac{7}{3}n\frac{dE}{dn}E=0.$$
(38)
This is a form of the Cauchy linear equation. The substitution $`n=e^z`$ results in

$$\frac{d^2E}{dz^2}+\frac{4}{3}\frac{dE}{dz}E=0,$$
(39)
which has constant coefficients and can be solved for $`E(z)`$ with standard methods. After the substitution $`z=ln(n)`$ the general solution to Eq. (37) results in

$`E_0(n)`$ $`=D_1\left({\displaystyle \frac{n}{\stackrel{~}{n}}}\right)^{\left(\sqrt{13}2\right)/3}+D_2\left({\displaystyle \frac{n}{\stackrel{~}{n}}}\right)^{\left(\sqrt{13}+2\right)/3}m`$ (41)
$`D_1\left({\displaystyle \frac{n}{\stackrel{~}{n}}}\right)^{0.535}+D_2\left({\displaystyle \frac{n}{\stackrel{~}{n}}}\right)^{1.869}m.`$
Equation (41) can be regarded as the maximally incompressible equation of state compatible with the causality condition Eq. (24). Equation (41) predicts the ground state energy increases more slowly with increasing baryon density than for Eq. (35), so the pressure, which is proportional to the rate of change of the ground state energy, will increase more slowly. Hence Eq. (41) is a softer maximally incompressible equation of state than Eq. (35) and will less readily support the star against gravitational collapse, hence yielding a lower maximum neutron star mass limit than Eq. (35).

## IV maximum mass of neutron stars

In this section the kinetic theory maximally incompressible equation of state Eq. (41) is used to determine a new least upper bound on the maximum neutron star mass stable against gravitational collapse to a black hole. An experimentally verified equation of state is used at low density, and the maximally incompressible equation of state Eq. (41) is used at high density and compared to using the sound speed equals light speed maximally incompressible equation of state Eq. (35). The low density equation of state consists of the BPS equation of state  below normal nuclear density, as in Ref. , and one of the recently developed equations of state which utilize fits of nucleon scattering data and light nuclei properties to cover the density range above normal nuclear density. Three of these experimentally based equations of state will be used here as representative examples: the beta-stable AV14 + UVII and the beta-stable UV14 + UVII models of Ref.  (hereafter WFF-A and WFF-U), and the beta-stable A18 + $`\delta v`$ \+ UIX\* model of Ref.  (hereafter APR). The maximum mass value is insensitive to the particular equation of state applied below normal nuclear density because little of the stellar mass is in the subnuclear crust when the star is near the mass limit. The low and high density equations of state are matched at a baryon density $`\stackrel{~}{n}`$ by choosing the integration constants $`C_i`$ and $`D_i`$ in Eqs. (35) and (41) so that the ground state energy and pressure values of the low and high density equations of state have the same value at the match point. For example, using the WFF-A equation of state at low density and choosing the match point to be twice normal nuclear energy density (304 MeV/fm3, 0.317 baryons/fm3), the ground state energy and pressure at the match point are $`E_0(\stackrel{~}{n})=21.9`$ MeV and $`p(\stackrel{~}{n})=8.67`$ MeV/fm3. The corresponding values of the $`C_i`$ and $`D_i`$ in the maximally incompressible equations of state are

$$C_1=\frac{\rho (\stackrel{~}{n})+p(\stackrel{~}{n})}{2\stackrel{~}{n}}494\text{ MeV},$$
(42)
$$C_2=\frac{\rho (\stackrel{~}{n})p(\stackrel{~}{n})}{2\stackrel{~}{n}}467\text{ MeV},$$
(43)
$$D_1=\frac{\left(2+\sqrt{13}\right)\rho (\stackrel{~}{n})+3p(\stackrel{~}{n})}{2\sqrt{13}\stackrel{~}{n}}759\text{ MeV},$$
(44)
$$D_2=\frac{\left(2\sqrt{13}\right)\rho (\stackrel{~}{n})+3p(\stackrel{~}{n})}{2\sqrt{13}\stackrel{~}{n}}203\text{ MeV}.$$
(45)
The Oppenheimer-Volkoff equation describes the interior of a static star in general relativity:

$$\frac{dp}{dr}=\frac{\left(\rho +p\right)\left[m(r)+4\pi r^3p\right]}{r\left[r2m(r)\right]},$$
(46)
where

$$m(r)=_0^r4\pi \rho \left(\xi \right)\xi ^2𝑑\xi .$$
(47)
The calculation of the maxmimum neutron star mass limit begins by choosing a central energy density value $`\rho _c`$ (or, equivalently, a central baryon density $`n_c`$) and integrating Eq. (46) radially outward to the edge of the star, marked by the radius $`R`$ where the pressure drops to zero. One of the maximally incompressible equations of state are used in the high density interior out to the radius where the energy density drops below the matching energy density $`\stackrel{~}{\rho }=\rho \left(\stackrel{~}{n}\right)`$. The experimentally verified low density equation of state is then used out to the edge of the star, and the mass of the star $`M(R)`$ determined from Eq. (47). The integration of Eq. (46) is repeated for a series of progressively increasing central densities. The maximum mass value is determined by the central density for which $`M\left(\rho _c\right)`$ has its stationary point ($`dM/d\rho _c=0`$). Above this central density $`dM/d\rho _c<0`$ and the star will be unstable to collapse to a black hole from radial perturbations (see, for example, Refs.  or ).

Shown in Fig. 1 is the neutron star mass as a function of the central baryon density, with the BPS equation of state used below normal nuclear density, the WFF-A equation of state used above normal nuclear density and below the matching energy density, and the kinetic theory maximally incompressible equation of state Eq. (41) used above the matching energy density. Mass models having a matching density corresponding to one, two, and four times normal nuclear energy density $`\rho _{nn}`$ are shown. Also shown for comparison are models with the sound speed equal to light speed maximally incompressible equation of state Eq. (35) of Refs.  and  used above the matching density, and models with the WFF-A equation of state used at all densities above normal nuclear density. Tables I-III list some specific numerical values of the maximum mass for a few different matching densities for the WFF-A, WFF-U, and APR equations of state. The softer kinetic theory maximally incompressible equation of state is seen to yield smaller maximum mass values than the sound speed equal to light speed maximally incompressible equation of state. If the WFF-A equation of state is regarded as experimentally verified up to twice normal nuclear energy density, as was assumed in Ref. , the result of using the kinetic theory maximally incompressible equation of state instead of the sound speed equal to light speed maximally incompressible equation of state is to reduce the maximum neutron star mass from 2.92 $`M_{}`$ to 2.63 $`M_{}`$. The average maximum mass reduction for these three low density equations of state, when using a matching energy density of twice normal nuclear density, is 0.29 $`M_{}`$, from an average of 2.93 $`M_{}`$ to 2.64 $`M_{}`$.

Figure 2 shows the maximum neutron star mass as a function of the matching energy density for the three low density equations of state and the two maximally incompressible high density equations of state. The three low density equations of state yield similar maximum mass values for matching densities up to twice normal nuclear energy density, so the value of the maximum mass derived from the kinetic theory maximally incompressible equation of state is substantially independent of the low density model for matching densities in this density range. Using the WFF-A, WFF-U, or APR equations of state for all densities above normal nuclear density yields a maximum mass of, respectively, 2.13, 2.19, or 2.21 $`M_{}`$. Examination of Fig. 2 shows when the matching energy density approaches four times normal nuclear energy density, the stellar mass value is largely determined by the low density equation of state. As was noted by Kalogera and Baym , once nuclear physics advances to where the experimentally based equation of state is verified up to four times normal nuclear density, the maximum neutron star mass problem will be solved.

## V discussion

The principal finding of this paper is stability and causality constraints resulting from applying kinetic theory to thermal and viscous dissipative processes in neutron star matter yields a softer maximally incompressible equation of state, and a lower maximum mass of neutron stars, than for the sound speed equal to light speed maximally incompressible equation of state. The reduction in the maximum mass averages from 2.93 $`M_{}`$ down to 2.64 $`M_{}`$ for the three low density equations of state considered if the low density equation of state is regarded as experimentally verified up to twice normal nuclear energy density.

The softening of the maximally incompressible equation of state is a generic prediction of kinetic theory, but the specific value of the resulting maximum neutron star mass derived here rests upon the accuracy of approximating the relaxation times and viscous-thermal couplings with their noninteracting degenerate Fermi gas forms, and the overall validity of dissipative relativistic kinetic theory for strongly interacting nuclear matter at high density. Future data from the Relativistic Heavy Ion Collider holds promise for experimental testing of the applicability of the kinetic theory approach to modeling high density nuclear matter.

Only static neutron star mass models have been constructed in this paper. Rotation can support larger mass stars against collapse, but only for stars rotating rapidly enough to be nearly shedding mass from the equator is the maximum mass value significantly changed from its nonrotating value. For these rapidly rotating stars the increase in the maximum mass is at most approximately 20% over the static value when undergoing uniform rotation . Differential rotation is strongly damped . The remnant of binary neutron star coalescence may undergo a short period of dynamically stable differential rotation, and thus briefly support a much larger mass neutron star remnant against prompt black hole collapse than is possible for a static or uniformly rotating star .

###### Acknowledgements.

This work was supported by the U.S. National Aeronautics and Space Administration through Grant No. NAG5-9541.