text:
# Lattice Dynamics and the High Pressure Equation of State of Au

## I Introduction

The elastic constants, phonon frequencies, and equation of state (EOS) are fundamental properties of matter. The values of these parameters under compression find application in geophysics manghnani  and in the prediction and interpretation of processes of dynamic compression. mcqueen ; bgp  For many materials, especially elemental metals, the principal Hugoniot curves, the set of states accessible via a single shock from ambient conditions, have been measured. mcqueen ; llnl ; lasl ; altshuler  Since the pressure, density, and internal energy are known along the Hugoniot from the jump conditions,mcqueen  this data is an important baseline for high pressure equations of state. The off-Hugoniot EOS is needed for the prediction of processes involving more complicated loading paths, such as multiple shocks or shock and release. nellis  Extrapolation from the Hugoniot has also been used to establish pressure standards for static high pressure experiments, whose room temperature isotherms are regarded as known. mcqueen ; jamieson ; mbss  Au has been used as a standard in this way, and has been used to calibrate the ruby $`R_1`$ line as a secondary standard.mbss ; mxb ; bell\_sccm  Recent studies have called into question the accuracy of the gold pressure standard. Akahama et al.akahama  compressed Au and Pt in the same cell and found two Au standardsjamieson ; heinz  to give pressures lower than Pt by 20 and 15 GPA, respectively, at 150 GPa. Shim et al.,shim  on the other hand, propose a new EOS for Au that gives pressures still lower than either of these earlier standards.

To relate the Hugoniot to the room temperature isotherm requires information on the Grüneisen parameter $`\gamma =V\left(P/E\right)_V`$ and the specific heat $`C_V`$.mcqueen  These are dominated by lattice vibrations in the regimes under consideration here. With increasing compression, the separation between the Hugoniot pressure and the room temperature pressure increases, resulting in greater dependence of the inferred room temperature isotherm on $`\gamma `$. Because $`\gamma `$ is not easily measured at high pressure, this introduces a non-negligible source of uncertainty in the pressure standards. In most cases, a model assumption of the form $`\gamma (V)=\gamma (V_0)(V/V_0)^q`$ has been used. The specific value $`q=1`$ has been used often in shock work.mcqueen ; jamieson  However, a power law dependence of $`\gamma `$ on $`V`$ is qualitatively incorrect at high compression, and extrapolation on this basis is inherently limited.

In principle, the various components of the EOS – the static lattice energy, the vibrational free energy, and the electronic excitation free energy – can be evaluated from electronic structure theory. Practical calculations based on approximate density functional theories typically have errors of several percent in the density at zero pressure, and 10 % or more in the bulk modulus. khein  These are unacceptably large errors for the purpose of high accuracy equations of state. These properties can be measured accurately, and are mainly determined by the cold energy curve, which is the largest contribution to the EOS in the regime considered here. On the other hand, we find that ab initio electronic structure calculations are capable of obtaining phonon frequencies of sufficient accuracy to strongly constrain the volume dependence of the Grüneisen parameter $`\gamma `$. We therefore propose that the most accurate EOS in the solid is obtained by combining an empirical cold energy with ab initio lattice vibration and electronic excitation free energies. This is analogous to procedures that have been used to reduce shock Hugoniot data and derive room temperature standards, but the present analysis gives a strong physical foundation for the volume dependence of $`\gamma `$, allowing confidence in the results at higher compression. We are not aware of any evidence of solid-solid phase transitions in Au, and we consider only the fcc solid. Recent calculationsboettger\_au  indicate a transition to the hcp structure for $`V/V_00.6`$, and we do not attempt to extend our solid EOS beyond this density.

This paper is arranged as follows. First we discuss lattice dynamics and its connection to the EOS. We emphasize the importance of the classical limit for defining the ion motion contribution to $`\gamma `$. We describe our procedures for calculating phonon frequencies and interpolating to the whole Brillouin zone. Next we present our results for elastic and phonon properties of Au, and our analysis to obtain $`\gamma ^{\mathrm{ion}}(V)`$. We then present the complete EOS for solid Au by giving a set of closed form expressions for the Helmholtz free energy with parameters. The resulting room temperature isotherm is compared with various proposed standards. Next we discuss shock melting and its effect on the Hugoniot within an approximate model, and finally give our conclusions.

## II Lattice Dynamics and the Equation of State

We write the Helmholtz free energy as

$$F(V,T)=\varphi _0(V)+F^{\mathrm{ion}}(V,T)+F^{\mathrm{el}}(V,T)$$
(1)
where $`\varphi _0`$ is the static lattice energy, $`F^{\mathrm{ion}}`$ is the ion motion free energy, and $`F^{\mathrm{el}}`$ is the electronic excitation free energy. In our applications, $`\varphi _0`$ gives the largest contribution to the pressure. The ion motion free energy gives the dominant temperature dependence of $`P=\left(F/V\right)_T`$. The electronic excitation term $`F^{\mathrm{el}}`$ is generally a small correction to the solid EOS, becoming important for the liquid Hugoniot at several hundred GPa.

It is, therefore, imperative to have an accurate $`F^{\mathrm{ion}}`$ for calculating the temperature dependence of $`P`$. The quasi-harmonic approximation has been found to have small errors in many cases.wallace\_therm  We have carried out Monte Carlo simulations using an embedded atom modelvoter  of Cu to investigate the importance of anharmonic corrections to the ion free energy. Here we define the anharmonic free energy to be the difference between the true ion free energy and the quasi-harmonic approximation, noting that the term has been used differentlyanderson\_au  by other authors.footnote01  We intend to publish details of our Cu simulations elsewhere. In summary, we find that up to 300 GPa on the melting curve, $`|P^{\mathrm{anh}}|<0.85`$ GPa, and is never more than 3% of the thermal pressure $`P(V,T)P(V,0)`$, a very small correction to the total pressure. Given the similarity of the bonding in Cu and Au, we expect these results to be relevant for Au also. Thus, in what follows we neglect anharmonicity and use the quasi-harmonic approximation for $`F^{\mathrm{ion}}`$,

$$F^{\mathrm{ion}}(V,T)=_0^{\mathrm{}}𝑑\omega g(\omega )\left[\frac{1}{2}\mathrm{}\omega +k_BT\mathrm{ln}(1e^{\mathrm{}\omega /k_BT})\right],$$
(2)
where we have introduced $`g(\omega )=\frac{1}{3N}_𝐤\delta (\omega \omega _𝐤)`$ , the normalized phonon density of states. The phonon frequencies $`\omega _𝐤`$ are functions of volume only, and the sum is over the $`3N`$ normal modes of the crystal.

Application of Eq. (2) requires the full phonon density of states $`g(\omega )`$ at all volumes. While in principle this information is available from our calculations, in practice the classical limit dominates our EOS, allowing for a substantial simplification. In the classical limit, which is the leading term in the high temperature expansion of Eq. (2), the free energy is given by

$$F^{\mathrm{cl}}(V,T)=3Nk_BT\mathrm{ln}\left(\frac{e^{1/3}\mathrm{}\omega _0}{k_BT}\right).$$
(3)
where we have introduced the moment

$$\omega _0=e^{1/3}\mathrm{exp}\left[_0^{\mathrm{}}𝑑\omega g(\omega )\mathrm{ln}\omega \right].$$
(4)
Other moments are conventionally defined as

$$\omega _n=\left[\frac{3+n}{3}_0^{\mathrm{}}𝑑\omega g(\omega )\omega ^n\right]^{1/n}n0,n>3,$$
(5)
where the normalizations in Eqs. (4) and (5) are chosen so that for a Debye spectrum, described by

$$g(\omega )=\frac{3}{\omega _D^3}\omega ^2\mathrm{\Theta }(\omega _D\omega ),$$
(6)
all $`\omega _n`$ are equal to $`\omega _D`$, the Debye frequency. We subsequently give results for $`\nu _n=\omega _n/(2\pi )`$ corresponding to the experimental convention of giving frequencies $`\nu `$ in Hertz, as opposed to angular frequencies, $`\omega `$.

In the classical limit, the ion pressure is linear in $`T`$ with $`\left(P^{\mathrm{ion}}/T\right)_V=\left(3Nk_B/V\right)\gamma _0`$, where $`\gamma _0=\frac{d\mathrm{ln}\omega _0}{d\mathrm{ln}V}`$, and the specific heat is constant, $`C_V=3Nk_B`$. The role of quantum ion motion in the EOS is illustrated in Figure 1, which shows the pressure along an isochore $`V/V_0=0.8`$. The solid curve is our full EOS, described below, and the dashed curve uses the classical limit of $`F^{\mathrm{ion}}`$ at all $`T`$. The arrows mark room temperature and the Hugoniot temperature, which is 1340 K at the given density. The melting temperature at this density is estimated to be 4900 K, beyond the range of the plot. It is clear that the classical limit dominates even at room temperature. The Hugoniot is well into the classical regime. At higher densities, the Hugoniot is still further above the quantum regime. The largest quantum effect on the pressure is at $`T=0`$, where the zero-point vibrations contribute a pressure of 0.8 GPa at this density. The classical limit is clearly dominant for interpolating between the Hugoniot and room temperature.

We have found that an interpolation based on the Debye model gives a very accurate approximation to the full quasi-harmonic free energy. The Debye free energy is a special case of the quasi-harmonic free energy, Eq. (2). Inserting the Debye density of states, Eq. (6), gives

$$F^\mathrm{D}(V,T)=N\left[\frac{9}{8}\mathrm{}\omega _D+3k_BT\mathrm{ln}\left(1e^{\mathrm{}\omega _D/k_BT}\right)k_BTD\left(\mathrm{}\omega _D/k_BT\right)\right],$$
(7)
where

$$D(x)=\frac{3}{x^3}_0^x𝑑z\frac{z^3}{e^z1}.$$
(8)
In light of the above remarks, it is important to capture the correct classical limit. This requires that we set

$$\omega _D(V)=\omega _0(V).$$
(9)
To emphasize the distinction between Eq. (9) and the standard definition of $`\omega _D`$ in terms of acoustic velocities, we refer to the Debye free energy together with Eq. (9) as the high-temperature Debye model. The high-temperature Debye model gives the same results as the full quasi-harmonic free energy in the classical regime, and in addition obeys the Nernst theorem at low $`T`$. We have checked that at low $`T`$, the error in the pressure compared to the full quasi-harmonic approximation is entirely negligible.

Thus, we can simplify the specification of the lattice vibration free energy from giving $`g(\omega )`$ at all $`V`$ to giving the single moment $`\omega _0`$ at all $`V`$. This is advantageous for numerical applications. It also allows us to express our full EOS in a few compact formulae so that it is generally accessible. These formulae are given below.

Calculating the moments $`\omega _n`$ requires the phonon frequencies for all $`𝐤`$ in the Brillouin zone. Direct ab initio calculations on a dense mesh in $`𝐤`$ would be quite expensive. As a result of another studygraf\_pdf  we found that the low order moments can often be accurately obtained from short-ranged force constant models. In particular, for Au, the moment $`\omega _0`$ is converged to less than 1 % with a 2nd neighbor (2NN) model. Therefore, our method is to calculate four zone boundary phonon frequencies corresponding to the transverse and longitudinal modes at the $`X`$ and $`L`$ points. These are calculated with standard frozen-phonon methods. In addition the three elastic moduli are calculated using the method described by Söderlind et al.soderlind\_elastic  We fit these results to a 2NN force constant model, which then allows the evaluation of $`\omega _𝐤`$ for arbitrary $`𝐤`$ for integration over the Brillouin zone.

The electronic structure calculations used the full-potential LAPW code WIEN97.Wien97  We used the LDA rather than the GGA, based on Boettger’s findingboettger\_au  that the LDA gives a better static lattice energy than the GGA for Au. Some numerical parameters used in the calculations were, in atomic units: muffin tin radius $`r_{\mathrm{MT}}=2.0`$; plane wave cut-off $`r_{\mathrm{MT}}k_{\mathrm{max}}=9.0`$; cut-off for expansion of density and potential $`g_{\mathrm{max}}=14`$. For elastic modulus calculations, Brillouin zone integrals used special points corresponding to $`18^3`$ points in the full zone, with Gaussian smearing of the energies by 20 mRy. The zone boundary phonons were found to be comparatively insensitive to the $`𝐤`$-point mesh, and smaller meshes of $`10^3`$ points were used. The $`5p`$, $`5d`$, and $`6s`$ shells were treated as valence states, and local orbital extensionssingh\_lo  were used in the $`p`$ and $`d`$ channels.

We used a generalized Born-von Kármán force model for our lattice-dynamical calculations, from which we were able to compute the phonon dispersions in the entire Brillouin zone. Since gold has a very simple phonon dispersion, we employed only a model with first (1NN) and second (2NN) nearest-neighbor interatomic shells of atoms. In an fcc lattice the 1NN and 2NN forces are determined by three and two independent parameters, respectively (for more details see e.g. ref. graf\_pdf, ). These five force constants were extracted by fitting simultaneously the phonon frequencies of gold at the $`X`$ and $`L`$ points of the Brillouin zone and the elastic constants C11, C12, and C44 near the zone center ($`\mathrm{\Gamma }`$ point). Thus we fit a total of seven independent data points. We gave equal weight to the zone boundary phonons and elastic constants in our $`\chi ^2`$-fit of the phonon dispersions. Although the 2NN Born-von Kármán force model is too simple to reproduce all phonon frequencies within less than approximately 10%, see Figure 2, it is sufficiently accurate to compute integrated quantities such as the phonon moments within less than approximately 3%. More accurate phonon dispersions can be obtained if needed, by computing ab initio frequencies at half and quarter distances in the Brillouin zone and fitting those to a 3NN or higher-order Born-von Kármán force model, or to phonon models with interatomic pseudo-potentials.

## III Elastic and Phonon Properties

Figure 2 shows the calculated phonon dispersion curves for Au at the density corresponding to ambient pressure and room temperature. The filled diamonds at the $`X`$ and $`L`$ points are the LDA frozen phonon results. These, together with the calculated elastic moduli are used to obtain the force constants. The solid curves are the interpolation to general wavevectors using the force constant model. The open circles are the experimental data of Lynn et al.Lynn1973  The force constant model is over-constrained, so it neither goes precisely through the calculated zone boundary frequencies, nor has it exactly the calculated elastic moduli. The shapes of the dispersion curves are simple enough for Au that they are generally well captured by the 2NN model.

Table 1 summarizes our results for the zone boundary phonons, elastic moduli, and moments $`\nu _n=\omega _n/2\pi `$, as functions of volume. The reference volume $`V_0`$ corresponds to $`T=298`$ K and $`P=1`$ bar $`=100`$ kPa. For Au, $`V_0=`$ 10.212 cm3/mol or 114.43 $`a_0^3`$/atom. Also shown are the experimental data at room temperature.Lynn1973 ; simmons  No attempt is made to account for the temperature dependence beyond comparing at the correct volume. There is generally good agreement between the calculated and experimental quantities. The main exception is $`C_{44}`$, which is calculated substantially lower than measured. Our calculation is in better agreement with an earlier LDA calculation.soderlind\_elastic  The main result of the present calculation is the value of $`\nu _0`$, which is within 3% of the measurement. The ratio $`\nu _1/\nu _0`$ is 1 for a Debye spectrum. Our calculations give $`\nu _1/\nu _0=1.03`$ at ambient density and 1.04 at twofold compression, so, by this measure, the departure from a Debye spectrum is small and nearly constant with volume.

Measurements of elastic moduli were reported by Duffy et al.duffy\_cij  to $`P=37`$ GPa, and ab initio calculations by Tsuchiya and Kawamuratsuchiya  to 100 GPa. Figure 3 shows these results along with our calculations. Our volumes are converted to room temperature pressures using the EOS described below. The highest density in table 1 corresponds to $`P780`$ GPa, and the graph is restricted to lower pressures to highlight the comparison with these other works. Our calculations of $`C_{11}`$ and $`C_{12}`$ are in good agreement with Tsuchiya and Kawamura, while our $`C_{44}`$ is systematically lower, in somewhat better agreement with the experiments. Our $`C_{11}`$ is in good agreement with the experiments while $`C_{12}`$ is slightly high. It should be noted that the experimental $`C_{ij}`$ depend on a model parameter $`\alpha `$ used in the analysis.duffy\_cij  The plotted points correspond to $`\alpha =1`$. It is encouraging that our calculations give the correct trends for the pressure dependence of $`C_{ij}`$.

Having determined the moments $`\nu _n`$ as functions of volume, we need to interpolate and differentiate them to obtain thermodynamic functions. To do this, we assumed the following functional form for $`\gamma _n(V)`$, which has been used in creating many wide ranging equations of state,grizz

$$\gamma _n(V)=\gamma ^{\mathrm{}}+A_n\left(\frac{V}{V_0}\right)+B_n\left(\frac{V}{V_0}\right)^2$$
(10)
where $`\gamma ^{\mathrm{}}`$ is the infinite density limit of $`\gamma `$, and $`A_n`$ and $`B_n`$ are parameters. Alternatively, we can express $`A_n`$ and $`B_n`$ in terms of $`q_n=d\mathrm{ln}(\gamma _n)/d\mathrm{ln}(V)`$ and $`\gamma _n(V_0)`$ as

$`A_n`$ $`=`$ $`\gamma _n(V_0)\left[2q_n(V_0)\right]2\gamma ^{\mathrm{}}`$
$`B_n`$ $`=`$ $`\gamma _n(V_0)\left[q_n(V_0)1\right]+\gamma ^{\mathrm{}}`$ (11)
Integrating $`\gamma _n=\frac{d\mathrm{ln}\nu _n}{d\mathrm{ln}V}`$, we have

$$\nu _n(V)=\nu _n(V_0)\left(\frac{V}{V_0}\right)^\gamma ^{\mathrm{}}\mathrm{exp}\left\{A_n\left[\left(\frac{V}{V_0}\right)1\right]\frac{B_n}{2}\left[\left(\frac{V}{V_0}\right)^21\right]\right\}.$$
(12)
By fitting this functional form to the calculated $`\nu _0(V)`$, we extract the parameters giving $`\gamma _0(V)`$. In our fits, we keep $`\gamma ^{\mathrm{}}`$ fixed. The value $`\gamma ^{\mathrm{}}=2/3`$ has been widely used,bushman  although arguments have been made for $`\gamma ^{\mathrm{}}=1/2`$.lb\_gamma  For our application and density range, we have found that the quality of the fit and the resulting $`\gamma _0(V)`$ are insensitive to this choice, and we use $`\gamma ^{\mathrm{}}=2/3`$. The fitting procedure is illustrated in Figure 4, where we show the fit for $`\nu _0(V)`$, the log moment and the resulting $`\gamma _0(V)`$. We did the fits with $`\nu _0(V_0)`$, $`\gamma _0(V_0)`$, and $`q_0(V_0)`$ as free parameters, and with $`\gamma _0(V_0)`$ constrained to the experimental value of 2.97. (This experimental value was determined so that the thermal expansion for the subsequent EOS overlies the recommended curve of Touloukiantouloukian  from 100 to 1200 K.) The unconstrained fit gives $`\gamma _0(V_0)=2.88`$. Constraining $`\gamma _0(V_0)`$ increases the RMS error of the fit from $`4.9\times 10^2`$ to $`5.7\times 10^2`$ THz, which is not a significant increase. Also, we note that constraining $`\gamma _0(V_0)`$ does not change $`\gamma _0(V)`$ at smaller $`V`$. For this reason, we believe that the constrained fit is the best overall approximation for $`\gamma ^{\mathrm{ion}}(V)`$ up to twofold compression. This is shown as the solid curve in Figure 4. The corresponding parameters are $`\nu _0(V_0)=3.46`$ THz, $`\gamma _0(V_0)=2.97`$, $`q_0(V_0)=1.37`$. We do not show the unconstrained fit in Figure 4 because it is too close to the constrained fit. The dot-dashed curves correspond to $`q=\mathrm{const}`$, $`\gamma (V)=\gamma (V_0)(V/V_0)^q`$. The value $`q=1`$ gives a $`\gamma `$ that is too high over this density range, while $`q=1.7`$heinz  is too low. Jamieson et al.jamieson  used $`q=1`$ and an ambient value $`\gamma (V_0)=3.215`$, which results in the dashed curve. We believe that their $`\gamma `$ is too high at all volumes considered here. Any constant value of $`q`$ will lead to a $`\gamma `$ that is too small at very high compression.

## IV Equation of State

The complete EOS is determined when the Helmholtz free energy $`F`$ is given as a function of volume and temperature. In this section we describe our form for $`F`$ and give the numerical parameters for fcc Au. First we write the free energy as in Eq. (1) as the sum of the static lattice energy $`\varphi _0`$, an ion motion free energy $`F^{\mathrm{ion}}`$, and an electronic excitation term $`F^{\mathrm{el}}`$. For the ion motion free energy we use the quasi-harmonic approximation, and further specialize to the high temperature Debye model. Thus we take the ion free energy to be given by Eqs. (7) and (8), where we identify $`\omega _D`$ with the log moment $`\omega _0`$ of the phonon frequencies. As discussed in section II, this gives the correct classical limit, and leads to extremely small errors at low $`T`$ compared to the full phonon spectrum. This way we are able to give a closed form expression for $`F^{\mathrm{ion}}`$, which is very convenient in numerical applications. The volume dependence of $`\omega _0`$ is that given in Eq. (12) with parameters given above. For the final EOS we have replaced the fitted value of $`\omega _0(V_0)`$ with the experimental one. This is essentially a small shift in the absolute entropy, which has negligible effect on the results.

Electronic excitations give a small contribution to the thermodynamics, which we approximate by

$$F^{\mathrm{el}}(V,T)=\frac{1}{2}N\mathrm{\Gamma }(V)T^2.$$
(13)
The Sommerfeld coefficient $`\mathrm{\Gamma }`$ is proportional to the electronic density of states at the Fermi energy, $`\mathrm{\Gamma }=(\pi ^2/3)k_B^2g(ϵ_f)`$. Our calculations show that $`d\mathrm{ln}g(ϵ_f)/d\mathrm{ln}V0.76`$ over the range of densities considered, so we take

$$\mathrm{\Gamma }(V)=\mathrm{\Gamma }(V_0)\left(V/V_0\right)^\kappa $$
(14)
with $`\kappa =0.76`$. The calculations give $`\mathrm{\Gamma }(V_0)=6.7\times 10^4`$ J/mol K2, while the measured value from the low-temperature specific heat ishultgren  $`\mathrm{\Gamma }(V_0)=7.28\times 10^4`$ J/mol K2. The measured low-temperature specific heat coefficient $`\mathrm{\Gamma }_{\mathrm{expt}}`$ is expected to be enhanced with respect to the bare value $`\mathrm{\Gamma }`$ as a result of electron-phonon interactions by a factor $`(1+\lambda )`$, where $`\lambda `$ is the dimensionless electron-phonon mass enhancement parameter. The value for gold is $`\lambda 0.050.15`$,allen  which agrees with our calculated $`\mathrm{\Gamma }/\mathrm{\Gamma }_{\mathrm{expt}}=1.09`$. At high temperatures, $`k_BT>\mathrm{}\omega _0`$, the phonon mass enhancement becomes ineffective. Hence we use our calculated density of states in the remainder of this paper.

The static lattice energy is needed to complete the free energy. We adopt the functional form due to Vinet et al.vinet

$`\varphi _0(V)`$ $`=`$ $`{\displaystyle \frac{4V^{}B^{}}{(B_1^{}1)^2}}\left[1(1+X)e^X\right]`$
$`X`$ $`=`$ $`{\displaystyle \frac{3}{2}}(B_1^{}1)\left[\left({\displaystyle \frac{V}{V^{}}}\right)^{1/3}1\right]`$ (15)
which is parameterized by the volume at the minimum, $`V^{}`$, the bulk modulus, $`B^{}`$ and its pressure derivative, $`B_1^{}`$. These parameters have been determined empirically by requiring that the EOS reproduce the ambient volume, and ultrasonic data for $`B_S`$ at ambient conditions.simmons  For the latter, we adopt the value $`B_S=173`$ GPa. Ultrasonic data gives values for $`\left(B_S/P\right)_T`$ from 5.2 to 6.4. Therefore, to complete the EOS, we require that it reproduce the measured slope of the Hugoniot, which we take from the fitllnl  $`U_s=3.12\mathrm{km}/\mathrm{s}+1.521U_p`$ relating the shock velocity $`U_s`$ to the particle velocity $`U_p`$. The resulting parameters give $`\left(B_S/P\right)_T=5.49`$, consistent with ultrasonic data.

The complete set of parameters for the free energy of fcc Au are summarized here. The static lattice energy is given by Eq. (15) with

$`V^{}`$ $`=`$ $`10.0834\mathrm{cm}^3/\mathrm{mol}`$
$`B^{}`$ $`=`$ $`180.0\mathrm{GPa}`$
$`B_1^{}`$ $`=`$ $`5.55.`$ (16)
The ion motion free energy is given by Eqs. (7) and (8) with the volume dependence of $`\omega _D`$ given by Eq. (12) (with $`\nu _0\omega _D`$, $`\gamma _0\gamma `$, etc.). The parameters $`A`$ and $`B`$ are given in terms of $`\gamma (V_0)`$ and $`q(V_0)`$ by Eq. (11). The numerical values are

$`V_0`$ $`=`$ $`10.212\mathrm{cm}^3/\mathrm{mol}`$
$`\omega _D(V_0)`$ $`=`$ $`22.9\times 10^{12}\mathrm{s}^1`$
$`\gamma (V_0)`$ $`=`$ $`2.97`$
$`q(V_0)`$ $`=`$ $`1.3677`$
$`\gamma ^{\mathrm{}}`$ $`=`$ $`2/3.`$ (17)
Finally, the electronic excitation free energy is given by Eqs. (13) and (14) with

$`\mathrm{\Gamma }(V_0)`$ $`=`$ $`6.7\times 10^4\mathrm{J}/\mathrm{mol}\mathrm{K}^2`$
$`\kappa `$ $`=`$ $`0.76.`$ (18)
The following calculations use these parameters for the free energy of fcc Au. Once the free energy is known, the pressure, internal energy, etc., can be evaluated. We calculate the Hugoniot by fixing a value of the volume and solving for the temperature such that the jump condition, $`EE_0=\frac{1}{2}(P+P_0)(V_0V)`$ , is solved. Here $`E_0`$, $`V_0`$ and $`P_0`$ correspond to the initial state, taken to be ambient temperature and pressure.

Figure 5 shows the Hugoniot and room temperature isotherm for the present EOS. The solid symbols are the Hugoniot data.lasl ; altshuler ; jones  The dashed line corresponds to the linear fit for $`U_s(U_p)`$. Also shown as the open circles are the room temperature data from Heinz and Jeanlozheinz  taken to 70 GPa with the ruby pressure scale. The open squares are the room temperature data of Bell et al.bell\_sccm  to 188 GPa. The pressure shown in the graph is taken from the ruby scale, which is extrapolated from a lower pressure calibration.mbss  Our calculations agree with these data to 3% at the highest pressure, indicating support for the extrapolation of the ruby scale.

For $`V>0.656V_0`$, the Hugoniot lies in the solid. By matching the measured Hugoniot with our theoretically based Grüneisen parameter, we confirm the accuracy of our room temperature isotherm, and we have high confidence in our EOS to this density. Combining our calculated Hugoniot temperatures with a Lindemann melting curve, we estimate that the Hugoniot enters a solid-liquid coexistence region at $`V/V_0=0.656`$, $`P=280`$ GPa. A more extensive discussion of shock melting is given in the following section. There we find that explicitly accounting for melting leads to good agreement with the high pressure Hugoniot data, indicating that our solid EOS is valid to densities of $`V/V_00.6`$.

Figure 6 shows our room temperature isotherm along with data from Heinz and Jeanloz,heinz  and Bell et al.bell\_sccm  Also shown are some of the proposed EOS standards. Jamieson et al.jamieson  used the Hugoniot as a reference, and calculated the room temperature isotherm using a Mie-Grüneisen EOS. They assumed $`\gamma /V=\mathrm{const}`$, and used $`\gamma (V_0)=3.215`$, so their $`\gamma `$ is always larger than ours. This results in lower room temperature pressures than ours, but by restricting their analysis to $`V/V_0>0.775`$, the impact of the assumed $`\gamma `$ is minimized. The Heinz and Jeanloz EOS is based on extrapolating their room temperature data, with some consideration of the shock data. Our room temperature isotherm is in good agreement with Heinz and Jeanloz to 200 GPa. The recently proposed EOS of Shim et al.shim  is significantly lower than all the other standards. It is lower than the present analysis by 10 GPa at $`V/V_0=0.7`$, corresponding to $`P=156`$ GPa.

Figure 7 shows $`P(T)`$ along two isochores in comparison with various equations of state. At $`V/V_0=0.8`$, Jamieson et al. give somewhat lower $`P`$ at 300 K with a larger $`dP/dT`$, in keeping with their higher $`\gamma `$. Heinz and Jeanlozheinz  are in the best overall agreement with the present pressures, with somewhat lower $`dP/dT`$. Shim et al.shim  give a somewhat higher $`dP/dT`$ with $`P`$ generally low due to their large offset at room temperature. Anderson et al.anderson\_au  adopted the room temperature isotherm from Heinz and Jeanloz.heinz  Their EOS gives $`dP/dT`$ substantially too low at compressed volumes.

We have tabulated $`P(V,T)`$ in table 2 in the same format adopted by other authors. As discussed, we have extended the range of densities to $`V/V_0=0.6`$ and the temperatures to 5000 K, except where such temperatures are thought to lie in the liquid phase. Boettgerboettger\_au  has proposed an extension of the Au 300 K standard to 500 GPa based on LDA calculations. He gives $`P=344`$ GPa at $`V/V_0=0.6`$ compared to 329 GPa from our semi-empirical EOS. Given the scatter in the high pressure Hugoniot data, and our approximate treatment of melting, this 4.6% difference is probably within the uncertainties of the present analysis at this high density.

## V Melting and the Hugoniot

In order to investigate the effects of melting on the Hugoniot, we have constructed a two-phase model free energy. The ion free energy in the liquid is based on the assumptions that $`C_V^{\mathrm{ion}}=3Nk_B`$, which is reasonable for temperatures near melting, and that $`\mathrm{\Delta }S_V^{\mathrm{ion}}`$, the entropy difference between solid and liquid at fixed volume is $`0.8Nk_B`$. These are empirically based model assumptions.grover  A statistical mechanical basis for these observations is discussed by Wallace,wallace\_ldyn  who argues for the universality of $`\mathrm{\Delta }S_V^{\mathrm{ion}}`$. Beyond these assumptions, we require the energy of the liquid with respect to the solid to fully determine the liquid free energy. We do this by imposing that the melting curve, obtained by equating the pressures and Gibbs free energies of liquid and solid, follow the Lindemann rule in the form

$$\frac{T_m}{\omega _0^2(V_s)V_s^{2/3}}=\mathrm{const},$$
(19)
where $`V_s`$ is the volume on the solidus. The electron excitation free energy is assumed to be the same in the liquid as in the solid. Application of this method to Cu leads to a shock melting threshold of 228 GPa, which compares well with 232 GPa obtained by Hayes et al.hayes\_cu  by analyzing sound speed data. The Au Hugoniot has been calculated for the two-phase model allowing for coexistence in the shocked state.greeff\_ti  Figure 8 shows the resulting Hugoniot in the pressure-volume plane. Also shown in the figure are the data points and the linear $`U_s(U_p)`$ fit. The dot-dashed curve is the Hugoniot for the solid only. The boundaries of the coexistence region are visible as kinks on the solid curve at 280 and 350 GPa. Above 350 GPa, the Hugoniot is in pure liquid phase. A similar anomaly is shown for Al in a two-phase calculation by Chisolm et al.chisolm\_al  A significant enhancement of the Hugoniot pressure of Au due to melting was also shown in calculations by Godwal et al.godwal , however they predict a much lower shock melting threshold than we do. Just above complete melting, the liquid Hugoniot lies above the solid by 30 GPa. The linear $`U_s(U_p)`$ curve goes smoothly through the melting region to intersect the liquid data.

Our two-phase EOS gives good agreement with the high pressure data, which suggests that our cold energy is valid to compressions of $`V/V_00.6`$. Electronic excitations have practically no effect on the $`P(V)`$ Hugoniot in the solid phase, although they significantly affect the temperature. At higher temperatures in the liquid, they act to soften the Hugoniot by absorbing energy. Coincidentally, neglecting both electronic excitations and melting gives a Hugoniot that agrees with the linear $`U_s(U_p)`$ curve quite well to 650 GPa. This is an accidental cancellation of errors. Ignoring melting and electronic excitations gives a temperature that is too high by $`10^4`$ K at this pressure. This offsets the neglect of the pressure enhancement due to melting. It is not recommended to ignore either melting or electronic excitations in this high temperature regime.

## VI Conclusions

By combining ab initio calculations of elastic moduli and zone boundary phonons with interpolations based on a force constant model, we have calculated moments of the phonon frequencies of fcc Au to twofold compression. This allows us to calculate the associated Grüneisen parameters. In particular, we have focused on $`\gamma _0`$ which corresponds to the classical limit. We emphasize that the classical limit dominates the thermal pressure in the EOS. We find that the frequently used form $`\gamma (V)=\gamma (V_0)(V/V_0)^q`$ does not represent the volume dependence of $`\gamma `$ well. We have used an expansion to $`2^{\mathrm{nd}}`$ order in $`V`$ with a physical asymptotic limit, which fits the calculated moments well and is consistent with the measured value of $`\gamma `$ at ambient pressure and temperature. Using the theoretical $`\gamma (V)`$ and electron excitation free energy, we have constructed a semi-empirical free energy for the solid, which we believe to be as accurate as possible. This free energy is given as a parameterized closed form expression and the resulting $`P(V,T)`$ is given in tabular form to $`V/V_0=0.6`$.

Our static lattice energy is empirical, and is verified by comparison with the measured Hugoniot, so that our EOS can be regarded as a generalization from the Hugoniot with a physically based $`\gamma `$. In the solid phase, electronic excitations have a small effect. This, together with the dominance of the classical limit in the vibrational free energy, means that the widely used Mie-Grüneisen approximation, that $`\left(P/E\right)_V`$ is independent of $`T`$, is accurate. At the highest compressions the physics affecting the Hugoniot becomes more complicated. Melting is predicted to begin at 280 GPa and $`V/V_0=0.656`$. Complete melting is estimated to lead to a 30 GPa increase in the Hugoniot pressure over the solid. At these densities and higher, the Hugoniot temperature is rising rapidly and electronic excitations are playing an increasing role. The highest Hugoniot data are at 580 GPa, where the temperature is calculated to be above $`2\times 10^4`$ K. The Mie-Grüneisen approximation is not expected to be valid for interpolating between the Hugoniot and room temperature at this high compression, because of the strong effects of melting and electronic excitations. These effects need to be taken explicitly into account, as has been done here.

Regarding EOS standards, our analysis gives a room temperature $`P(V)`$ curve that agrees well with that of Heinz and Jeanlozheinz  to 200 GPa. Our EOS gives $`\left(P/T\right)_V`$ generally somewhat larger than theirs, and we believe that in this regard our EOS is preferred. The EOS of Jamieson et al.jamieson  is based on reduction of shock data, and was originally given to 80 GPa. It gives somewhat lower pressures than ours, which is due to their use of a Grüneisen parameter that is too large at all volumes. Extrapolation of their room temperature isotherm to higher pressuresakahama  is not recommended. Anderson et al.anderson\_au  adopted the room temperature isotherm of Heinz and Jeanlozheinz , while giving a different thermal pressure. Their EOS gives $`\left(P/T\right)_V`$ which is substantially too low under compression, and is not recommended for high temperatures. Simultaneous compression of Au and Pt showedakahama  that the Pt standard of Holmes et al.holmes\_pt  gave a pressure higher by 15 and 20 GPa than the Au standards of Heinz and Jeanlozheinz  and Jamieson,jamieson  respectively, at a pressure of 150 GPa. While the extrapolated Jamieson isotherm is expected to be somewhat low, our analysis agrees with the room temperature isotherm of Heinz and Jeanloz, suggesting that the discrepancy between the Au and Pt pressures is due to errors in the Pt standard. Given the importance of accurate pressure standards, the status of the Pt EOS seems to warrant further investigation.

###### Acknowledgements.

We thank J. C. Boettger, J. D. Johnson, E. D. Chisolm and S. Crockett for helpful discussions. This work was supported by the U. S. Department of Energy under Contract No. W-7405-ENG-36.


format: text

text:
# 𝐵_gg revisited: The environments of low-excitation radio galaxies and unified models

## 1 Introduction

In unified models for powerful radio galaxies and radio-loud quasars (Scheuer 1987, Barthel 1987, 1989) isotropic properties of the two classes of objects, such as their clustering environments, should be statistically identical. However, the classical unified models work best at high redshift (there are few low-redshift objects classified as radio-loud quasars) and at these redshifts the clustering environments of the sources are hard to measure. Some attempts have been made to test the predictions of unified models using X-ray observations of extended emission, assumed to be cluster-related (Hardcastle & Worrall 1999) but sample sizes are generally small and ROSAT-based work is generally not reliable because of the contamination of the X-ray emission by AGN- and radio-source-related components. Chandra should provide significantly better constraints from the X-ray once sufficient objects are in the archive, though it will be important to take into account the significant evolution undergone by the ICM of a radio-source host cluster between $`z1`$ and 0. In the meantime, there is still an important rôle for the traditional method of calculating the amplitude of the galaxy-galaxy or quasar-galaxy spatial covariance function, hereafter referred to as $`B_{\mathrm{gg}}`$ (e.g. Longair & Seldner 1979; Yee & Green 1987; Prestage & Peacock 1988; Yates, Miller & Peacock 1989; Ellingson, Yee & Green 1991; Hill & Lilly 1991; Yee & Ellingson 1993). Although it is recognised that selection in more than one optical filter is desirable to achive good statistics at high redshifts (Barr et al. 2003), at lower redshifts simple $`B_{\mathrm{gg}}`$ determinations seem likely to be able to provide a statistical measure of the environments of different classes of source.

Recently Harvanek et al. (2001: hereafter HESR) have carried out a large systematic study of the clustering properties of 3CR (Spinrad et al. 1985) radio sources. By making new observations, and collating and cross-calibrating data from the literature, they were able to obtain a nearly complete set of $`B_{\mathrm{gg}}`$ values in the redshift range $`0.15<z<0.65`$. 3CR objects in this redshift range are almost all luminous FRII-type radio galaxies or lobe- or core-dominated quasars. Yee & Green (1987) had shown that few low-redshift, radio-loud quasars are located in rich environments. HESR aimed to test the unified-model prediction that radio galaxies at these redshifts would also tend to lie in poor environments. They obtained the surprising result that the $`B_{\mathrm{gg}}`$ values of their 3CR radio galaxies were not consistent with those of a sample of radio-loud quasars compiled by Yee & Ellingson (1993: hereafter YE93) in the redshift range $`0.15<z<0.4`$, in the sense that a number of the radio galaxies were found in significantly richer environments than any quasar. An immediate objection to this result is that the YE93 quasars were not in general as radio-luminous in extended emission as the 3CR comparison sample (otherwise they would all have been 3CR objects themselves on the basis of their low-frequency extended emission, whereas only a fraction of them are in fact in 3CR). The effects of this bias are hard to quantify, but since (for equal jet powers and source ages) simple radio-galaxy physics suggests that a more luminous object will lie in a richer environment, it might tend to produce a difference in the sense of the one observed by HESR. In what follows I shall ignore this potential problem and concentrate instead on the effects on HESR’s conclusions of using a more sophisticated variant of the unified model.

## 2 Unification at low redshifts and the low-excitation objects

At low redshifts (or, more generally, low luminosities) the simple unified model of Barthel (1989) breaks down, most obviously because of the absence of low-redshift FRII quasars. I have argued elsewhere (Hardcastle et al. 1998) that the broad-line radio galaxies (BLRG) take the place of quasars at low redshifts/luminosities, although that does not rule out the possibility that more luminous BLRG are objects intermediate in viewing angle between quasars and narrow-line radio galaxies (NLRG) (e.g. Dennett-Thorpe et al. 2000). However, low-excitation radio galaxies (LERG) (Hine & Longair 1979; Laing et al. 1994), which are a significant population at low FRII luminosities, cannot be unified with either BLRG or quasars if the narrow-line emission is isotropic. Instead, it has been suggested (Barthel 1994, Laing et al. 1994), that they form a separate population, possibly unified at small angles to the line of sight with FRII-like BL Lac objects. They have other anomalous properties: Hardcastle et al. (1998) noted that the LERG in their $`z<0.3`$ 3CR/3CRR sample were on average significantly smaller than the high-excitation objects, and that they included a number of objects with prominent jets and relaxed or absent hotspot structure. Hardcastle & Worrall (1999) found that the LERG in their ROSAT-observed subsample of 3CRR tended to have relatively luminous X-ray environments, whereas at least some of the NLRG at similar redshifts were in poorer environments. This motivates a re-analysis of the results of HESR taking into account the emission-line type of the radio galaxies.

## 3 Data and analysis

The data for 3CR objects, including the emission-line classifications, are taken from HESR, with emission-line classifications supplemented by data from the literature (Table 1). HESR took most of their emission-line classifications from Jackson & Rawlings (1997). Lacking new spectroscopic observations, I have classified objects that were not classified in Jackson & Rawlings mostly on the basis of qualitative statements in the literature, such as whether they show ‘strong’ or ‘weak’ emission lines and what species are present in the spectrum. These classifications must be regarded as best guesses only, but, as discussed by Hardcastle et al. (1998), the approach tends to give the same results as a more quantitative determination. To check this quantitatively, I examined all galaxies that are classified both in Jackson & Rawlings and in Spinrad et al. (1985), comparing Spinrad et al.’s qualitative spectral classification with the Jackson & Rawlings class. The two papers were taken to agree if an object classed as Spinrad et al. as ‘WE’ or ‘ABS’ was classed by Jackson & Rawlings as a low-excitation object, or if an object classed by Spinrad et al. as ‘SE’ was classed by Jackson & Rawlings as a high-excitation object or ‘weak quasar’ (effectively a broad-line radio galaxy). Of the 97 overlapping objects, the two papers disagreed on only 8. This gives us some confidence that a qualitative approach will produce acceptable results. One object, 3C 142.1, remains unclassified, but lies outside the redshift range of particular interest in any case. Following HESR, 3C 275 and 3C 435A were omitted from the analysis. The data for HESR’s comparison sample of radio-loud quasars were taken directly from YE93; HESR’s $`B_{\mathrm{gg}}`$ values were taken in preference to YE93’s where the sample had objects in common.

Fig. 1 shows the distribution of $`B_{\mathrm{gg}}`$ values for the objects in the 3CR sample with $`0.15<z<0.4`$, classified by emission-line properties. This may be compared with the left-hand panel of HESR’s figure 6. It is clear that the low-excitation objects have a different distribution of $`B_{\mathrm{gg}}`$ from the high-excitation radio galaxies and quasars; in fact, all the 3CR objects in the richest environments found by HESR are LERG.

HESR used a Kolmogorov-Smirnov test to investigate the differences in $`B_{\mathrm{gg}}`$ distribution between the radio galaxies and quasars in their sample with $`0.15<z<0.4`$. To compare with HESR’s approach, I use a quasar sample consisting of 21 YE93 quasars (the 24 objects in the redshift range, omitting the 3 quasars in the HESR 3CR sample) plus the 5 quasars from the 3CR sample, giving 26 quasars in all. In the HESR redshift range there are 31 radio galaxies (including BLRG and N-galaxies; here I class these with the NLRG for consistency with HESR). The K-S test finds only a 5% probability that these two samples are drawn from the same parent distribution, consistent with HESR’s value of 3%. But if the 8 LERG are excluded from the radio galaxy sample, the K-S probability that the two are from the same parent population rises to 50%; there is no significant difference between the distributions of radio galaxies and quasars in these samples in this redshift range.

We can restrict ourselves to the 3CR-derived sample with $`0.15<z<0.4`$ in order to test unified models without including the lower-radio-power YE93 quasars. Here, in order to obtain a sufficiently large sample of sources that make a small angle to the line of sight to apply statistical tests, I include the BLRG with the quasars; this is a valid procedure whether or not some of the BLRG are intermediate-angle objects. In this sample, the $`B_{\mathrm{gg}}`$ distribution of the 12 broad-line objects (quasars and BLRG) has only a 5% probability of being drawn from the same population as the 24 LERG and NLRG combined (cf. Fig. 1), but a 72% probability of being drawn from the same population as the 16 NLRG; thus, if we neglect the difference between LERG and NLRG, there is a problem for low-redshift unification as described by Hardcastle et al. (1998), while if we remove the LERG there is no problem. On the other hand, the $`B_{\mathrm{gg}}`$ distribution of the LERG with $`0.15<z<0.4`$ has a 0.04% probability of being drawn from the same population as the other 3CR objects, a highly significant result. To summarize, all the differences between the distributions of radio galaxies and quasars/broad-line objects in the redshift range $`0.15<z<0.4`$ become statistically insignificant if the LERG are excluded from the radio galaxy sample; the LERG have significantly different properties from the other objects.

## 4 Discussion and conclusions

The key result of this work is that 3CR LERG lie in significantly richer environments than the HERG/quasar population at comparable redshifts. This strong difference between the environmental properties of LERG and HERG supports suggestions by Hardcastle et al. (1998) and Hardcastle & Worrall (1999) that the 3CR LERG population is physically different from the HERG. If we believe that the emission-line properties of the radio source reflect the power of the jet (e.g. Rawlings & Saunders 1991) then a plausible explanation is that at least some LERG are intrinsically low-jet-power radio galaxies which owe their comparatively high overall radio luminosity and their other properties (small size, dissipative jets, relaxed appearance – cf. Hardcastle et al. 1998, Harvanek & Stocke 2002) to a rich environment. This recalls, and may be related to, the arguments of Barthel & Arnaud (1996).

The fact that the narrow-line objects and quasars have similar $`B_{\mathrm{gg}}`$ distributions, when the LERG are excluded, supports the version of the orientation-based unified model for radio galaxies and quasars outlined in section 2.

###### Acknowledgements.

I thank the Royal Society for a Research Fellowship, and Michael Harvanek for discussions of $`B_{\mathrm{gg}}`$ values for radio galaxies in the mid-1990s which helped to alert me to the interesting properties of LERGs. I am grateful to Michael Harvanek, John Stocke and Judith Croston for their helpful comments on an earlier draft of this paper, and to Philip Best for a rapid and constructive referee’s report.


text:
# Cooling by Heat Conduction Inside Magnetic Flux Loops and the Moderate Cluster Cooling Flow Model

## 1 Introduction

Recent Chandra and XMM-Newton observations of cooling flow (CF) clusters of galaxies have lead to a shower of papers examining reheating processes of the intracluster medium (ICM), which together with pre-Chandra and XMM-Newton reheating papers make a long list of papers with different ideas (e.g., Binney & Tabor 1995; Tucker & David 1997; Ciotti & Ostriker 2001; David et al. 2001; Quilis, Bower, & Balogh 2001; Brüggen, M. & Kaiser 2001; Ruszkowski & Begelman 2002; Nulsen et al. 2002; Fabian 2003; Fabian et al. 2003; Binney 2003; Mathews et al. 2003; for more references see Peterson et al. 2003b). These recent observations, although in contradiction (Fabian 2003) with the cluster-CF models of 10 years ago (Fabian 1994), are compatible with models having low mass cooling rates (see Binney 2003, for why the new results are not a surprise), such as in the moderate CF model (Soker et al. 2001). The moderate CF model, which was proposed before the new results from the two X-ray telescopes, is different from many early proposed processes which aim is to prevent the CF in cluster of galaxies altogether. The main ingredient of the moderate CF model is that the effective age, i.e., the time period since the last major disturbance of the ICM inside the cooling radius $`r_c100\mathrm{kpc}`$, is much shorter than the cluster age (e.g., Binney & Tabor 1995; see Binney 2003; Soker et al. 2001). The cooling radius is defined as the place where the radiative cooling time equals the cluster age. Originally the heating in the moderate CF model was proposed to be intermittent (Soker et al. 2001), but the basic idea can hold for a steady heating, or heating in short intervals (Binney 2003).

Kahn et al. (2003) and Peterson et al. (2003a), among others, find that there is less gas than that predicted by old versions of the steady-state isobaric CF model at all temperatures below the ambient gas temperature. Also, the discrepancy increases with decreasing temperature. There is not just a deficit in gas below $`1/3`$ of the ambient cluster temperature, there is a deficit of gas at all temperatures below the ambient temperature (Peterson et al. 2003b). Soker & David (2003) show that these XMM-Newton results are consistent with the expectations of the moderate CF model. Soker & David (2003) claim that the X-ray emitting gas in CF clusters can only be in steady-state if there exists a steady heating mechanism that scales as $`H(T)T^\alpha `$ where $`\alpha =12`$; this is unlikely. They show then that in the moderate CF model most of the gas within CF resides in the hottest gas, which is prevented from cooling continuously and attaining a steady-state configuration by being reheated. This results in a mass cooling rate that decreases with decreasing temperature, with a much lower mass cooling rate at the lowest temperatures, as discovered by Peterson et al. (2003a). Although there is no indication yet for a violent heating, I note the following. An intermittent reheating can work for an AGN kinetic energy of $`10^{47}\mathrm{erg}\mathrm{s}^1`$, and its strong activity should last $`12\times 10^7`$ yr and occur every $`10^9`$ yr (Soker et al. 2001). A violent heating event by AGN jets was proposed by Kawano, Ohto, & Fukazawa (2003) to have occurred $`>10^9`$ yr ago in the CF clusters 2A 0335 + 096 and A 2199. Mazzotta, Edge, & Markevitch (2003) claim that the structure of the cluster 2A 0335+096 can be accounted for if strong AGN activity periods are followed by relative quiescence periods. O’Sullivan & Vrtilek (2003) propose that a heating event occured $`2\times 10^8\mathrm{yr}`$ ago in the poor cluster MKW 4, which consequently is presently not in a steady state cooling. Such a scenario is in accord with the explanation of Soker & David (2003) to the observations of Peterson et al. (2003a) and Kahn et al. (2003). Only $`12\%`$ of all CF clusters should be found during that stage. In addition, Heinz & Enßlin (2003) argue for bias against detecting strong shocks. Therefore, from the well studied CF clusters, only one, or at most two, are expected to be in the stage of violent heating. In Cygnus A there is a strong radio source which heats the ICM (Smith et al. 2002). Wilson, Young, & Smith (2003) estimate the mechanical power of the jets in Cygnus A to be $`L_{\mathrm{jet}}6\times 10^{46}\mathrm{erg}\mathrm{s}^1`$. This is $`100`$ times larger than the radio emission, and it is in the range required by the moderate CF model. Therefore, present observations do not rule-out a violent intermittent reheating process acting for a short time every $`10^9`$ yr.

The limit on the cooling rate below a temperature $`T_{\mathrm{min}}`$ inferred from X-ray observations is $`<20\%`$ of the mass cooling rates cited in the past (Fabian 2002). That is, X-ray emission from gas below $`T_{\mathrm{min}}`$ is absent in present X-ray telescopes. This minimum temperature is different for different clusters, but it is $`1/3`$ of the background temperature (Kaastra et al. 2003; Peterson et al. 2003b). Despite the missing X-ray emission reported by these studies, there are other indications that cooling to temperatures much below $`T_{\mathrm{min}}`$ does occur, but on a much lower rate than that inferred by taking the CF age to be equal to the cluster age. Wise, McNamara, & Murray (2003) and McNamara, Wise, & Murray (2003), for example, find the cooling rate within $`r30\mathrm{kpc}`$ of the CF cluster A 1068 to be about equal to the star formation rate there $`2070M_{}\mathrm{yr}^1`$. Although the limits on the mass cooling rates cited above are within the values expected in the moderate CF model, a much better agreement with observations can be achieved if some non-radiative cooling of the X-ray gas is added to the moderate CF model (Soker, Blanton, & Sarazin 2003, hereafter SBS). The thermal energy released in cooling the hot gas is delivered to the cold clouds, and is radiated as optical/UV line emission at $`10^410^5`$ K. This energy of the hot gas can be delivered to the cool gas by mixing (e.g., Oegerle et al. 2001; Fabian et al. 2001, 2002; Johnstone et al. 2002; Bayer-Kim et al. 200), or via heat conduction (SBS)111This local heat conduction model is completely different from large scale ($`100`$ kpc) heat conduction models aiming at heating, rather than cooling, the X-ray emitting gas in the center of CF clusters.

In the scenario proposed by SBS the heat conduction is regulated by reconnection between the magnetic field lines in cold ($`10^4`$ K) clouds and the field lines in the ICM. A narrow conduction front is formed, which, despite the relatively low temperature, allows efficient heat conduction from the hot ICM to the cold clouds. The reconnection between the field lines in cold clouds and those in the ICM occurs only when the magnetic field in the ICM is strong enough. This occurs only in the very inner regions of CF clusters, at $`r1030`$ kpc. While in the scenario of SBS the non-radiative cooling starts with the reconnection of magnetic field lines from the pre-existing cold clouds to the X-ray emitting gas, in the present paper I examine heat conduction within pre-existing magnetic flux loops of X-ray emitting gas. The non-radiative cooling starts when small portion of the loop is radiatively cooling to $`10^4\mathrm{K}`$.

Magnetic flux loops in CF clusters were studied before, e.g., in relation to global heat conduction (Norman & Meiksin 1996), in the amplification of magnetic fields in the center of CF clusters (Godon, Soker, & White 1998), and the formation of multiphase ICM (Zoabi, Soker, & Regev 1996, 1998; hereafter ZSR). In a recent paper Kaastra et al. (2003) note the similarity of the differential emission measure in the clusters they study to that of coronal loops in the sun. The behavior of loops found by ZSR, the finding of Kaastra et al. (2003), and the dynamically important magnetic field in the inner regions of cluster CFs (Eilek & Owen 2002; Taylor, Fabian, & Allen 2002), motivate me to study non-radiative cooling within magnetic flux loops. This paper is not aiming at completely suppressing CF, as there are many evidence for cooling to low temperatures, e.g., Kaastra et al. (2003) find radiation from gas at temperatures of $`10^6\mathrm{K}`$ in some clusters, Edge et al. (2002; see also Edge & Frayer 2003; Salome & Combes 2003a,b; Edwards 2003) find massive cold reservoir in the form of molecular gas, probably from cooling ICM, and the star formation rate cited above (Wise et al. 2003). The main goal is to study another aspect, i.e., non radiative cooling of X-ray emitting gas in magnetic flux loops, and incorporate it into the moderate cluster CF model. As such, this paper is complementary to SBS.

## 2 The Rapid Radiative Cooling Phase of the Loop

ZSR consider an elliptically shaped magnetic flux loop moving radially inward in a cluster CF. They follow only the two end points of the loop, assuming that they evolve in series of hydrostatic equilibria, i.e., magnetic tension $`T_B`$ balances buoyancy $`F`$ at each of the two end points of the loop. This dictates that the density in the upper segment of the loop $`\rho _u`$, is lower than that of its surroundings, and the density of the lower segment of the loop $`\rho _d`$, is higher than that of its surroundings. ZSR consider radiative cooling, but neglect heat conduction, even along the magnetic field lines. They base this assumption on a possible reduction in the heat conduction coefficient along magnetic fields lines (Pistinner & Eichler 1998). This reduction in heat conduction is not efficient when the magnetic pressure becomes comparable to the thermal pressure (Pistinner & Eichler 1998). ZSR take the magnetic pressure to increase as density increases, and take a pressure equilibrium between the loop and its surroundings. As the loop evolves and falls toward the cluster center it cools with the lower-segment temperature $`T_d`$ falling faster than the upper-segment temperature $`T_u`$. The loop is stretched in the radial direction and gets narrow on its sides (see fig. 1). ZSR find that because the magnetic tension prevents the dense lower segment of the loop from falling toward the cluster center, this lower segment reaches a phase of catastrophic cooling. For the parameters used by ZSR the catastrophic cooling occurs when the lower segment is inside a cluster radius of $`r<20\mathrm{kpc}`$, but the upper regions of the loop extend to much larger distances, depending on its size.

In the present paper I examine possible effects of heat conduction along magnetic field lines inside the loop. The two relevant time-scales are the conduction time scale along the loop, from the upper hot segment to the cooling lower segment, and the radiative cooling time at the denser lower segment. Heat conduction proceeds on a time scale of $`\tau _{\mathrm{con}}(5n\mathrm{k}T/2)(\chi T^{7/2}L^2)^1`$, where the first parenthesis is the enthalpy per unit volume, and the second parenthesis is the heat loss per unit volume per unit time because of heat conduction. The temperature $`T`$ is some average inside the loop, which at early stages is close to both $`T_u`$ and $`T_d`$. Before the lower segment of the flux loop rapidly cools to $`T_d10^7\mathrm{K}`$, the temperature gradient inside the loop is shallower than $`T_u/L`$, and the heat conduction time is longer than the one given above. Here, $`\chi `$ is (nearly) a constant defined such that $`\chi T^{5/2}`$ is the heat conduction coefficient, $`n`$ is the total number density, and $`L`$ the width of the conduction front, taken here to be the length of the magnetic flux loop (see figure 1).

Scaling with typical values gives

$$\tau _{\mathrm{con}}4\times 10^9\left(\frac{n_e}{0.04\mathrm{cm}^3}\right)\left(\frac{T}{10^7\mathrm{K}}\right)^{5/2}\left(\frac{L}{10\mathrm{kpc}}\right)^2\mathrm{yr},$$
(1)
where $`n_e`$ is the electron density. At early stages, the time $`\tau _{\mathrm{con}}`$ is the time required to transfer heat from the upper segment to the lower one to prevent it from cooling to very low temperatures.

The isobaric radiative cooling time of the lower segment in the temperature range $`2\times 10^5T4\times 10^7`$ K is (see cooling curve in, e.g., Gaetz, Edgar, & Chevalier 1988)

$$\tau _{\mathrm{rad}}2\times 10^8\left(\frac{n_e}{0.04\mathrm{cm}^3}\right)^1\left(\frac{T_d}{10^7\mathrm{K}}\right)^{3/2}\mathrm{yr}.$$
(2)
The heat conduction will not prevent the lower segment from cooling to lower temperatures when the condition $`\tau _{\mathrm{con}}>\tau _{\mathrm{rad}}`$ is met. Using equations (1) and (2), and substituting $`T_d`$ for $`T`$ (this is reasonable before the catastrophic cooling of the lower segment starts), the condition for rapid cooling of the lower segment reads

$$T_d1.7\times 10^7\left(\frac{Tn_e}{4\times 10^5\mathrm{K}\mathrm{cm}^3}\right)^{1/3}\left(\frac{L}{10\mathrm{kpc}}\right)^{1/3}\mathrm{K},$$
(3)
where I took a constant pressure (expressed as $`Tn_e`$, and scaled with $`P=10^{10}\mathrm{erg}\mathrm{cm}^3`$), such that the condition is on the temperature of the lower loop’s segment. It is evidence that the temperature when rapid cooling starts depends weakly on the external pressure and on the loop size. For the typical parameters used here, it is in the range of the minimum temperature of X-ray emitting gas deduced from observations (Peterson et al. 2003a; Kaastra et al. 2003).

To summarize this section, there are two conditions for the catastrophic cooling of the lower segment of a flux loop to occur. First, the magnetic field must be strong, close to equipartition, in order to prevent the dense segment from falling in too fast. This condition is also true for the scenario of SBS, where the strong magnetic field is required to initiate reconnection between cold clouds and the ICM. Second, the lower segment temperature must reach the value given by equation (3). Both conditions are likely to be met only near the cluster’s center.

## 3 The Quasi-Static Heat Conduction Front Phase

Borkowski, Balbus, & Fristrom (1990) numerically study the evolution of a conduction front between a gas at $`T_u10^6\mathrm{K}`$ and cooler gas at $`10^4\mathrm{K}`$. They start with a sharp boundary between the two media. After a short period of time, during which the conduction front widens, a quasi-static front is developed. In the third phase$``$the condensation phase$``$the entire front radiatively cools down. During the quasi-static front phase Borkowski et al. (1990) find the mean temperature of different ions to be $`T_{\mathrm{ions}}14\times 10^5\mathrm{K}`$, i.e., much below the temperature of their hot phase which is $`T_u510\times 10^5K`$. I consider a similar effect here, where the mean temperature is low, with most radiation emitted at much lower temperatures than $`T_u`$. The initial situation here is not a sharp boundary between the hot and cold gases, but rather the lower segment of the loop cools from a high temperature. Therefore, I expect the quasi-static front phase to be established quite early.

An accurate solution for the temperature distribution along a static magnetic flux tube was derived analytically by Rosner, Tucker, & Vaiana (1978) and numerically by Aschwanden & Schrijver (2002). Both works consider heating as well as heat conduction and radiative cooling. Kaastra et al. (2003) when comparing their spectra to the results of Aschwanden & Schrijver (2002) assume the existence of a heating source, but they don’t specify its nature. I here neglect heating, because I assume sporadic heating events in the cluster, which most loops do not encounter during their catastrophic cooling. The results of Rosner et al. (1978) can be reasonably fitted along most of the tube by $`Ts^n`$ with $`n=1/3`$, where $`s`$ is the coordinate along the flux tube, i.e., along the magnetic field lines, and $`s=0`$ is the location of the already cold gas; here it is the lower segment of the loop. Only at the highest temperature segment the temperature pofile becomes shallower. Practically the temperature of the lower segment after its catastrophic cooling can be set to $`T=0`$, although it is $`10^4\mathrm{K}`$. The numerical results of Aschwanden & Schrijver (2002) depend on the assumptions used, but are crudely fit by $`n0.3`$. These results fit the quasi-static heat conduction phase discussed here, as I now show.

A simple balance between radiative cooling and heat conduction at each point along the flux loop reads (in the temperature range $`2\times 10^5T4\times 10^7`$)

$$\frac{d}{ds}\left(\chi T^{5/2}\frac{dT}{ds}\right)=n_e^2\mathrm{\Lambda }_e,$$
(4)
where $`\mathrm{\Lambda }_e`$ is the cooling function, with $`\mathrm{\Lambda }T^{1/2}`$ in the relevant temperature range. To obtain a simple analytical solution I take a constant pressure along the flux tube, such that $`n_eT^1`$. This allows the last equation to be written as

$$\frac{d^2(T^{7/2})}{ds^2}=A_1T^{5/2},$$
(5)
where$`A_1`$ is a constant composed of the different coefficients in equation (4). This can be integrated once to yield

$$T^{5/2}\frac{dT}{ds}=2\left(\frac{A_1}{7}\right)^{1/2}(T+C_1)^{1/2},$$
(6)
where $`C_1`$ is a constant of integration to be determined by boundary conditions. The left hand side of the last equation is the heat flux. Rosner et al. (1978) take this flux to be equal to zero at the higher temperature end of the loop. However, we see that such a boundary condition here requires $`C_1=T_u`$, which is not possible for lower temperatures. A constant heating along the flux tube can make this solution possible, by balancing the slow cooling of the high temperature (and low density) segment. Otherwise heat conduction balances radiative cooling and the heat flux can’t vanish, and on the hottest end of the loop no real steady state situation exists.

Equation (4) can analytically be integrated in principle. However, the expression gives $`s`$ as function of $`T`$ and it is too complicated for practical uses. Instead I look for a solution in the form $`Ts^n`$. I start by taking a minimum heat flux, which vanishes when $`T=0`$. This gives $`C_1=0`$. The solution reads (with $`s=0`$ where $`T=0`$)

$$T=A_2s^{1/3},$$
(7)
where $`A_2=6^{1/3}(A_1/7)^{1/6}`$. As this is close to the solution obtained by Rosner et al. (1978), I consider it a satisfactory approximation. $`C_1>0`$ gives a lower value of $`n`$. For the extreme case of $`C_1T`$ in the entire range (which implies a huge heat flux, hence it is not appropriate here where no heating is considered), the solution gives $`n=2/7`$. The variation of $`T`$ with $`s`$ is somewhere between these two extremes, i.e., $`(2/7)<n<(1/3)`$, but very likely close to $`Ts^{1/3}`$. For a constant pressure along the loop, the relative emission per unit volume is given by $`\dot{e}=n_e^2\mathrm{\Lambda }_eT^{5/2}`$. For a temperature given by $`Ts^n`$, and a constant pressure along the loop, the luminosity fraction of gas at temperature $`<T`$ is given by (assuming constant loop’s cross section)

$$\eta (T)=\left(_0^TT^{5/2}𝑑s\right)\left(_0^{T_u}T^{5/2}𝑑s\right)^1=\left(\frac{T}{T_u}\right)^{\frac{1}{n}\frac{5}{2}},$$
(8)
where in the second equality I substituted $`ds=T^{(\frac{1}{n}1)}dT/n`$, and preformed the integration. Note that $`T_u`$ is not the maximum temperature in the CF cluster $`T_{\mathrm{max}}`$, but rather the maximum temperature inside the specific magnetic flux loop under consideration, and $`T_u<T_{\mathrm{max}}`$.

For $`n=1/3`$ the luminosity of gas at temperature $`<T`$ goes as $`\eta T^{1/2}`$. Lower values of $`n`$ result in more emission from high temperature gas, but the higher value of $`n1/3`$ seems appropriate for the following reason. I assumed above that the pressure along the flux tube is constant. However, in clusters the upper segment of the loop will be further out in a lower pressure environment. Outside $`r25\mathrm{kpc}`$ in A 2052 (Blanton et al. 2001) and $`r10\mathrm{kpc}`$ in Hydra A (David et al. 2001), the pressure drops as $`Pr^1`$. This implies a pressure drop by a factor of $`1.5`$ between the two ends of a loop with $`L10\mathrm{kpc}`$. The lower pressure, hence density, at the upper loop segments implies less emission from the hot gas compared to the assumed constant pressure case. Therefore, taking the higher value of $`n`$ in the range given above is justified. For the $`n=1/3`$ case, hence $`\eta =(T/T_u)^{1/2}`$, $`30\%`$ of the luminosity comes from gas at $`T<0.1T_u`$, and $`50\%`$ from, gas at $`T<0.25T_u`$. For example, if $`T_u=1\mathrm{k}\mathrm{e}\mathrm{v}`$ (recall that $`T_u`$ is not the maximum temperature in the CF cluster), then $`30\%`$ ($`50\%`$) of the luminosity comes from gas at temperatures $`T<0.1\mathrm{kev}`$ ($`T<0.25\mathrm{kev}`$). Overall, heat conduction within a magnetic flux loop can reduce the fraction of energy radiated in the X-ray band by a factor of $`1.52`$.

## 4 Discussion and Summary

I found above that heat conduction within a magnetic flux loop can reduce the fraction of energy radiated in the X-ray band by a factor of $`1.52`$, i.e., the X-ray luminosity is $`5070\%`$ of the value expected without heat conduction. Cooling of X-ray emitting gas by conduction to cold gas was proposed to occur via mixing (e.g., Oegerle et al. 2001; Fabian et al. 2001, 2002; Johnstone et al. 2002; Bayer-Kim et al. 200) and by local heat conduction regulated by magnetic field reconnection (SBS). In the later two processes the temperature gradient is much higher than the one inside the loop studied here, and hence these mechanism are more efficient in reducing X-ray emission. However, the question is whether these processes occur, and how often. The magnetic flux loop model studied here is motivated by previous suggestions that many such loops exist in the ICM of CF clusters and that they play an active role in cooling and/or dynamics of the gas (e.g., Zoabi et al. 1996; ZSB; Kaastra et al. 2003).

SBS estimate that in their heat conduction scenario via magnetic field reconnection the total line emission expected for gas cooling from $`10^7`$ K to $`10^6`$ K might be factor of $`30100`$ lower than that expected in older versions of the CF model, where only line emission is considered for the cooling process, and where the cooling rates are much higher. This factor comes from a $`10`$ times lower mass cooling rate of the hot gas in the moderate CF model (Soker et al. 2001), and from another factor of $`310`$ for cooling via heat conduction. This much lower expected radiative cooling rate at low X-ray emitting temperatures is compatible with limits inferred from observations. The observational limit, which is $`1020\%`$ of the expected cooling rate in old versions of the CF model (Fabian 2002; Peterson et al. 2003a), can be met if the reduction due to heat conduction is only a factor of $`2`$ as found here, with another factor of $`10`$ from the moderate CF model. Possibly, the more significant role of magnetic flux loops is to stir the ICM in the center of CF clusters, by reconnection and motion caused by magnetic tension (Zoabi et al. 1996; ZSR). This stirring, where some portions of loops falls fast and others buoyant outward, can lead to more reconnection between magnetic flux loops (Zoabi et al. 1996), enhancing the efficiency of non-radiative cooling via local heat conduction as discussed by SBS.

If the general scenario discussed here and/or in SBS holds, then a CF model with the following ingredients emerges.
(1) CF does occur. The evidences for CF in many clusters were summarized by Soker et al. (2001) prior to the new results of the Chandra and XMM-Newton X-ray telescopes. The new results show that the cooling rate to low temperatures is much lower than the expected values in the old versions of the CF model, but they do show that some cooling occurs (e.g., a few clusters in Kaastra et al. 2003). This is compatible with the moderate CF model. On the other hand, in most CF clusters the new X-ray observations didn’t reveal any violent heating events, e.g., major mergers or strong shocks, and the ICM structure is not consistent with global heat conduction from the outer hot massive clusters regions. The X-ray observations do show heating, e.g., bubbles and sound waves (e.g., Fabian et al. 2003), but these can’t compensate for the cooling in many cases (e.g., Kawano et al. 2003). Observations using the high spatial resolution of Chandra show most activity to occur in the very inner region, on a scale much smaller than the cooling radius (where radiative cooling equals the cluster age), compatible with much lower mass cooling rate. All these suggest that the new X-ray observations actually strengthen the CF model, although a moderate one, rather than “kill” it.
(2) Non-steady state evolution. Soker & David (2003) show how a non-steady state situation can account for the new observations. Such a non-steady state evolution is expected in intermittent-heating models, with long enough time intervals. The gas is heated up, and then cools. Most of the heated gas does not cool to low temperatures before the next major heating event. Soker & David claim that no fine tuning is required. More accurate analysis may show that this non-steady state evolution can hold for heating events with shorter intervals as well.
(3) Non-radiative cooling of X-ray emitting gas. Although not strictly a necessary ingredient, non-radiative cooling of the X-ray emitting gas after it cooled already by a factor of $`25`$ from its initial temperature, makes the agreement of the moderate CF with observations much better (SBS; this paper). The initial radiative cooling is required, most likely, to enhance by a factor of $`10`$ the magnetic pressure relative to the thermal pressure (SBS). The magnetic fields play a crucial role in the scenario discussed here, and in the more efficient non-radiative cooling scenario studied by SBS.
(4) A corona-like ICM. The outflow (buoyant bubbles and magnetic flux loops) and inflow, the strong magnetic fields, the possible importance of local heat conduction, radiative cooling in optically thin plasma, the presence of X-ray emitting gas together with optically emitting gas, all these and more, make the ICM in CF clusters similar in some aspects to the gas above the solar (and similar stars) photosphere, i.e., photosphere to corona regions. The similarity with respect to some processes relevant to magnetic field reconnection and optical line emission was mentioned by Jafelice & Friaca (1996) and Godon et al. (1998). Sparks & Donahue (2003) note the morphological similarity between filaments in M 87 and the solar atmosphere. Godon et al. (1998) discuss magnetic field amplification in CF clusters and the possibility of “cluster flares”, and use hints from solar flares to project on the expected nature of violent magnetic reconnection in CF clusters. Some aspects of magnetic flux loop dynamics where studied by Zoabi et al. (1996) and ZSR. Kaastra et al. (2003) extend the comparison of ICM in CF cluster with stellar coronae to explain the X-ray spectra they observed in CF clusters.

This research was partially supported by a grant from the Israel Science Foundation.


text:
# Shock in a Branching-Coalescing Model with Reflecting Boundaries

## Abstract

A one-dimensional branching-coalescing model is considered on a chain of length $`L`$ with reflecting boundaries. We study the phase transitions of this model in a canonical ensemble by using the Yang-Lee description of the non-equilibrium phase transitions. Numerical study of the canonical partition function zeros reveals two second-order phase transitions in the system. Both transition points are determined by the density of the particles on the chain. In some regions the density profile of the particles has a shock structure.

journal: Physics Letters A

One-dimensional driven lattice gases are models of particles which diffuse, merge and separate with certain probabilities on a lattice with open, periodic or reflecting boundaries. In the open boundaries case the particles are allowed to enter or leave the system from both ends or only one end of the chain. In the reflecting boundaries or periodic boundary cases the number of particles will be a conserved quantity provided that no other reactions other than the diffusion of particles take place. In the stationary state, these models exhibit a variety of interesting properties such as non-equilibrium phase transitions and spontaneous symmetry breaking which cannot be found in equilibrium models (see  and references therein). Different applications are also found for such models which include the kinetics of biopolymerization  and traffic flow modelling . These models have also allowed the study of shocks i.e. discontinuities in the density of particles over a microscopic region. Over the last decade people have studied the shocks in one-dimensional driven-diffusive models with open and periodic boundary conditions. A prominent example of such models with open boundary is the Asymmetric Simple Exclusion Process (ASEP) in which the particles enter the system from the left boundary, diffuse in the bulk and leave the chain from the right boundary with certain rates . For specific tuning of injection and extraction rates a shock might appear in the system which moves with a constant velocity towards the boundaries. The ASEP has also been studied on a ring in the presence of a second class particle called the impurity . In this case the impurity will track the shock front with a constant velocity which is determined by the reaction rates of the model. The shocks in the models with reflecting boundaries have not been studied yet.
In the present letter we study the phase transitions in a one-dimensional branching-coalescing model with reflecting boundaries in which the particles diffuse, coagulate and decoagulate on a lattice of length $`L`$. The reaction rules are specifically as follows:

$$\begin{array}{cccc}\text{Diffusion to the left:}\hfill & & \mathrm{}+AA+\mathrm{}\hfill & \text{with rate}q\hfill \\ \text{Diffusion to the right:}\hfill & & A+\mathrm{}\mathrm{}+A\hfill & \text{with rate}q^1\hfill \\ \text{Coalescence to the left:}\hfill & & A+AA+\mathrm{}\hfill & \text{with rate}q\hfill \\ \text{Coalescence to the right:}\hfill & & A+A\mathrm{}+A\hfill & \text{with rate}q^1\hfill \\ \text{Branching to the left:}\hfill & & \mathrm{}+AA+A\hfill & \text{with rate}\mathrm{\Delta }q\hfill \\ \text{Branching to the right:}\hfill & & A+\mathrm{}A+A\hfill & \text{with rate}\mathrm{\Delta }q^1\hfill \end{array}$$
(1)
in which $`A`$ and $`\mathrm{}`$ stand for the presence of a particle and a hole respectively. It is assumed that there is no injection or extraction of particles from the boundaries. We will also assume that the number of particles on the chain is a conserved quantity. This model was first introduced and treated in continuum approximation in . It was then studied using the Empty Interval Method (EIM) in . In this formalism the physical quantities such as the density of particles are calculated from the probabilities to find empty intervals of arbitrary length. Later this model was studied using so-called the Matrix Product Formalism (MPF) . According to this formalism the stationary probability distribution function of the system is written in terms of the products of non-commuting operators $`E`$ and $`D`$ and the vectors $`|V`$ and $`W|`$ as follows

$$P(\tau _1,\mathrm{},\tau _L)=\frac{1}{Z_L}W|\underset{i=1}{\overset{L}{}}(\tau _iD+(1\tau _i)E)|V.$$
(2)
Each site of the lattice is occupied by a particle ($`\tau _i=1`$) or is empty ($`\tau _i=0`$). The factor $`Z_L`$ in (2) is a normalization factor. The operators $`D`$ and $`E`$ stand for the presence of particles and holes respectively and besides the vectors $`|V`$ and $`W|`$ should satisfy the following quadratic algebra

$$\begin{array}{c}[E,\overline{E}]=0\hfill \\ \overline{E}DE\overline{D}=q(1+\mathrm{\Delta })ED\frac{1}{q}DE\frac{1}{q}D^2\hfill \\ \overline{D}ED\overline{E}=qED+\frac{1+\mathrm{\Delta }}{q}DEqD^2\hfill \\ \overline{D}DD\overline{D}=q\mathrm{\Delta }ED\frac{\mathrm{\Delta }}{q}DE+(q+\frac{1}{q})D^2\hfill \\ W|\overline{E}=W|\overline{D}=0\hfill \\ \overline{E}|V=\overline{D}|V=0.\hfill \end{array}$$
(3)
The operators $`\overline{D}`$ and $`\overline{E}`$ are auxiliary operators and do not enter into calculating (2). Having a representation for the quadratic algebra (3) one can easily compute the steady state weights of any configuration of the system using (2). It has been shown that (3) has a four-dimensional representation . For $`q^21+\mathrm{\Delta }`$ we have

$$\begin{array}{c}D=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& \frac{\mathrm{\Delta }}{1+\mathrm{\Delta }}& \frac{\mathrm{\Delta }}{1+\mathrm{\Delta }}& 0\\ 0& 0& \mathrm{\Delta }& 0\\ 0& 0& 0& 0\end{array}\right),E=\left(\begin{array}{cccc}\frac{1}{q^2}& \frac{1}{q^2}& 0& 0\\ 0& \frac{1}{1+\mathrm{\Delta }}& \frac{1}{1+\mathrm{\Delta }}& 0\\ 0& 0& 1& \frac{1}{q^2}\\ 0& 0& 0& \frac{1}{q^2}\end{array}\right),|V=\left(\begin{array}{c}a\\ 0\\ q^2\\ q^21\end{array}\right)\\ W|=(1q^2,1,0,b)\end{array}$$
(4)
in which $`a`$ and $`b`$ are arbitrary constants. The matrix representations for $`\overline{D}`$ and $`\overline{E}`$ are also given in . Both EIM and MPA approaches showed that the model has two different phases: a low-density phase for $`q^2>1+\mathrm{\Delta }`$ and a high-density phase for $`q^2<1+\mathrm{\Delta }`$ and a phase transition takes place at the critical point where $`q^2=1+\mathrm{\Delta }`$. In the low-density and the hight-density regions the density profile of the particles on the chain is an exponential function while on the coexistence line it changes linearly along the lattice.
Recently this model has been studied under the open boundary conditions on a specific manifold of the parameters of the system . It is shown that if the particles are injected and extracted from the left boundary with the rates $`\alpha `$ and $`\beta `$ respectively then the model has the same phase structure provided that $`\alpha =(q^1q+\beta )\mathrm{\Delta }`$. In the latter case the operators $`D`$ and $`E`$ and also the vectors $`W|`$ and $`|V`$ have two-dimensional representations. The only difference is that for the reflecting boundary conditions the system involves three different length scales while for the open boundary conditions it is characterized by one length scale. Moreover, it has been shown that in the open boundary case the probability distribution function of the system can be written in terms of superposition of Bernoulli shock measures . For the open boundary case if we fix the density of particles, for example by working in a canonical ensemble, we can see the real shock structures in the density profile of particles. In this case the system will also have two different phases: a low-density phase and a jammed-phase where the shocks evolve in the system. These phases are specified by the density of particles and are separated by a second-order phase transition .
A natural question that might arise is whether or not we can see the shocks in our branching-coalescing model defined by (1) with the reflecting boundaries. To answer this question we will study the model with reflecting boundaries in a canonical ensemble where the number of particles on the chain is equal to $`M`$ so that the density of particles $`\rho =\frac{M}{L}`$ remains constant. We will then investigate the phase transitions and the density profile of particles on the chain. Recently it has been shown that the classical Yang-Lee theory  can be applied to the one-dimensional out-of-equilibrium systems in order to study the possible phase transitions of these models . According to this theory in the thermodynamic limit, the zeros of the canonical or grand canonical partition function, as a function of an intensive variable of the system, might approach the real positive axis of that parameter at one or more points. Depending on how these zeros approach the real positive axis the system might have one or more phase transition of different orders. If the zeros intersect the real positive axis at a critical point at an angle $`\frac{\pi }{2n}`$, then $`n`$ will be the order of phase transition at that point .
Let us define the canonical partition of our model using the MPF as follows

$$Z_{L,M}=\underset{\{\tau _i=0,1\}}{}\delta (M\underset{i=1}{\overset{L}{}}\tau _i)W|\underset{i=1}{\overset{L}{}}(\tau _iD+(1\tau _i)E)|V$$
(5)
in which $`M`$ and $`L`$ are the number of particles and the length of the system respectively and $`\delta (\mathrm{})`$ is the ordinary Kronecker delta function $`\delta _{x,0}`$. Using the matrix representations $`D`$, $`E`$, $`W|`$ and $`|V`$ given by (4) we have been able to calculate the canonical partition of this model (5) and its zeros numerically. One can use MATHEMATICA to calculate $`W|(E+xD)^L|V`$ for arbitrary $`q`$ and $`\mathrm{\Delta }`$ and finite $`L`$ in which $`x`$ is a free parameter. The result will be a polynomial of $`x`$. The coefficient $`x^M`$ in this polynomial gives the canonical partition function of the system. Formally we can write

$$Z_{L,M}=Coefficient[W|(E+xD)^L|V,M]$$
(6)
in which $`Coefficient[Expr,n]`$ gives the coefficient of $`x^n`$ in the polynomial $`Expr`$. In Fig. 1 we have plotted the numerical estimates for the zeros of $`Z_{L,M}`$ obtained from (6) on the complex-$`q`$ plane for $`L=80`$, $`M=42`$. The canonical partition function (6) has $`4(LM)`$ zeros in the complex-$`q`$ plane. We have found that for large $`L`$ and $`M`$ the locations of these zeros are not sensitive to the value of $`\mathrm{\Delta }`$. We have also calculated the numerical estimates for the roots of (6) as a function of $`\mathrm{\Delta }`$ for fixed values of $`q`$. It turns out that (6), as a function of $`\mathrm{\Delta }`$, does not have any positive root; therefore, we expect that the phase transition points do not depend on $`\mathrm{\Delta }`$.

As can be seen in Fig. 1 the zeros lie on two different curves and accumulate towards two different points on the positive real-$`q`$ axis. By extrapolating the real part of the nearest roots to the positive real-$`q`$ axis for large $`L`$ and $`M`$, we have found that the transition points are $`q_c=\frac{1}{\sqrt{1\rho }}(1<q_c<\mathrm{})`$ and $`q_c^{}=\sqrt{1\rho }(0<q_c^{}<1)`$. As $`\rho 0`$ the two curves lie on each other and we will find only one transition point at $`q_c=q_c^{}=1`$. It appears also that the zeros on both curves approach the real-$`q`$ axis at an angle $`\frac{\pi }{4}`$ (the smaller angle). This predicts two second-order phase transitions at $`q_c`$ and $`q_c^{}`$. The reason that the system has two phase transitions can easily be understood. The parameter $`q`$ determines the asymmetry of the system and for any $`q`$ the system is invariant under the following transformations

$`q`$ $``$ $`q^1`$
$`i`$ $``$ $`Li+1.`$ (7)
Therefore, one can expect to distinguish two critical points which are related according to the symmetry of the system.
Let us now study the density profile of particles on the chain $`\rho (i)`$ in each phase. The density of particles at site $`i`$ is defined as

$$\rho (i)=\frac{_𝒞\tau _iP(\tau _1,\mathrm{},\tau _L)}{_𝒞P(\tau _1,\mathrm{},\tau _L)}$$
(8)
in which $`𝒞`$ is any configuration of the system with fixed number of particles $`M`$ and $`P(\tau _1,\mathrm{},\tau _L)`$ is given by (2). It can be verified that the density profile of the particles $`\rho (i)`$ can be written as

$$\begin{array}{c}\rho (i)=\frac{1}{Z_{L,M}}\\ _{k=0}^{M1}W|Coefficient[C^{i1},k]DCoefficient[C^{Li},Mk1]|V\end{array}$$
(9)
where we have defined $`C:=E+xD`$. Now one can use the matrix representation (4) to calculate (9) using MATHEMATICA. In Fig. 2 we have plotted (9) for two different values of $`q`$ with $`L=60`$ and $`M=36`$.

For this choice of the parameters the transition points are $`q_c=1.581`$ and $`q_c^{}=0.633`$. The density of particle has two general behaviors for $`q>1`$. For $`q>q_c`$ and in the thermodynamic limit ($`L\mathrm{},M\mathrm{},\rho =\frac{M}{L}`$) the density profile of particles is a shock in the bulk of the chain; while in the close vicinity of the left boundary, it increases exponentially. The density of particles in the hight-density region of the shock is equal to $`\rho _{High}=1q^2`$. This region is separated by a rather sharp interface from the low-density region in which the density of particles is equal to $`\rho _{Low}=0`$. The low-density region is extended over $`(1\frac{\rho }{1q^2})L`$ sites. This can be seen in Fig. 2 for $`q=2.5`$; however, the reason that the shock interface is not sharp is that our calculations are not done in real thermodynamic limit. One should expect that the shock front becomes sharper and sharper as the length of the system $`L`$ and also the number of particles on the chain $`M`$ increase. For $`1<q<q_c`$ the density of particles in the bulk of the chain is constant equal to $`\rho `$, it drops near the right boundary exponentially and increases exponentially in the close vicinity of the left boundary. The exponential behavior of the density profile of particles near the boundaries in this phase is due to the finiteness of the representation of the algebra (3). It is known that if the associated quadratic algebra of the model has finite dimensional representations, the density-density correlation functions cannot have algebraic behaviors . At $`q=1`$ one finds $`\rho (i)=\rho `$. The density profile of particles in the region $`q<1`$ is related to that of $`q>1`$ through (Shock in a Branching-Coalescing Model with Reflecting Boundaries) that is

$$\rho (i,q)=\rho (Li+1,q^1).$$
(10)
One can also study this model on a ring of length $`L`$ with periodic boundary conditions. Let us assume that the number of particles and therefore their density fluctuates and is not a constant. In this case the probability of finding the system in a specific configuration $`𝒞=\{\tau _1,\mathrm{},\tau _L\}`$ should be obtained from

$$P(\tau _1,\mathrm{},\tau _L)=\frac{1}{Z_L}tr[\underset{i=1}{\overset{L}{}}(\tau _iD+(1\tau _i)E)].$$
(11)
in which $`tr[\mathrm{}]`$ is the trace of the products of matrices. The normalization factor $`Z_L`$ can be obtained from the fact that $`_𝒞P(𝒞)=1`$. This function plays the role of the grand canonical partition function of the system. One can easily check that the quadratic algebra associated to the periodic boundary condition case is the same as (3) except the boundary terms which contain the vectors $`|V`$ and $`W|`$; therefore, we can still use the algebra (3) and its representation (4) to calculate (11). As we mentioned the grand canonical partition function of the system can be defined as

$$Z_L=\underset{\{\tau =0,1\}}{}tr[\underset{i=1}{\overset{L}{}}(\tau _iD+(1\tau _i)E)]=tr[(E+D)^L].$$
(12)
This can easily be calculated and we find

$$Z_L=1+q^{2L}+q^{2L}+(1+\mathrm{\Delta })^L.$$
(13)
The study of the zeros of (13) as a function of $`q`$ in the thermodynamic limit $`L\mathrm{}`$ shows that they approach the positive real-$`q`$ axis at two different points $`q_c=\sqrt{1+\mathrm{\Delta }}(1<q_c<\mathrm{})`$ and $`q_c^{}=\frac{1}{\sqrt{1+\mathrm{\Delta }}}(0<q_c^{}<1)`$. Moreover, the zeros approach the positive real-$`q`$ axis at angle $`\frac{\pi }{2}`$; therefore, unlike the reflecting boundaries case both phase transitions in this case are of first-order. In the steady state the density profile of the particles on the ring is flat and equal to $`\rho (i)=\frac{\mathrm{\Delta }}{1+\mathrm{\Delta }}`$ for $`q_c^{}<q<q_c`$ and $`\rho (i)=0`$ for $`q>q_c`$ and $`q<q_c^{}`$. One can take the total density of particles on the ring as an order parameter and since it changes discontinuously over the transition points, as the Yang-Lee theory predicts, the phase transitions are of first-order. One should note that in comparison to the reflecting boundaries case not only the geometry of the model is changed but also the number of particles in this case in not a conserved quantity. This means that we should not expect the transitions in these models to be similar even if we take the thermodynamic limit $`L\mathrm{}`$. The reason for the existence of two phase transition points is again the symmetry of the model.
In this letter we have studied a branching-coalescing model in which particles hop, coagulate and decoagulate on one-dimensional lattice of length $`L`$. We have restricted ourself to the case where the total number of particles on the chain is constant. For this we have worked in the canonical ensemble in which the number of particles $`M`$ is fixed; therefore, the density of particles on the chain $`\rho =\frac{M}{L}`$ is constant. The Yang-Lee theory predicts that the model has two second-order phase transitions. Both phase transition points are determined by the density of particles on the system $`\rho `$. The study of the mean particle concentration at each site of the chain for $`q>1`$ shows that the density profile of the particles has a shock-like structure in the region $`q>q_c=\frac{1}{\sqrt{1\rho }}`$. The exception is near the left boundary where the density of particles increases exponentially. This is the first time that shocks are seen in one-dimensional reaction-diffusion models with reflecting boundaries. In the region $`1<q<q_c=\frac{1}{\sqrt{1\rho }}`$ the density profile of the particles is constant in the bulk of the chain; however, near the left (right) boundary it increases (decreases) exponentially. The exact form of the correlation lengths is under investigation. Our numerical investigations also show that the width of the shock scales as $`L^\nu `$ with $`\nu =\frac{1}{2}`$. In the thermodynamic limit $`L\mathrm{}`$ the shock width goes to zero and one finds a very sharp shock interface. Since the system is invariant under the transformation (Shock in a Branching-Coalescing Model with Reflecting Boundaries), the density profile of the particles for $`q<1`$ can be obtained from (10). We also studied the periodic boundary case and found that the system possess two first-order phase transitions which are determined by the values of $`\mathrm{\Delta }`$ for fixed value of $`q`$. In this case the density profile of particles is flat everywhere on the lattice and is either equal to $`\frac{\mathrm{\Delta }}{1+\mathrm{\Delta }}`$ or zero. The formulas (6) and (9) provide us with a simple and general tool for numerical study of the phase transitions and also the particle concentration behaviors of one-dimensional stochastic models in canonical ensemble for which a finite- or infinite-dimensional representation of the associate algebra exists.