format: text

text:
# Liquid-gas phase behaviour of an argon-like fluid modelled by the hard-core two-Yukawa potential

## I Introduction

Much attention has been paid in recent years to the hard core Yukawa (HCY) potential as a model for the pair interactions of fluids . Interest in the potential is motivated on the one hand by its relevance to solvent averaged interactions in polyelectrolytic and colloidal particles and, on the other hand, by its analytical tractability in the context of liquid state theories such as the Mean Spherical Approximation (MSA) and the Self Consistent Ornstein-Zernike Approximation (SCOZA). Recent studies of the HCY fluid can be found in  and references therein.

For simple fluids interacting via dispersion forces, however, the bare HCY potential fails to provide a realistic representation of the interactions, which are much better modelled by the Lennard-Jones (LJ) potential. Unfortunately owing to its mathematical structure, the latter is less amenable to direct study by the MSA and SCOZA than is the HCY potential. Nevertheless, it has long been appreciated that a LJ-like simple-fluid potential can be replaced by a potential with a hard core plus a linear combination of two Yukawa tails  with no appreciable loss of agreement between the experimental and model equations of state. Such a potential permits the theoretical study of simple fluids whilst retaining the convenient mathematical properties of the HCY potential.

Several different representations of a LJ-like fluid in terms of a hard-core plus two Yukawa fluid (HC2YF) have been proposed in the literature . The one we use here is very similar to that recently given by Kalyuzhnyi and Cummings :

$$v(r)=\{\begin{array}{cc}\mathrm{}\hfill & r\sigma ,\hfill \\ & \\ \frac{A_1ϵ}{r}\mathrm{exp}[z_1(r\sigma )]\frac{A_2ϵ}{r}\mathrm{exp}[z_2(r\sigma )]\hfill & r>\sigma .\hfill \end{array}$$
(1)
Here $`\sigma `$ and $`ϵ`$ set the length and energy scales of the model and coincide with the zero and the well depth of the LJ potential respectively. The parameters $`A_1=1.6438\sigma `$, $`z_1=14.7\sigma ^1`$, $`A_2=2.03\sigma `$ and $`z_2=2.69\sigma ^1`$ are chosen to fit the LJ potential and match the second virial coefficient to that of the LJ fluid. Eq. (1) differs slightly from the original parameterisation by Kalyuzhnyi and Cummings insomuch as these authors determined the hard-sphere diameter $`\sigma _{\mathrm{HS}}`$ as a function of temperature via the Barker-Henderson procedure  in order to represent the soft-core repulsive interaction which results from the Weeks-Chandler-Henderson separation of the LJ potential . We instead hold the hard-core diameter fixed by setting $`\sigma _{\mathrm{HS}}=\sigma `$ for conceptual simplicity and clarity in making comparison with our simulation results. A comparison of Eq. (1) with the LJ potential is shown in Fig. 1. Kalyuzhnyi and Cummings  have used the MSA, Percus-Yevick approximation and reference hypernetted chain theories to obtain the liquid-gas phase behaviour of the HC2YF specified by their parameterisation. Their results showed good agreement with existing literature data for the LJ fluid.

In the present work, we extend the work of Kalyuzhnyi and Cummings  by determining the liquid-gas phase behaviour of the HC2YF using the SCOZA of Høye and Stell. Unlike those considered in ref. , this theory is neither mean-field nor mean-spherical like and on the basis of a similar study already performed for the HCY fluid, it can be expected to give a superior account of fluid properties in the near-critical regime. We compare the predictions of the SCOZA for the the critical point and the coexistence curve to those obtained by other approaches, including the renormalization-group based hierarchical reference theory (HRT) of Parola and Reatto , and to new high resolution simulation results for the HC2YF. We find that one version of the SCOZA considered provides a remarkably accurate critical point and coexistence curve. The critical density and temperature predicted by the theory agree with the simulation results to within 1%.

## II Theory

The SCOZA deals with two-body potentials which, like that of Eq. (1), consist of a singular hard-sphere repulsion with diameter $`\sigma `$ and a longer-ranged tail $`w(r)`$. As is customary in integral-equation theories, this approach is based upon the Ornstein-Zernike (OZ) equation linking the two-body radial distribution function $`g(r)`$ to the direct correlation function $`c(r)`$. A closed theory is obtained by supplementing the OZ equation with an approximate relation involving $`g(r)`$ and $`c(r)`$. In its simplest form, the SCOZA amounts to setting:

$$\{\begin{array}{cc}g(r)=0\hfill & r<1,\hfill \\ & \\ c(r)=K(\rho ,\beta )w(r)\text{ }\hfill & r>1,\hfill \end{array}$$
(2)
where $`\rho `$ is the number density of the system, $`\beta =1/(k_\mathrm{B}T)`$ is the inverse temperature, and the hard-sphere diameter has been set equal to one. This closure resembles that adopted in the MSA, except that the amplitude $`K`$ of the direct correlation function outside the repulsive core is regarded as an unknown state-dependent quantity, to be determined in such a way that consistency between the compressibility and the energy route to thermodynamics is enforced. This constraint amounts to requiring that the reduced compressibility $`\chi _{\mathrm{red}}`$ and the excess internal energy per unit volume $`u`$ satisfy the condition:

$$\frac{}{\beta }\left(\frac{1}{\chi _{\mathrm{red}}}\right)=\rho \frac{^2u}{\rho ^2},$$
(3)
where it is understood that $`\chi _{\mathrm{red}}`$ is obtained by the compressibility sum rule as the structure factor $`S(k)`$ evaluated at $`k=0`$, while $`u`$ is obtained by the energy equation as the spatial integral of the tail interaction $`w(r)`$ weighted by the radial distribution function $`g(r)`$. If the closure (2) is used for the correlations, the consistency condition (3) yields a closed partial differential equation (PDE) for the function $`K(\rho ,\beta )`$. In the present case, implementing this scheme is made simpler by taking advantage of the analytical results obtained for the OZ equation with the closure (2) when the tail potential $`w(r)`$ has a two-Yukawa form like that of Eq. (1). These enable one to obtain $`\chi _{\mathrm{red}}`$ as a function of $`\rho `$ and $`u`$. The function $`\chi _{\mathrm{red}}(\rho ,u)`$ can then be used in Eq. (3) by taking $`u`$ instead of $`K`$ as the unknown quantity. The algebraic manipulations are similar (although not identical) to those we performed in Ref. , and are based on the results obtained in Ref. . We start by introducing the quantity:

$$f=(1\xi )\sqrt{\frac{1}{\chi _{\mathrm{red}}}},$$
(4)
where $`\xi =\pi \rho /6`$ is the packing fraction. This can be written as:

$$f=\frac{(z_1^2z_2^2)+4\sqrt{q}(\gamma _2\gamma _1)}{4\left[(z_1/z_2)\gamma _2(z_2/z_1)\gamma _1\right]}\frac{z_1^2z_2^2}{z_1z_2}\frac{\gamma _1\gamma _2(\gamma _2\gamma _1)}{\left[(z_1/z_2)\gamma _2(z_2/z_1)\gamma _1\right]^2},$$
(5)
where $`q=(1+2\xi )^2/(1\xi )^2`$, and $`\gamma _1`$, $`\gamma _2`$ are state-dependent functions. By using Eqs. (4) and (5), Eq. (3) becomes:

$$\frac{2f}{(1\xi )^2}\left[\left(\frac{f}{\gamma _1}\right)_{\rho ,\gamma _2}\left(\frac{\gamma _1}{u}\right)_\rho +\left(\frac{f}{\gamma _2}\right)_{\rho ,\gamma _1}\left(\frac{\gamma _2}{u}\right)_\rho \right]\left(\frac{u}{\beta }\right)_\rho =\rho \left(\frac{^2u}{\rho ^2}\right)_\beta ,$$
(6)
where the partial derivatives of $`f`$ with respect to $`\gamma _1`$, $`\gamma _2`$ are straightforwardly obtained from Eq. (5). In order to obtain a closed PDE for $`u`$, one then needs to express $`\gamma _1`$, $`\gamma _2`$ as a function of $`\rho `$ and $`u`$. To this end, we recall the definitions of $`\gamma _1`$, $`\gamma _2`$ :

$$\gamma _i=2\sqrt{q}\frac{4+2z_iz_i^2}{2(2+z_i)}\frac{\tau _iI_i1}{\sigma _iI_i1}\text{ }(i=1,2),$$
(7)
where the quantities $`\sigma _i`$, $`\tau _i`$ are known and depend only on the inverse range $`z_i`$, and $`I_i`$ is given by:

$$I_i=4\pi \rho _1^+\mathrm{}𝑑rr\mathrm{exp}[z_i(r1)]g(r)\text{ }(i=1,2).$$
(8)
For the two-Yukawa tail potential of Eq. (1), the integrals $`I_1`$, $`I_2`$ are straightforwardly related to the internal energy per unit volume $`u`$:

$$u=\frac{1}{2}\rho ϵ(A_1I_1A_2I_2).$$
(9)
If Eq. (7) is used to express $`I_i`$ as a function of $`\gamma _i`$ and the result is substituted into Eq. (9), one obtains an expression for $`u`$ as a function of $`\rho `$, $`\gamma _1`$, $`\gamma _2`$. This can be inverted algebrically to obtain $`\gamma _2`$ as a function of $`\rho `$, $`\gamma _1`$, and $`u`$. One is then left with the task of expressing $`\gamma _1`$ in terms of $`\rho `$ and $`u`$. If we set:

$`x`$ $`=`$ $`\sqrt{q}{\displaystyle \frac{z_1^2}{4\gamma _1}},`$ (10)
$`y`$ $`=`$ $`\sqrt{q}{\displaystyle \frac{z_2^2}{4\gamma _2}},`$ (11)
we find by applying the results of Ref.  that the following equation hold:

$`A_2z_2^4\sigma _1^2(2+z_1)^2[4(2\sqrt{q}\alpha _1)(\sqrt{q}x)z_1^2]^2\{4(z_2^24y^2)(yx)^2`$ (13)
$`(z_1^2z_2^2)[z_1^2z_2^2+4(y^2x^2)]\}`$
$``$ $`A_1z_1^4\sigma _2^2(2+z_2)^2[4(2\sqrt{q}\alpha _2)(\sqrt{q}y)z_2^2]^2\{4(z_1^24x^2)(yx)^2`$ (15)
$`(z_1^2z_2^2)[z_1^2z_2^2+4(y^2x^2)]\}=0,`$
where $`\alpha _i`$, $`i=1,2`$, are known functions of the inverse-range parameters $`z_i`$. Note that this equation does not contain the unknown amplitude $`K(\rho ,\beta )`$. By writing $`x`$ and $`y`$ as functions of $`\gamma _1`$, $`\gamma _2`$ via Eqs. (10), (11), and subsequentely $`\gamma _2`$ as a function of $`\rho `$, $`\gamma _1`$, $`u`$, we finally obtain an equation for $`\gamma _1`$ as a function of $`\rho `$ and $`u`$. If we indicate the l.h.s. of Eq. (15) by $`F(x,y,\rho )`$, the equation for $`\gamma _1`$ has the form:

$$F\{x(\gamma _1),y[\gamma _2(\rho ,\gamma _1,u)],\rho \}=0.$$
(16)
Eq. (6) then becomes:

$$B(\rho ,u)\left(\frac{u}{\beta }\right)_\rho =C(\rho ,u)\left(\frac{^2u}{\rho ^2}\right)_\beta ,$$
(17)
where $`B(\rho ,u)`$ and $`C(\rho ,u)`$ are given by:

$`B(\rho ,u)`$ $`=`$ $`{\displaystyle \frac{2f}{(1\xi )^2}}{\displaystyle \frac{\gamma _2}{u}}\left[{\displaystyle \frac{f}{\gamma _2}}{\displaystyle \frac{F}{x}}{\displaystyle \frac{x}{\gamma _1}}{\displaystyle \frac{f}{\gamma _1}}{\displaystyle \frac{F}{y}}{\displaystyle \frac{y}{\gamma _2}}\right],`$ (18)
$`C(\rho ,u)`$ $`=`$ $`\rho \left[{\displaystyle \frac{F}{y}}{\displaystyle \frac{y}{\gamma _2}}{\displaystyle \frac{\gamma _2}{\gamma _1}}+{\displaystyle \frac{F}{x}}{\displaystyle \frac{x}{\gamma _1}}\right],`$ (19)
where all the partial derivatives are performed at constant density $`\rho `$.

We note that in Eq. (2) there is no hard-sphere contribution to the direct correlation function outside the repulsive core. As a consequence, the treatment of thermodynamics and correlations of the hard-sphere gas that come out of Eq. (2) in the high-temperature limit coincide with that of the Percus-Yevick (PY) integral equation. An improved description of the hard-sphere thermodynamics at the level of the well known Carnahan-Starling (CS) equation of state is desirable in this context, as the slight inaccuracy in the PY treatment of the hard-sphere gas affects not only the high-density behavior of the system in study, but also the location of its critical point and of the coexistence curve. A non-vanishing contribution to $`c(r)`$ for $`r>1`$ due to hard-core part of the interaction could be taken into account by replacing the expression for $`c(r)`$ of Eq. (2) with

$$c(r)=c_{\mathrm{HS}}(r)+K(\rho ,\beta )w(r)\text{ }r>1,$$
(20)
where the hard-sphere direct correlation function $`c_{\mathrm{HS}}(r)`$ can be determined, for instance, by the Waisman parameterisation . This is what has been done in Ref.  for the HCY fluid. However, in the present case of a two-Yukawa tail potential, this would require one to deal with a direct correlation function of three-Yukawa form for $`r>1`$. Although such an extension of our treatment appears to be feasible , in order to minimize the complexity of the computation, we have not pursued it here. Instead, we assumed for $`c_{\mathrm{HS}}(r)`$ outside the core a Yukawa form whose range coincides with that of the repulsive contribution to the tail potential, and whose density-dependent amplitude $`H`$ is set so as to give CS thermodynamics in the high-temperature limit. This procedure does not aim at the most accurate description of the hard-core contribution to the correlations possible in the context of the SCOZA. But it has the advantage of taking into account the hard-sphere thermodynamics beyond the PY level without going beyond the two-Yukawa form of $`c(r)`$. An analogous procedure involving only a single Yukawa term was already implemented for the HCY fluid , where it was found to give results very similar to those of the more reliable closure (20) both for the thermodynamics and the phase diagram. According to this prescription, for $`r>1`$, $`c(r)`$ has the form:

$$c(r)=\frac{(H+KA_1)}{r}\mathrm{exp}[z_1(r1)]\frac{KA_2}{r}\mathrm{exp}[z_2(r1)]\text{ }r>1,$$
(21)
while for $`r<1`$ the core condition $`g(r)=0`$ holds as before. Here $`K`$ is the unknown, state-dependent amplitude of Eq. (2), which vanishes in the high-temperature limit, while $`H`$ is a known function of the density. As the inverse range of the hard-sphere contribution to $`c(r)`$ is locked to $`z_1`$, the treatment of the hard-sphere gas that comes out of Eq. (21) clearly lacks the virial-compressibility consistency of the Waisman parameterisation. However, this does not affect the consistency between internal energy and compressibility route upon which SCOZA hinges. The manipulations that lead to the PDE (17) remain unchanged, the only difference being that Eq. (15) will contain an extra term related to $`H`$, which will affect the partial derivatives $`F/x`$, $`F/y`$ that appear in Eqs. (18), (19).

The PDE (17) has been integrated numerically. The initial condition at $`\beta =0`$ and the boundary conditions at the ends of the density interval are the same as in Ref.  and will not be detailed here. Once $`u`$ has been obtained by solving Eq. (17), integration of $`u`$ with respect to $`\beta `$ yields the Helmholtz free energy and hence all the other thermodynamic quantities.

## III Simulation procedure

The principal aspects of the simulation and finite-size scaling techniques employed in this work have previously been detailed elsewhere in the context of a similar study of the Lennard-Jones fluid. Accordingly we restrict ourself to a brief summary of the methodology and refer the reader to reference  for a fuller account.

Grand canonical Monte-Carlo (MC) simulations were performed for the HC2YF model of Eq. (1). The algorithm used had a Metropolis form  and comprised only particle transfer (insertion and deletion) steps, leaving particle moves to be performed implicitly as a result of repeated transfers. The potential was cut at a radius $`r_c=3.0\sigma `$, and a standard correction term was applied to the internal energy to compensate for the trunction . To simplify identification of particle interactions the periodic simulation space of volume $`L^3`$ was partitioned into $`m^3`$ cubic cells, each of side the cutoff $`r_c`$. This strategy ensures that interactions emanating from particles in a given cell extend at most to particles in the $`26`$ neighbouring cells.

System sizes having $`m=3,4,5,6,7`$ and $`8`$ were studied, corresponding (at coexistence) to average particle numbers of approximately $`230,540,1050,1750,2900`$ and $`4500`$ respectively. For the $`m=3,4,5`$ and $`6`$ system sizes, equilibration periods of $`10^5`$ Monte Carlo transfer attempts per cell (MCS) were utilised, while for the $`m=7`$ and $`m=8`$ system sizes up to $`2\times 10^6`$ MCS were employed. Sampling frequencies ranged from $`20`$ MCS for the $`m=3`$ system to $`250`$ MCS for the $`m=8`$ system. The total length of the production runs was also dependent upon the system size. For the $`m=3`$ system size, $`1\times 10^7`$ MCS were employed, while for the $`m=8`$ system, runs of up to $`1\times 10^8`$ MCS were necessary.

In the course of the simulations, the observables recorded were the particle number density $`\rho =N/V`$ and the energy density $`u=E/V`$. The joint distribution $`p_L(\rho ,u)`$ was accumulated in the form of a histogram. In accordance with convention , we express all thermodynamic quantities in reduced units: $`\rho ^{}=\rho \sigma ^3`$, $`u^{}=u\sigma ^3/ϵ`$, $`T^{}=k_\mathrm{B}T/ϵ`$.

Efficient exploration of the phase space was facilitated through use of the histogram reweighting technique . This method allows histogram accumulated at one set of model parameters to be reweighted to provide estimates appropriate to another set of not-too-distant model parameters. Use of the method permits large areas of phase space to be mapped using only a few simulations performed at strategic state points. To facilitate study of the subcritical coexistence region, the multicanonical preweighting technique  was employed. This method employs a biased sampling technique to overcome the free energy barrier separating the coexisting phases and thus allow both to be sampled in a single simulation run. When combined with histogram reweighting in the manner described in ref. , multicanonical preweighting permits an extremely efficient accumulation of coexistence curve data.

The task of estimating the critical point parameters of the model was performed using finite-size scaling techniques . In brief, the strategy is to match the measured “ordering operator distribution” to an independently known universal critical point form appropriate to the Ising universality class. The ordering operator itself is defined as $`(\rho ^{}+su^{})`$, where $`s`$ is a non-universal “field mixing” parameter, which is finite in the absence of particle-hole symmetry, and which is chosen to ensure that $`p()`$ is symmetric in $``$. For sufficiently large $`L`$, the matching should occur at the critical point parameters. For small $`L`$, however, the apparent critical temperature obtained by this matching procedure is subject to systematic errors associated with corrections to finite-size scaling. To deal with this, we extrapolate to the thermodynamic limit using the known scaling properties of the corrections, which are expected to diminish (for sufficiently large system sizes) like $`L^{\theta /\nu }`$ , where $`\theta `$ is the correction to scaling exponent and $`\nu `$ is the correlation length exponent. The extrapolation has been performed using a least squares fit to the data for the four largest system sizes. The results of the extrapolation are shown in Fig. 2, from which we estimate $`T_c^{}=1.295(10)`$. The associated estimate for the critical density is $`\rho _c^{}=0.310(1)`$ and for the reduced chemical potential is $`\mu _c^{}=3.588(30)`$.

In addition to the phase coexistence data, we have also measured the form of the radial distribution function $`g(r)`$ for a number of state points, corresponding to reduced temperatures $`T^{}=2`$, $`1.5`$, $`1`$, and reduced densities $`\rho ^{}=0.4`$, $`0.6`$, $`0.8`$. In the following Section, the simulation results for the coexistence curve, critical point properties and radial distribution function of the HC2YF are compared with the theoretical predictions.

## IV Results and Discussion

The SCOZA coexistence curve has been determined by equating the pressure $`P`$ and the chemical potential $`\mu `$ on the low- and high-density branch of the subcritical isotherms. As we specified above, these quantities were obtained via the energy route by integrating $`u`$ with respect to $`\beta `$. The advantage of doing so is that one does not have to circumvent the forbidden region bounded by the spinodal curve, i.e. the locus of diverging compressibility, in order to obtain $`P`$ and $`\mu `$ on the high-density branch. On the other hand, because of the compressibility–internal energy consistency of the theory, this is fully equivalent to using a mixed path combining integration of the inverse compressibility $`1/\chi _{\mathrm{red}}`$ with respect to $`\rho `$ and integration of $`u`$ with respect to $`\beta `$, as is often done in calculations based on integral equations . The coexistence curve in the density-temperature and in the temperature-chemical potential plane is shown in Figs. 3 and 4 respectively. Reduced units have been used throughout. The SCOZA results obtained both by the closure (2) and by the modified version (21) are compared with the MC data obtained in this work. We have also plotted the coexistence curve predicted by HRT and by the energy route of the lowest-order gamma-ordered approximation (LOGA) , also known as optimized random-phase approximation (ORPA) , together with the spinodal curve given by LOGA/ORPA compressibility route. It appears that the SCOZA yields the most accurate determination of the coexistence curve among the theories considered here, provided Eq. (21) is used in order to describe hard-sphere thermodynamics at the CS level. The predictions for the critical point are compared in Tab. I. The error on the SCOZA critical density and temperature is below $`1\%`$. Actually, the critical temperature predicted by SCOZA with closure (2) and PY hard-sphere thermodynamics is even closer to the MC result, but this is most likely to be accidental since, as stated in Sec. II, on the basis of our previous calculations on the HCY potential  we expect that a treatment of the HC2YF based on closure (20) and a $`c(r)`$ of three-Yukawa form will in fact deliver results almost identical to those found here by Eq. (21). We note that, like the HRT or the energy-route LOGA/ORPA, the SCOZA gives a coexistence curve that goes right up to the critical point. This feature is not shared by other integral-equation approaches, such as the MHNC or the HMSA theories . Both the SCOZA below the critical temperature and the HRT yield non-classical critical exponents; for the exponent $`\beta _{\mathrm{coex}}`$ which describes the curvature of the coexistence curve the SCOZA gives $`\beta _{\mathrm{coex}}=7/20=0.35`$ , while according the HRT $`\beta _{\mathrm{coex}}0.345`$ , the best theoretical estimate being $`\beta _{\mathrm{coex}}0.327`$ . Moreover, in the SCOZA, as well as in the HRT, the critical point, identified as the top of the coexistence curve, coincides with that determined by locating the divergence of the isothermal compressibility as given by the compressibility sum rule. As a consequence, at the critical point one has both coalescence of the vapor and liquid phases and occurrence of long-range correlations, as expected. Because of the lack of thermodynamic consistency, this is not the case with LOGA/ORPA, in which the critical point obtained by the compressibility route corresponds to the top of the spinodal curve shown in Fig. 3.

We now consider the radial distribution function $`g(r)`$. The SCOZA $`g(r)`$ has been determined using both closure (2) and (21). In Fig. 5 the SCOZA and LOGA/ORPA results for two different states are compared with MC simulations. We see that, while the overall agreement is satisfactory, the behavior near contact is not reproduced very well. It is unlikely that this can be traced back to the slight inaccuracy in the treatment of the hard-sphere contribution to $`c(r)`$ for $`r>1`$ entailed by Eq. (2) or (21). In fact, such a contribution is accurately taken into account in LOGA/ORPA, which nevertheless does not appear to perform better than SCOZA. In particular, using closure (21) in SCOZA produces results for $`g(r)`$ which are undistinguishable from those of LOGA/ORPA, at least for the states we investigated. Moreover, we checked that for the densities considered here the hard-sphere $`g(r)`$ obtained from closure (21) in the high-temperature limit is nearly superimposed to that given by the Waisman parameterisation. On the other hand, both in SCOZA and in LOGA/ORPA the contribution to $`c(r)`$ due to the tail potential $`w(r)`$ is bound to follow the profile of $`w(r)`$ itself for all $`r`$’s. Such a form is best suited for interactions that are slowly varying on a lengthscale of the order of the particle size, while the two-Yukawa tail of Eq. (1) changes quite steeply near its minimum. As is well known, such a problem is usually dealt with by splitting the potential according to the WCA prescription , but implementing this procedure in the SCOZA seems very artificial, as it would hinder the analytical tractability of the system, which is the very reason why the HC2YF interaction (1) was considered in the first place instead of the LJ potential.

Since on the basis of the comparison with MC results the SCOZA proves to accurately reproduce the coexistence curve of the HC2YF, it could be worthwhile using it to assess the ability of different HC2YF parameterisations of the LJ potential to reproduce the LJ coexistence curve. To this end, in Fig. 6 we have compared the coexistence curve of the LJ potential given by MC simulations  with the SCOZA results using closure (21) for the parameterisations proposed by Sun  and by Foiles and Ashcroft , together with that of Eq. (1). We see that the HC2YF form of Eq. (1), closely resembling that of Ref. , gives the best representation of the LJ potential among those considered here, at least as far as the liquid-vapor phase diagram is concerned.

In summary, we have presented SCOZA results for the phase diagram and the correlations of a two-Yukawa parameterisation of the LJ potential. These have been compared with the predictions of other theories and with new MC simulation data. The comparison shows that the SCOZA provides a very accurate coexistence curve, with a critical point that differs from the prediction of MC supplemented by finite-size scaling techniques by less than $`1\%`$. At the same time, the two-Yukawa potential we studied, closely resembling that previously proposed by Kalyuzhnyi and Cummings , reproduces satisfactorily the LJ coexistence curve. In view of the very good performance of SCOZA in predicting the liquid-vapor phase diagram of simple fluids, we think that it could be worthwhile to generalize the solution procedure, so that it will not be confined anymore to tail potentials of Yukawa form. This could also allow one to implement SCOZA with a more sophisticated closure than the forms (2) and (21) used here.

###### Acknowledgements.

N.B.W. thanks the Royal Society (grant number 19076), the Royal Society of Edinburgh and the EPSRC (grant no. GR/L91412) for financial support. D.P. thanks the National Science Foundation for financial support while visiting Stony Brook. G.S. gratefully acknowledges the support of the Division of Chemical Sciences, Office of Basic Energy Sciences, Office of Energy Research, U.S. Department of Energy.


text:
# 1 The assumed ranges for 𝑎_{𝑑𝑖𝑟}⁰⁰, 𝐵⁰⁰/𝐵⁺⁻ and 𝐵⁺⁰/𝐵⁺⁻ for five (hypothetical) sets of experimental measurements.

## Acknowledgments

We thank the organizers of BCP4 for an exciting, fascinating conference. The work of D.L. was financially supported by NSERC of Canada.

## References


text:
# Relativistic close coupling calculations for photoionization and recombination of Ne-like Fe XVII

## I INTRODUCTION

Laboratory, astrophysical, and theoretical studies of Fe XVII are of considerable interest as it is a prime constituent in high temperature plasmas responsible for strong X-ray emission . A number of atomic processes need to be considered in detail, primarily: electron impact excitation, photoionization, (e + ion) recombination, and radiative transitions. Large-scale atomic calculations are in progress for all of these processes in Fe XVII under the Iron Project and related works , in extended energy ranges suitable for practical applications. While electron impact excitation is an independent part of this effort, in this work we describe photoionization and (e + ion) recombination of Fe XVII. The coupled channel approximation including relativistic effects for many channel systems can be very involved owing to many infinite series of resonance structures converging on to the various excited levels of the core ion. Whereas the relativistic and coupling effects have been studied previously, all such theoretical and experimental studies of Fe XVII photoionization (e.g. Haque et al.  and Mohan et al. ) have been limited to the ground state and the relatively small energy range spanned by core excitations within the $`n`$ = 2 complex of the residual ion Fe XVIII comprising of 3 fine structure levels up to about 132 eV, i.e.

$$h\nu +\mathrm{Fe}\mathrm{XVII}(2s^22p^6{}_{}{}^{1}S_{0}^{})e+\mathrm{Fe}\mathrm{XVIII}(2s^22p^5{}_{}{}^{2}P_{3/2,1/2}^{o}),2s2p^6(^2S_{1/2}).$$
(1)
Although the near-threshold behavior of photoionization and recombination cross sections (Pradhan et al. ) is physically interesting, it is inadequate for practical applications that require the cross sections to be calculated up to high energies typical of the variety of conditions where L-shell ions are abundant. Purely photoionized plasmas, such as in H II regions, planetary nebulae, novae etc., are typically low temperature, but coronal plasmas cover a much wider range . For example, the temperature of maximum abundance of Fe XVII in the coronal ionization equilibrium is about 4 $`\times 10^6`$ K . Furthermore, in astrophysical objects such as the warm absorber ionized gas thought to surround the central black hole in active galactic nuclei, the plasma is most likely of a composite nature since most ionization states of several elements are observed (e.g. ). High accuracy throughout the energy range of practical importance is therefore essential.

As the (e + ion) recombination is unified in nature, it is theoretically desirable to consider the non-resonant and resonant processes (RR and DR) together. A unified theoretical formulation has been developed  including relativistic fine structure , and used to compute cross sections and rates for many atomic systems, such as the K-shell systems C IV – C V and Fe XXIV – Fe XXV of interest in X-ray spectroscopy . The unified results may be directly compared with experimental results, without the need to separate RR and DR. In this paper we present details of the low-energy results for Fe XVII and show that not only are the unified cross sections and rates in good agreement with experiment, but also illustrate how the unified calculations avoid the basic inconsistency and incompleteness of photoionization and recombination data for modeling of laboratory and astrophysical plasma sources.

The present report describes in detail the 3-level and the 60-level close coupling calculations, with a discussion of the relativistic and coupling effects and comparison with earlier 3-level theoretical and experimental data for photoionization and recombination. While photoionization and recombination are usually considered separately, we exemplify and emphasize the underlying physical unity, via detailed balance, between the two processes as naturally treated in the close coupling method.

## II THEORY

Photoionization and (e + ion) recombination may both be considered using identical coupled channel wavefunction expansion,

$$\mathrm{\Psi }(E)=A\underset{i}{}\chi _i\theta _i+\underset{j}{}c_j\mathrm{\Phi }_j,$$
(2)
where $`\mathrm{\Psi }`$ represents a ($`N`$+1)-electron bound or continuum state depending on $`E<`$0 or $`E>`$0, expressed in terms of the $`N`$-electron residual core ion eigenfunctions. The $`\chi _i`$ is the target wavefunction in a specific state $`S_iL_i\pi _i`$ or $`J_i\pi _i`$, and $`\theta _i`$ is the wavefunction for the ($`N`$+1)-th electron in a channel labeled as $`S_iL_i(J_i)\pi _ik_i^2\mathrm{}_i(SL\pi \mathrm{or}J\pi ))`$; $`k_i^2`$ being its incident kinetic energy. $`\mathrm{\Phi }_j`$’s are the correlation functions of the ($`N`$+1)-electron system that account for short range correlation and the orthogonality between the continuum and the bound orbitals. The R-matrix method , and its relativistic Breit-Pauli extension , enables a solution for the total $`\mathrm{\Psi }`$, with a suitable expansion over the $`\chi _i`$. The Breit-Pauli R-matrix (BPRM) method has been extensively employed for electron impact excitation under the Iron Project . The extension of the BPRM formulation to unified electronic recombination , and theoretically self-consistent calculations of photoionization and recombination is sketched below.

Recombination of an incoming electron to the target ion may occur through non-resonant, background continuum, usually referred to as radiative recombination (RR),

$$e+X^{++}h\nu +X^+,$$
(3)
which is the inverse process of direct photoionization, or through the two-step recombination process via autoionizing resonances, i.e. dielectronic recombination (DR):

$$e+X^{++}(X^+)^{}\{\begin{array}{c}e+X^{++}\\ h\nu +X^+\end{array},$$
(4)
where the incident electron is in a quasi-bound doubly-excited state which leads either to (i) autoionization, a radiation-less transition to a lower state of the ion and the free electron, or to (ii) radiative stabilization predominantly via decay of the ion core, usually to the ground state, and the bound electron.

In the unified treatment the photoionization cross sections, $`\sigma _{\mathrm{PI}}`$, of a large number of low-$`n`$ bound states – all possible states with $`nn_{\mathrm{max}}10`$ – are obtained in the close coupling (CC) approximation as in the Opacity Project . Coupled channel calculations for $`\sigma _{\mathrm{PI}}`$ include both the background and the resonance structures (due to the doubly excited autoionizing states) in the cross sections. The recombination cross section, $`\sigma _{\mathrm{RC}}`$, is related to $`\sigma _{\mathrm{PI}}`$, through detailed balance (Milne relation) as

$$\sigma _{\mathrm{RC}}(ϵ)=\frac{\alpha ^2}{4}\frac{g_i}{g_j}\frac{(ϵ+I)^2}{ϵ}\sigma _{\mathrm{PI}}$$
(5)
in Rydberg units; $`\alpha `$ is the fine structure constant, $`ϵ`$ is the photoelectron energy, and $`I`$ is the ionization potential.

Resonant and non-resonant electronic recombination takes place into an infinite number of bound levels of the (e + ion) system. These are divided into two groups: (A) the low-$`n`$ ($`nn_o`$ 10) levels, considered via detailed close coupling calculations for photorecombination, with highly resolved delineation of autoionizing resonances, and (B) the high-$`n`$ ($`n_on\mathrm{}`$) recombining levels via DR, neglecting the background. In previous works (e.g. ) it has been shown that in the energy region corresponding to (B), below threshold for DR, the non-resonant contribution is negligible. The DR cross sections converge on to the electron impact excitation cross section at threshold ($`n\mathrm{}`$), as required by unitarity, i.e. conservation of photon and electron fluxes. This theoretical limit is an important check on the calculations, and enables a determination of field ionization of rydberg levels of resonances contributing to DR.

The ab initio method outlined above is a theoretically and computationally unified treatment based on the close coupling approximation. Recombination involves an infinite number of recombined bound states, and several infinite series of resonances. In principle, the unified method may be used for photoionization/photorecombination of arbitrarily high $`n,l,J`$ levels. However, in practice approximations may be made for sufficiently high quantum numbers. Background recombination is negligible, and DR dominates, usually for $`n_{max}`$ 10. Similarly, background (non-resonant) cross sections may be accurately obtained using hydrogenic approximation for $`n,l`$ levels with $`n>10`$. But there is nothing particular about $`n_{max}`$ = 10, and any larger or lower value may be used provided the approximations are verifiably valid, as has been shown in our previous works and is done in the present calculations. For example, $`n`$ = 18 – 22+ resonances in the present work are fully delineated using group (A) photorecombination calculations, and not DR. Thus the use of these approximations does not result in any significant error, or loss of generality, and does not detract from the main part of the calculations that are a unified representation of the non-resonant and resonant recombination (RR and DR), including any interference effects between the two. The present DR calculations use an extension of the precise theory by Bell and Seaton , based on multi-channel quantum defect theory, that is very accurate for high-$`n`$ (correspondence between photorecombination and DR is established in our previous work ). Finally, all close coupling scattering and photoionization calculations employ a “top-up” procedure for high partial waves, and approximations for high-$`n`$ resonances below Rydberg series limits as $`n\mathrm{}`$ (e.g. “Gailitis averaging”). Such procedures are routinely implemented in large-scale calculations in the Opacity Project ( and references therein) and the Iron Project R-matrix calculations () that the unified method for (e + ion) recombination is based upon.

## III COMPUTATIONS

The complete wavefunction expansion entails the 60 fine structure levels of Fe XVIII given in Table 1. These are obtained from an optimized configuration-interaction (CI) type calculation using the code SUPERSTRUCTURE . The configuration set is divided into the 5 spectroscopic configurations $`2s^22p^5,2s2p^6,2s^22p^43s,3p,3d`$ that dominate the 60 core level wavefunctions, and correlation configurations, $`2sp^53s,3p,3d`$, $`2p^63s,3p,3d`$. Calculated Fe XVIII eigenenergies are compared with experimental data from the National Institute for Standards and Technology (NIST). The accuracy of the eigenfunctions is also ascertained by comparing the Fe XVIII oscillator strengths for dipole transitions with available data from NIST in Table 2. Photoionization and recombination calculations both employ the Fe XVIII eigenfunctions with the same CI. We carry out two sets of calculations, (i) a 3-level calculations including only the $`n`$ = 2 levels, and (ii) the 60-level eigenfunction expansion including most of the $`n`$ = 3 complex. The inner $`2s`$-shell excitations are not considered owing to computational constraints and possibly weaker couplings (only the allowed $`2s3p`$ core excitations are likely to be of importance, and not the $`2s3s`$ or $`2s3d`$).

For the 3-level case, since we calculate both photoionization and photo-recombination cross sections, we include many $`LS\pi `$ symmetries to obtain 15 total $`J\pi `$ symmetries. Specifically these are $`J=07`$ for the even parity and $`J=06`$ for the odd parity, and the $`LS\pi `$ symmetries used are $`L=08`$ for both singlets and triplets in both parities. For the 60-level case, presently we only calculate photoionization cross section for the ground level and some selected excited levels. Therefore, we only include $`J=0`$ for the even ($`0^e`$) and $`J=1`$ for the odd parity ($`1^o`$) at this moment. All the $`LS\pi `$ symmetries that contribute to these two $`LS\pi `$’s are included. Of course, for obtaining photo-recombination cross sections we need to include an many $`J\pi `$ symmetries as in the 3-level case.

## IV RESULTS AND DISCUSSION

The following sections present a sample of the extensive results from the two sets of calculations for photoionization and (e + ion) recombination. The 3-level calculations are compared with earlier theoretical and experimental works. The present close coupling calculations for the 3-level and the 60-level cases are labeled 3CC and 60CC respectively.

### A Photoionization

Fig. 1a presents the BPRM photoionization cross section for the ground level $`2s^22p^6(^1S_0)`$ of Fe XVII from the 60CC calculation (solid line), showing the series of resonances converging on to the $`n`$ = 2 thresholds $`2s^22p^5(^2P_{1/2}^o)`$, $`2s2p^6(^2S_{1/2})`$, and the $`n`$ = 3 thresholds. For comparison, the non-resonant cross sections from a relativistic distorted wave (RDW) calculation (e.g. ) are also shown (dashed line). The resonance pattern, and the background cross sections, are essentially similar to the 3-level calculations in  in the relativistic random phase approximation (RRPA), with resonances included using multi-channel quantum defect theory (MQDT), and the $`LS`$ and Breit-Pauli $`JK`$-coupled R-matrix calculations also reported by Haque et al. . While there are no significant differences in magnitude or detail with the earlier calculations, it might be noted that the near-threshold region spanning the 3 levels of the $`n`$ = 2 complex, is rather simple in terms of structure and coupling effects. Fig. 1b shows an expanded view of this region with series of resonances $`2s^22p^5(^2P_{1/2}^o)n\mathrm{}`$ and the stronger series, connected to the ground level via a dipole core transition, $`2s2p^6(^2S_{1/2})n\mathrm{}`$. Fig. 2 shows a comparison between the 60CC and the 3CC calculations, with similar resolution, indicating that below the $`n`$ = 2 thresholds there is no significant difference between the two.

#### 1 Channel coupling effects

The situation is considerably more complicated above the $`n`$ = 2 complex. Although the ground level photoionization cross section of Fe XVII is not hugely affected by the $`n`$ = 3 complexes of resonances (Fig. 1a), the excited level cross sections are, as seen in Fig. 3a,b,c. This is of considerable importance in (e + ion) recombination work where photorecombination to group (A) levels is considered explicitly. The dense and detailed resonance structures converging on to the 57 $`n`$ = 3 levels, and in between, would enhance the effective photoionization and recombination cross sections and rates far above the background. It might be noted that the resonances in cross sections below the $`n`$ = 2, in the energy range covered by the 3-level calculations shown in Figs. 1 and 2, are much smaller than in the 60CC cross sections. This implies, in particular, that resonant recombination into the $`n`$ = 3 series of resonances converging on to a large number of excited states will be important for (e + ion) recombination, discussed in the next section.

As the above results show, the 3-level calculations in the present and earlier works (e.g. ) are inadequate for the entire energy range of interest in practical applications for photoionization and recombination. Also, a 3-level calculation gives little indication of the complexity of the cross sections, particularly for the excited states, since it covers only $`\mathrm{\Delta }n=0`$ core excitations and couplings that are responsible for resonances. The $`\mathrm{\Delta }n>0`$ couplings can be much stronger and give rise to more extensive resonances as in Fig. 3. It is clear that although Fe XVII is a highly charged ion the electron correlation effects are not weak in excited state photoionization, or in the near-threshold region. Finally, Fe XVII is a closed shell system where simpler approximations (e.g. Haque et al. ) can be readily applied without explicit consideration of detailed multiplet and fine structure that is more involved in open-shell systems. Thus photoionization of other highly charged ions may not be amenable to the approximations described in . In fact the atomic structure of open-shell Fe ions isoelectronic with the third row of the periodic table present considerable difficulties owing to strong coupling effects among up to several hundred fine structure levels. Such is the case in a number of scattering calculations carried out under the Iron Project for open-shell Fe ions where extensive resonance structures manifest themselves .

#### 2 Photoexcitation-of-core (PEC) resonances

Giant resonances manifest themselves at photon frequencies associated with strong dipole transitions in the core ion. These are a particularly important example of the coupling effects and are called the photoexcitation-of-core (PEC) resonances (e.g. ). The PEC resonances have the following properties,: (i) they are at the photon frequency of the dipole transition in the core, (ii) they are present in photoionization cross sections of the entire Rydberg series of bound levels of the (e + ion) system, (iii) their width and height are orders of magnitude larger than individual Rybderg resonances, and (iv) they are related to the inverse resonant recombination process DR (discussed in the next section). The PEC features are most evident in photoionization cross sections, as a function of photon energy, of several members of a Rydberg series of bound levels where, in analogy with the DR process, the outer electron is weakly bound and may be considered a ‘spectator’ interacting but weakly with the core excitation(s). Fig. 4 shows the large PEC feature in the 60CC photoionization cross sections of the series $`2p^5np^3P_0,n=310`$ levels of Fe XVII. The PEC resonances in Fig. 4 are associated with not just one dipole transition, but several transitions belonging to the transition arrays $`2p^52p^43s,3d`$ at about 63 Ry corresponding to all such levels included in the 60CC expansion of the core ion Fe XVIII (Eq. 1.1, Table 1). The PEC resonances rise order of magnitude above the background, and are much wider than all other resonances. The different threshold ionization energies in Fig. 4 approach $`0`$ as $`n`$ increases.

#### 3 Radiation damping of autoionizing resonances

Radiation damping of resonances has been addressed in many previous works. It was pointed out in Ref.  that it is likely to be of practical importance only for H-like and He-like core ions, when the core radiative transition rates are of the same order as the autoionization rates, typically $`10^{1314}sec^1`$, but not for other ionic systems. Fe XVII was explicitly mentioned in  as the next possible candidate (other than H-like and He-like ions) for an investigation of the radiation damping effects in toto. We discuss here the radiation damping involving the resonances associated with the $`n=2`$ levels in the context of low-energy recombination, and leave the discussion on the $`n=3`$ resonances in future when the much more extensive photo-recombination calculations with the 60CC target are presented. It might be noted however that although the radiative decay rates of the $`n`$ = 3 levels are higher, they also have additional autoionization modes of decay into excited $`n`$ = 2 levels.

The associated core transition rates in Fe XVIII for $`2s^12p^6(^2S_{1/2})`$ to $`2s^22p^5(^2P_{3/2}^o)`$ and $`2s^22p^5(^2P_{1/2}^o)`$ are $`9.13\times 10^{10}`$ and $`3.31\times 10^{10}`$ respectively. Fig. 5 shows an enhanced view of the radiatively damped (RD) and undamped (NRD) resonances in photoionization of Fe XVII up to resonances complexes with $`n`$ = 16. The PR calculations in the unified formulation generally employ cross sections up to $`n`$ = 10 only. No significant effect is discernible between Fig. 5a and 5b, and it is concluded that radiation damping of resonances in $`n`$-complexes up to $`n`$ = 10 (e.g. group (A)) is not likely to affect any practical applications of the computed photoionization and recombination cross sections. That is not to say that resonances with sufficiently high $`n`$ and $`\mathrm{}`$ will not be damped significantly (or completely); since the autoionization rates decrease as $`n^3`$, they must. However, such resonances are extremely narrow (not, for example, evident in Fig. 5), and do not affect effective photoionization or recombination cross sections. Therefore it is unlikely that radiation damping of group (A) resonances (an integral part of the unified (e + ion) calculations) in any other ionic system up to the iron-peak elements will be important, since the dipole transition probabilities of resonance transitions in all ions up to Fe ions are less than or equal to those in Fe XVIII, of the order of 10<sup>12</sup> sec<sup>-1</sup>, with the already noted exception of H-like and He-like ions .

### B Electron-ion recombination

Salient features of the (e + Fe XVIII) $``$ Fe XVII recombination are described within the unified formulation, and with reference to experimental data from the ion storage ring at Heidelberg, Germany .

#### 1 Comparison with experiment

Both the unified cross sections and experimental measurements naturally measure the combined non-resonant and resonant (RR and DR) contributions to (e + ion) recombination and should in principle be compared directly. Fig. 6 from Ref.  shows a comparison of the present unified cross sections as computed in detail, and averaged over a gaussian function for comparison with experiment, together with the experimental cross sections . The agreement is generally very good over the entire range, and in both the detail and the magnitude of resonances. This includes the dominant DR contribution below the $`{}_{}{}^{2}S_{1/2}^{}`$ threshold connected via a dipole transition; in this range the background non-resonant contribution (RR) is small. On the other hand the weaker series of resonances, $`{}_{}{}^{2}P_{1/2}^{}n\mathrm{}`$ lie in the near-threshold region dominated by the non-resonant contribution that rises steeply as $`E0`$, and is therefore not a major contributor to (e + ion) recombination rate (discussed below). Thus the unified theoretical (and experimental) results shown in Fig. 6 display the three related, but discernible, types of contributions to total (e + Fe XVIII) $``$ Fe XVII recombination cross section: the overlapping non-resonant (RR) contribution and the DR contribution from the $`{}_{}{}^{2}P_{1/2}^{o}n\mathrm{}`$ series, and the mainly DR contribution from the $`{}_{}{}^{2}S_{1/2}^{}n\mathrm{}`$ series.

#### 2 Resonance strengths and rate coefficients

The experimental data in Fig. 6 do not precisely delineate or identify the resonances, and the background contribution (RR-type) is not ascertained from the measurements, possibly owing to contribution from charge transfer . The blended resonance features are fitted to a beam response function to eliminate the background, and their energies are determined approximately according to $`n`$ and $`\mathrm{}`$, using non-relativistic $`\mathrm{}`$-dependent quantum defects in the Rydberg formula. However, $`\mathrm{}`$ is not a good quantum number and the number of $`\mathrm{}`$-resonances within an $`n`$-complex is not exactly known. The theoretical resonances on the other hand are uniquely identified with the intermediate coupling spectroscopic designation $`(S_iL_iJ_i)n\mathrm{}J\pi `$. Therefore a 1–1 correspondence between the experimental measurements and relativistic cross sections can not be established. Further, since the background contribution, although dominant at low energies as $`E0`$ (Fig. 6), is not considered, a direct comparison with the unified cross sections and the experimental data  is not possible.

Nonetheless, for $`n`$-complexes where the background contribution is small compared to the resonant part, we may compare the average ‘resonance strengths’ , although these are not exactly defined (see  for a definition of the resonance oscillator strength in terms of the integrated $`\frac{df}{dϵ}`$, the differential oscillator strength per unit energy). Figs. 7b,c,d show a detailed view at high resolution of the first three $`{}_{}{}^{2}S_{1/2}^{}n\mathrm{}`$ complexes, with $`n`$ = 6,7,8. In order to ensure complete resolution of resonances an energy mesh of up to 10<sup>-7</sup> eV was used before numerical integration. The integrated, and summed, resonance strengths for the $`n`$ = 6,7,8 complexes are 1201.2, 421.8, and 221.1 cm<sup>2</sup> eV, compared to experimental values of 1240.2, 412.0, and 253.9 . The present value for $`n`$ = 7 complex is higher than reported in Ref.  as it is recalculated with higher resolution. The value for the $`n`$ = 8 complex has been complemented by the contributions from $`J>`$ 7 symmetries, $`J\pi `$ = 8<sup>o</sup> and 9<sup>e</sup>; without which the value is 200.5 cm<sup>2</sup> eV. Although the theoretical resonance strengths were checked to have converged with respect to the energy mesh, they seem to be somewhat systematically lower than the reported experimental data (the theoretical MCBP and the MCDF values in  also showed the same trend).

Fig. 7a is rather different, in that it shows the much narrower $`{}_{}{}^{2}P_{1/2}^{o}n\mathrm{}`$ resonances with high $`n`$ = 18 – 22. While it is not apparent from the figure, nor from the experimental data in Fig. 6 (bottom most panel), there are 2,985 resonances found in the small 0-5 eV range just above ionization threshold of Fe XVII. The set of ($`E_o,\mathrm{\Gamma }_a,\mathrm{\Gamma }_r)`$ for all these resonances have been computed. The integrated resonance strengths for the $`n`$ = 18,19 and 20 complexes are 2290.9, 521.6, 290.1 cm<sup>2</sup> eV, compared to experimental values  of 2452.8, 605.7, and 336.5 respectively. Again, we find the integrated values to be up to about 10% lower than measured, although here the uncertainties are much greater than for the lower-$`n`$ $`{}_{}{}^{2}S_{1/2}^{}n\mathrm{}`$ resonances in Figs. 7b,c,d since they may not have been completely resolved and because the integration energy ranges are very closely spaced. Also, the precise range of angular momenta $`(\mathrm{},J)`$ in the ion storage rings is not known , and $`J`$ 7 may not be quite sufficient for perfect agreement. Detailed comparison of resonances beyond this level is neither feasible nor necessary. Although not all resonances are experimentally identified, the theoretical cross sections, resonance strengths, and rates (see below) can be compared, as shown, to within experimental and theoretical uncertainties.

Finally we evaluate quantities of practical interest, the maxwellian averaged unified (e + ion) recombination rate coefficients $`\alpha _R(T)`$ shown in Fig. 8. These are compared with the sum of the available experimental DR rate coefficients , and the background RR-type contribution extracted from the present theoretical cross sections. The agreement is of the order of 20%, the estimated uncertainty in both the experimental and theoretical datasets. However, we believe that the agreement could be slightly better since the theoretical results might be somewhat enhanced if the few $`J\pi `$ symmetries with $`J>7`$, for resonances with $`n10`$, are also included. These were omitted to reduce the complexity of calculations (except to gauge their effect on the completeness of the $`n`$ = 8 complex mentioned above), and since their resonant contribution is small.

#### 3 $`\mathrm{\Delta }n=1`$ resonances and high energy recombination

Thus far only the low energy $`\mathrm{\Delta }n=0`$ resonances due to the $`n=2`$ levels of Fe XVIII have been considered. However, as demonstrated in this work, the high energy recombination cross sections due to the $`n`$ = 3 levels are much larger and will be more important at high temperatures close to the temperature of maximum abundance of Fe XVII in collisional equilibrium, around $`T=45\times 10^6`$ K. This would especially be the case owing to the huge PEC resonances shown in Fig. 4 that in fact correspond to the peak values of DR due to resonances converging on to the series limits of strong dipole $`\mathrm{\Delta }n=1`$ transitions. Therefore, in addition to the DR bump corresponding to the $`n`$ = 2 resonances in the total $`\alpha _R(T)`$ shown in Fig. 8, we expect a much larger bump at higher temperatures from the $`n`$ = 3 resonances. However, the 60CC calculations are orders of magnitude more expensive in terms of computational and other resources and, although they are in progress, would require a considerable amount of time to be completed.

### C Bound states and transition probabilities of Fe XVII

In addition to photoionization and recombination the BPRM calculations also enable unprecedented quantities of accurate bound state and transition probability datasets in intermediate coupling (e.g. ). The accuracy of these results is comparable to the most elaborate configuration-interaction atomic structure calculations since the wavefunction in Eq. 2.1 entails a large configuration expansion, with each closed channel as a bound configuration of the (e + ion) system for $`E<0`$ and given symmetry $`J\pi `$. These calculations are in progress for all fine structure levels of Fe XVII up to $`n`$ = 10, and all associated E1 $`A`$\- and $`f`$-values. Table 3 provides a brief sample of the bound level energies computed for some of the levels of interest in this work. Together with the electron impact excitation collision strengths for Fe XVII in progress , these results should help complete the radiative and collisional data for Fe XVII needed for most plasma modeling applications.

## V CONCLUSION

The most extensive relativistic close coupling calculations for photoionization and recombination of an atomic species are reported for the astrophysically important ion Fe XVII. Based on this work we may note the following conclusions: (i) self-consistent datasets may be obtained for photoionization and recombination within the close coupling formulation, (ii) the coupling to the $`n`$ = 3 thresholds strongly manifests itself in excited state photoionization; resonances enhance the effective cross sections by orders of magnitude particularly below thresholds coupled via dipole photoexcitation-of-core levels, (iii) unified (e + ion) recombination cross sections are in good agreement with experimental data in terms of both detailed resonance strengths and rates; the resonances have been delineated at very high resolution with considerably more structure than experimentally observed (in principle all quantum mechanically allowed resonances in intermediate coupling may be obtained), (iv) it is necessary to consider the higher $`n`$ = 3 levels in high energy recombination that would dominate the DR part of the (e + Fe XVIII) $``$ Fe XVII recombination at temperatures close to maximum abundance of Fe XVII in coronal (collisional) equilibrium, and (v) it is necessary to account for not only relativistic fine structure but also the strong coupling among those levels in order to accurately reproduce the results for photoionization and recombination over the entire range of practical interest.

## VI ACKNOWLEDGEMENTS

We would like to thank Dr. Werner Eissner for several contributions, and Mr. Guoxin Chen for assistance with some calculations. This work was supported in part by the NSF and NASA. The computational work was carried out at the Ohio Supercomputer Center in Columbus, Ohio.


text:
# Bifurcation Diagram for Compartmentalized Granular Gases

## I Introduction

One of the key features of a granular gas is the tendency to spontaneously separate into dense and dilute regions goldhirsch93 ; mcnamara94 ; du95 ; jaeger96 ; kudrolli97 ; kadanoff99 . This clustering phenomenon manifests itself in a particularly clear manner in a box that is divided in a series of $`N`$ connected compartments, with a hole (at a certain height) in the wall between each two adjacent compartments. The system is vibro-fluidized by shaking the box vertically. With vigorous shaking the granular material is observed to be distributed uniformly over the compartments as in any ordinary molecular gas. Below a certain driving level however, the particles cluster in a small subset of the compartments, emptying all the others.

For $`N=2`$ the transition from the uniform to the clustered state is of second order, taking place through a pitchfork bifurcation eggers99 . For $`N=3`$ it was recently found that the transition is hysteretic. It is a first order phase transition, involving saddle-node bifurcations vdweele00 . This difference has been explained by a flux model. In the present paper we will use the same flux model to construct the bifurcation diagrams for arbitrary $`N`$.

The main ingredient of this model is a flux function $`F(n)`$, which gives the outflow from a compartment to one of its neighbors as a function of the fraction of particles ($`n`$) contained in the compartment eggers99 . The function $`F(n)`$ starts out from zero at $`n=0`$ and initially increases with $`n`$. At large values of $`n`$ it decreases again because the particles lose energy in the non-elastic collisions, which become more and more frequent with increasing particle density. So $`F(n)`$ is non-monotonic, and that is why the flux from a well-filled compartment can balance that from a nearly empty compartment.

Assuming that the granular gas in each compartment is in thermal equilibrium at any time (in the sense of the granular temperature mcnamara98 ) the following approximate form for $`F(n)`$ can be derived eggers99 :

$$F(n_k)=An_k^2e^{BN^2n_k^2},$$
(1)
which is a one-humped function, possessing the features discussed before (See Fig. 1). In the above equation $`n_k`$ is the fraction of particles in the $`k`$-th compartment, normalized to $`n_k=1`$. The factors $`A`$ and $`B`$ depend on the number of particles and their properties (such as the radius, and the restitution coefficient of the interparticle collisions), on the geometry of the system (such as the placement and form of the aperture between the compartments), and on the driving parameters (frequency and amplitude). The factor $`A`$ determines the absolute rate of the flux, and will be incorporated in the time scale, which thus becomes dimensionless. The clustering transition is governed only by $`B`$.

The time rate of change $`\dot{n}_k`$ of the particle fraction in the $`k`$-th compartment is given by the inflow from its two neighbours minus the outflow from the compartment itself,

$$\dot{n}_k=F(n_{k1})2F(n_k)+F(n_{k+1}),$$
(2)
with $`k=1,2,..,N`$. Here we have assumed that the interaction is restricted to neighboring compartments only.

For a cyclic arrangement the above equation is valid for all $`N`$ compartments (with $`k=N+1`$ equal to $`k=1`$). If we take non-cyclic boundary conditions, by obstructing the flux between two of the compartments, the equation has to be modified accordingly for these compartments.

The total number of particles in the system is conserved ($`_kn_k=1`$), so

$$\underset{k}{}\dot{n}_k=0.$$
(3)
Statistical fluctuations in the system would add a noise term to Eq. (2), but we will not consider such a term here. So the present analysis has to be interpreted as a mean field theory for the system.

Equation (2) can also be written in matrix-form, as $`\dot{𝐧}=𝐌𝐅`$, or more explicitly:

$`\dot{n}_k`$ $`={\displaystyle \underset{l}{}}M_{kl}F(n_l)`$ (4)
$`=\left(\begin{array}{ccccccc}2& 1& 0& 0& \mathrm{}& 0& 1\\ 1& 2& 1& 0& \mathrm{}& 0& 0\\ 0& 1& 2& 1& \mathrm{}& 0& 0\\ :& :& :& :& & :& :\\ :& :& :& :& & :& :\\ 1& 0& 0& 0& \mathrm{}& 1& 2\end{array}\right)\left(\begin{array}{c}F\left(n_1\right)\\ F\left(n_2\right)\\ F\left(n_3\right)\\ :\\ :\\ F\left(n_N\right)\end{array}\right)`$
The given matrix $`𝐌`$ corresponds to a cyclic arrangement of the $`N`$ compartments. A similar matrix can be written down for the case of a non-cyclic arrangement. We will come back to this later, when we will see that most of the results for the cyclic arrangement carry over to the non-cyclic case.

It is easily seen, from the fact that the elements of each row of $`𝐌`$ sum up to zero, that $`\mathrm{𝟏}=(1,1,\mathrm{},1)`$ is an eigenvector. The corresponding eigenvalue $`\lambda =0`$ physically reflects the fact that the compartments cannot all be filled (or emptied) simultaneously: $`_k\dot{n}_k=0`$ or $`_lM_{lk}=0`$. For future reference we note that all the other eigenvalues of $`𝐌`$ are negative (see Appendix).

The remainder of the paper is set up as follows. In Section II we show how to construct the bifurcation diagram, on the basis of Eq. (4), for an arbitrary number of compartments. In Section III we discuss the stability of the various branches in the diagram. Section IV discusses the physical consequences resulting from the diagram, in particular in the limit for $`N\mathrm{}`$. Finally, Section V contains concluding remarks. The paper is accompanied by a mathematical Appendix, in which some essential results concerning the stability analysis are derived.

## II Constructing the bifurcation diagram

To calculate the bifurcation diagram, we have to find the fixed points of Eq. (4) as a function of the parameter $`B`$, i.e., those points for which $`\dot{n}_k=𝐌𝐅=0`$. So $`𝐅`$ must be a multiple of the zero-eigenvalue vector $`\mathrm{𝟏}=(1,1,\mathrm{},1)`$. This tells us that, in a stationary situation, all components of the flux vector $`𝐅`$ are equal: there is a detailed balance between all pairs of neighboring compartments. This rules out, for instance, the possibility of stable standing-wave-like patterns with equal but non-zero net fluxes throughout the system. The fixed point condition now becomes

$$\{\begin{array}{cc}\hfill F(n_k)& =\text{constant}\hfill \\ \hfill n_k& =1\hfill \end{array}$$
(5)
Since $`F`$ is a one-humped function, $`F(n_k)=constant`$ has two solutions, which will be called $`n_{}`$ and $`n_+`$ (see Fig. 1). Every fixed point can be represented as a vector with elements $`n_{}`$ and $`n_+`$ (in any order, and summing up to 1) corresponding to a row of nearly empty and well-filled compartments. Let us call the number of well-filled compartments $`m`$. Apart from the ordering of the elements, every fixed point is then specified by only two numbers: $`n_+`$ and $`m`$.

Before actually calculating the bifurcation diagram, it is convenient to replace the fraction $`n`$ by the (also dimensionless) variable $`z=Nn\sqrt{B}`$, as then the flux (1) simplifies to $`F(z_k)z_k^2\mathrm{exp}(z_k^2)`$. The fixed point condition Eq. (5) then reads:

$$\{\begin{array}{cc}\hfill F(z_k)& =\text{constant}\hfill \\ \hfill z_k& =N\sqrt{B}\hfill \end{array}$$
(6)
So the $`B`$-dependence has been transferred from $`𝐅`$ to the particle conservation, and this enables us to determine the entire bifurcation diagram from one single graph. This is illustrated in Fig. 2 for the case of $`N`$=5 compartments.

First, the one-humped function $`F(z)`$ is inverted separately on both sides of the maximum, yielding the functions $`z_{}(F)`$ and $`z_+(F)`$. Then, we construct the sumfunctions:

$$S_m(F)=mz_+(F)+(Nm)z_{}(F)$$
(7)
Now, from Eq. (6), the fixed points are found by intersecting the horizontal line $`z=N\sqrt{B}`$ with the sumfunctions $`S_m(F)`$. In Fig. 2 this is done for $`B`$=1.08. Each intersection point yields a pair $`\{z_{},z_+\}`$, or equivalently $`\{n_{},n_+\}`$. Repeating the procedure for all $`B`$, we obtain the bifurcation diagram depicted in Fig. 3.

It contains several branches. First, a horizontal line (from the sumfunctions $`S_0`$ and $`S_5`$) corresponding to the equal distribution $`n_+=n_{}=0.2=1/N`$. Second, the branches corresponding to the $`m=1`$ clustered state (from $`S_1`$), which at $`B=1`$ goes over into the $`m=4`$ state (from $`S_4`$). And third, the branch of the $`m=2`$ clustered state (from $`S_2`$), which at $`B=1`$ becomes the $`m=3`$ state. The physical appearance of these solutions is sketched in the small diagrams. Note that only the $`m=0`$ branch (i.e., the uniform solution up to $`B=1`$) and the outer $`m=1`$ branch are stable. All the other branches are unstable, as will be discussed in the next section.

At $`B=1`$, where the branches intersect with the uniform distribution $`n_+=n_{}=1/N`$, we have a critical point. In the flux function one passes the maximum here. This means that $`n_+`$ and $`n_{}`$ are switched (relatively empty compartments become relatively filled, and vice versa), so $`m`$-branches change into $`(Nm)`$-branches. From a physical point of view, the most important thing that happens at the $`B=1`$ intersection point is the destabilization of the uniform distribution.

The saddle-node bifurcations of the $`m=1`$ and $`m=2`$ branches correspond to the minima of the sumfunctions $`S_1`$ and $`S_2`$ respectively, which in Fig 2 can be seen to occur at $`F0.014`$ for $`S_1`$ and $`F0.202`$ for $`S_2`$. In general, if a sumfunction $`S_m(F)`$ has a minimum for a certain $`B`$, the associated $`m`$ branch will have a bifurcation. So the bifurcation condition is that the derivative $`dS_m(F)/dF`$ equals zero, or equivalently:

$$\frac{\left({\displaystyle \frac{dz_{}}{dF}}\right)}{\left({\displaystyle \frac{dz_+}{dF}}\right)}=\frac{m}{Nm}$$
(8)
Not surprisingly, the quantity on the left hand side ($`dz_{}/dz_+\sigma `$) will play an important role in the stability analysis of the next section.

## III Stability of the branches

The stability of the branches (i.e., of the fixed points) is determined by the eigenvalues of the Jacobi matrix $`𝐉`$ corresponding to Eq. (4), with components:

$$J_{jk}=\frac{\dot{n}_j}{n_k}=\underset{l}{}M_{jl}F^{}(n_l)\frac{n_l}{n_k}=M_{jk}F^{}(n_k)$$
(9)
Here $`F^{}`$ denotes the derivative of $`F`$ with respect to $`n`$. Note that the Jacobi matrix can also be written as the product of $`𝐌`$ and the diagonal matrix $`𝐃=\mathrm{diag}(F^{}(n_1),\mathrm{},F^{}(n_N))`$, see also Eq. (18) in the Appendix. For a fixed point the only diagonal elements that occur are $`F^{}(n_+)`$ ($`m`$ times) and $`F^{}(n_{})`$ ($`Nm`$ times), in any order. The ratio between these two functions is precisely the quantity we encountered earlier in the bifurcation condition Eq. (8), namely $`\sigma `$:

$$\sigma =\frac{F^{}(n_+)}{F^{}(n_{})}=\frac{dn_{}}{dn_+}$$
(10)
The Jacobi matrix $`𝐉`$ has $`N`$ eigenvalues, one of which is always zero. The other $`N1`$ eigenvalues depend on $`m`$ and the value of $`\sigma `$.

For $`m=0`$ (the equipartitioned state) all non-trivial eigenvalues are negative, up to the point $`B=1`$. This can be seen either by direct numerical calculation, or analytically (see Appendix). At $`B=1`$, the $`m=0`$ state becomes the $`m=N`$ state. Here, the functions $`F^{}`$ in the Jacobi matrix (9) change sign, and so do all of its eigenvalues. So suddenly the uniform state has $`N1`$ positive eigenvalues, which implies a high degree of instability. Only in the limit $`B\mathrm{}`$ does the uniform state regain some of the lost terrain: the magnitude of all positive eigenvalues tends to zero here. Physically speaking, in this limit the vibro-fluidization is too weak to drive the particles out of the boxes anymore.

As for the other values of $`m`$, in Fig. 4 we have plotted the numerically evaluated eigenvalues (as functions of $`\sigma `$) for the system with $`N=5`$ compartments.

For $`m=1`$, we see that there are three eigenvalues that are always negative. The fourth non-trivial eigenvalue changes sign at $`\sigma =0.25`$. This corresponds to the saddle-node bifurcation of the $`m=1`$ branch in the bifurcation diagram (Fig. 3), and the bifurcation value of $`\sigma `$ is in agreement with Eq. (8). The region to the right of $`\sigma =0.25`$ (where all non-trivial eigenvalues are negative) belongs to the stable outer branch. The left part $`\sigma <0.25`$ belongs to the unstable inner branch, up to the point $`B=1`$ (at $`\sigma =1`$), where the $`m=1`$ branch goes over into the $`m=4`$ branch. That is, the state $`\{+\}`$ now switches to $`\{++++\}`$. At the same time all eigenvalues change sign, so suddenly we have $`3`$ positive eigenvalues, which is only one less than for the uniform $`m=5`$ state. (Indeed, the only stable manifold of the $`m=4`$ fixed point comes from the direction of the completely unstable $`m=5`$ state). The positive eigenvalues never cross zero anymore (there are no bifurcations beyond $`B=1`$) but, as before, in the limit $`B\mathrm{}`$ ($`\sigma 0`$) they go to zero.

For $`m=2`$ there are two possible configurations: $`\{++\}`$ and $`\{++\}`$. Due to the cyclic symmetry, all other combinations are equivalent to these two. The eigenvalues of the first configuration are given by the dotted lines, and those of the second by the solid lines. Although they are very similar (and are represented by exactly the same branch in the bifurcation diagram), it is clear that the second configuration is the more stable of the two. Apparently the two well-filled compartments prefer to keep a distance.

The saddle-node bifurcation of the $`m=2`$ branch takes place at $`\sigma =2/3`$ \[cf. Eq. (8)\], where the third non-trivial eigenvalue goes through zero. The fourth non-trivial eigenvalue always remains positive, indicating that the $`m=2`$ branch never becomes completely stable. (As a matter of fact, only the $`m=0`$ branch and part of the $`m=1`$ branch can be completely stable). Note that for $`\sigma 0`$ (large $`B`$) the positive eigenvalue tends to zero, so the degree of instability is quite weak there.

At $`B=1`$ the $`m=2`$ branch becomes the $`m=3`$ branch, with the two configurations $`\{+++\}`$ and $`\{+++\}`$, and with all eigenvalues switching sign. As we see, the more dispersed configuration is again the less unstable one. Also the phenomenon of all positive eigenvalues going to zero as $`\sigma `$ approaches zero (the weak driving limit $`B\mathrm{}`$) is again apparent.

In the present example for $`N=5`$, and in fact for all odd values of $`N`$, the branches in the bifurcation diagram are all born by means of a saddle-node bifurcation. But for even values of $`N`$ this is different: in that case there is one branch-pair that springs from the uniform distribution, at $`B=1`$, by a pitchfork bifurcation. This is illustrated in Fig. 5 for $`N=6`$. Here one sees all the branches that were present already for $`N=5`$, only slightly shifted towards the left, plus an additional pair of branches ($`m=3`$) bifurcating in the forward direction from $`B=1`$.

The special status of the branch $`m=N/2`$ is also evident from Eq. (8), which tells us that the bifurcation condition for this branch is $`\sigma =1`$. This condition is fulfilled only by $`n_+=n_{}=1/\sqrt{B}=1/N`$. So, unlike all other branches, this one originates at $`B=1`$ from the (until then stable) uniform state. Related to this, the branch is the only one that is symmetric for interchanging $`n_+`$ and $`n_{}`$.

## IV Physical aspects

The bifurcation analysis from the previous section can also be understood from a more physical point of view. To this end, let us first have a closer look at a 2-box system. In the equilibrium situation the net flux between the two boxes is zero, with one filled ($`n_+`$) and one nearly empty ($`n_{}`$) box. Suppose the level of the empty box is raised by an amount $`\delta n`$. The level of the filled box then decreases by an equal amount and the net flux $`\varphi _+`$ from the empty to the filled box becomes (see also Fig. 1):

$`\varphi _+`$ $`=F(n_{}+\delta n)F(n_+\delta n)`$ (11)
$`=\left({\displaystyle \frac{dF}{dn_{}}}+{\displaystyle \frac{dF}{dn_+}}\right)\delta n=\left(1+\sigma \right){\displaystyle \frac{dF}{dn_{}}}\delta n`$
where we have used that $`\sigma =dn_{}/dn_+`$ and neglected the higher order terms in the Taylor expansion. There are two different regimes. If $`\sigma >1`$, the net flux is positive (as $`F^{}(n_{})`$ is always positive), so particles are flowing from $`n_{}`$ to $`n_+`$, restoring the equilibrium position. This is actually the situation along the entire $`m=1`$ branch, for all $`1<B<\mathrm{}`$. For $`\sigma <1`$ (a situation which does not occur for our choice of $`F`$), the net flux would be negative, raising the level of the emptier box even further, away from the equilibrium position. In the borderline case, $`\sigma =1`$ (at $`B=1`$), the system is indifferent to infinitesimal changes.

This argument is readily generalized to the $`N`$-compartment system, for an equilibrium with $`m`$ filled boxes. Now we raise the level of all $`Nm`$ nearly empty boxes simultaneously by $`\delta n`$. This is done by lowering all levels in the $`m`$ filled boxes by an equal amount, which by particle conservation must be equal to $`\delta n(Nm)/m`$. The equivalent of Eq. (11) for the flux between any of the empty boxes to a neighbouring filled box then reads:

$`\varphi _+`$ $`=F(n_{}+\delta n)F(n_+{\displaystyle \frac{Nm}{m}}\delta n)`$ (12)
$`=\left({\displaystyle \frac{dF}{dn_{}}}+{\displaystyle \frac{Nm}{m}}{\displaystyle \frac{dF}{dn_+}}\right)\delta n`$
$`=\left({\displaystyle \frac{m}{Nm}}+\sigma \right){\displaystyle \frac{dF}{dn_{}}}\delta n`$
From this expression it follows that the transition between a (relatively) stable ($`\sigma >m/(Nm)`$) and a (relatively) unstable ($`\sigma <m/(Nm)`$) configuration is marked by the bifurcation condition Eq. (8). So, by straightforward physical reasoning we have reproduced the exact result obtained earlier from an eigenvalue analysis.

The pitchfork bifurcation discussed at the end of Section III is especially important for $`N=2`$. In this case it is the only non-uniform branch. To be specific, it is a stable $`m=1`$ branch. This $`N=2`$ case eggers99  is the only one without any saddle-node bifurcations, and consequently it is the only case where the change from the uniform to the clustered situation takes place via a second order phase transition without any hysteresis. For all $`N>2`$ the transition is of first order vdweele00 , and shows a hysteretic effect that becomes more pronounced for growing $`N`$.

In the limit $`N\mathrm{}`$ the hysteresis is maximal: the first saddle-node bifurcation takes place immediately after $`B=0`$, and this means that there exists a stable $`m=1`$ solution over the entire range $`B>0`$. So, if one starts out from this solution (at a certain value of $`B`$) and then gradually turns down $`B`$, one will never witness the transition to the uniform distribution. Vice versa, also the transition from the uniform solution to the $`m=1`$ state will not occur in practice, even though the uniform distribution becomes unstable at $`B=1`$. If one starts out from the uniform solution (at a certain value of $`B`$ below 1) and increases $`B`$, one will witness the transition to a clustered state, but in practice this will always be one with a number of peaks. That is, the system gets stuck in a transient state with $`m>1`$, even though such a state is not stable (it has one or more positive eigenvalues).

The fact is that its lifetime may be exceedingly large, since the flux in the neighborhood of a peak and its adjacent boxes (which are practically empty) is very small. Furthermore, the communication between the peaks is so poor that usually (even for moderate values of $`N`$) the dynamics comes to a standstill in a state with peaks of unequal height.

Another point we would like to address is that practically the transition to a clustered state will take place already before $`B=1`$, because the solution is kicked out of its basin of attraction by the statistical fluctuations in the system vdweele00 . An example is shown in Fig. 6. Here we see a snapshot for the cyclic system with $`N=80`$ compartments, which were originally filled almost uniformly, at $`B=0.90`$. The small random fluctuations in the initial condition are sufficient to break away from the (still stable) uniform distribution, and one witnesses the formation of a number of isolated clusters. In the further evolution these clusters deplete the neighbouring compartments and indeed the whole intermediate regions. But the peaks themselves, once they are well-developed, do not easily break down anymore.

## V Concluding remarks

In this paper we have constructed the bifurcation diagram for a vibro-fluidized granular gas in $`N`$ connected compartments. Let us now comment upon the result.

Starting out from $`B=0`$, i.e, vigorous shaking, the equi-partitioned state is for some time the only (and stable) fixed point of the system. For increasing $`B`$ we first come upon the $`m=1`$ bifurcation, where the single-cluster state is born. For all $`N>2`$ this happens by means of a saddle-node bifurcation, creating one completely stable state and one unstable state (with $`1`$ positive eigenvalue). The one with the largest difference between $`n_+`$ and $`n_{}`$ is the stabler one of the two states. Strictly speaking, there are $`N`$ equivalent single-cluster states, since the cluster can be in any of the $`N`$ compartments.

For further growing $`B`$ we come across the $`m=2`$ bifurcation, where two unstable 2-peaked states are created. The state with the largest difference between $`n_+`$ and $`n_{}`$ has $`1`$ positive eigenvalue, and the other one $`2`$. The two peaks can be distributed in $`\left(\genfrac{}{}{0pt}{}{N}{2}\right)`$ ways over the $`N`$ compartments, but as we have seen they are not all equivalent. When the peaks are situated next to each other we have a more unstable situation (the positive eigenvalues are larger in magnitude) than when the peaks are further apart. This is generally true for $`m`$-peaked solutions: of the $`\left(\genfrac{}{}{0pt}{}{N}{m}\right)`$ ways in which $`m`$ peaks can be distributed, the ones in which the peaks are next to each other are the least favorable of all.

For increasing $`B`$ we encounter more and more bifurcations, where unstable $`m`$-clustered states come into existence (each with $`1`$ more positive eigenvalue than the previous one), and for large $`N`$ the bifurcation diagram is covered by a dense web of branches. In Fig. 7 this is shown for N=80. The last saddle-node bifurcation takes place shortly before $`B=1`$ and, for this even value of $`N`$, is followed by a final pitchfork bifurcation (creating the $`m=N/2`$ branch) at $`B=1`$.

The uniform solution (or $`m=0`$ state) is stable until $`B=1`$, with $`N1`$ negative eigenvalues and $`1`$ zero. For $`B>1`$ all its negative eigenvalues become positive, making it suddenly the most unstable state of all. Also, it now formally becomes the $`m=N`$ state. Moving away from this uniform solution one encounters first the $`m=N1`$ branch with $`N2`$ positive eigenvalues, then the $`m=N2`$ branch with $`N3`$ positive eigenvalues, etc. Finally, one arrives at the outermost $`m=1`$ branch, which has no positive eigenvalues. This is the only solution that is completely stable for $`B>1`$. But as we have seen in the previous section, on its way from the uniform distribution to the single peaked state, the system can easily get stuck in one of the transient states (especially for large $`N`$) even though these are not strictly stable.

Throughout the paper, we have concentrated on the case where the $`N`$ compartments are arranged in a cyclic manner. But in doing so, we have in fact also solved the non-cyclic case. Here we close the hole in the wall between the $`1`$st and $`N`$th compartment, and consequently the flux between them is zero. The matrix $`𝐌`$ then takes the following form \[differing from the cyclic one only in the first and last row, cf. Eq. (4)\]:

$$𝐌^{(nc)}=\left(\begin{array}{ccccccc}1& 1& 0& 0& \mathrm{}& 0& 0\\ 1& 2& 1& 0& \mathrm{}& 0& 0\\ 0& 1& 2& 1& \mathrm{}& 0& 0\\ :& :& :& :& & :& :\\ :& :& :& :& & :& :\\ 0& 0& 0& 0& \mathrm{}& 2& 1\\ 0& 0& 0& 0& \mathrm{}& 1& 1\end{array}\right)$$
(13)
The eigenvalue problem for this matrix is treated in the Appendix. One eigenvalue is identically zero, and the other $`N1`$ eigenvalues are negative, just like for the cyclic system. This leads to a bifurcation diagram that is indistinguishable from that of the cyclic case. Even the stability along the branches is the same; only the magnitude (not the sign) of the eigenvalues of the Jacobi-matrix $`𝐉`$ is slightly different for the two cases.

Finally, it should be emphasized that the results of the present paper do not depend on the precise form of the flux function. We have concentrated on the form given by Eq. (1), but virtually everything remains true for other choices of this function, as long as it is a non-negative, one-humped function, starting out from zero at $`n=0`$ (no flux if there are no particles) and going down to zero again for very many particles (no flux also in this limit, since - due to the inelastic collisions - the particles form an inactive cluster, unable to reach the hole in the wall anymore). Any function with these properties will produce a bifurcation diagram similar to that of Eq. (1).

In the likely case that the range of $`\sigma =dn_{}/dn_+`$ is the same, extending from $`1`$ (this value is attained in the maximum) to zero (in the outer regions of the flux function, for $`n_{}\sqrt{B}0`$, $`n_+\sqrt{B}\mathrm{}`$), the bifurcation diagram will have the same number of saddle-node bifurcations and the same number of branches. The only things that change are the exact position of the bifurcation points, and the magnitude of the eigenvalues along the branches.

Slight differences in the diagram would occur if the slope of $`F`$ on the $`n_+`$ side was to become steeper than on the $`n_{}`$ side. In that case, the bifurcation condition Eq. (8) would also have solutions for $`m>N/2`$, thus allowing saddle-node bifurcations for branches with $`m>N/2`$. These branches, however, would certainly be quite unstable.

## APPENDIX: On the eigenvalues of $`𝐌`$ and $`𝐉`$

In this appendix we present the analytical eigenvalues of the flux matrix $`𝐌`$ \[introduced in Eqs. (4) and (13)\] and discuss the eigenvalue problem for the Jacobian matrix $`𝐉`$ \[see Eq. (9)\], thereby determining the stability of the branches in the bifurcation diagram.

First, we briefly treat the eigenvalues of $`𝐌`$. After that, we turn to $`𝐉`$. In Subsection 2 we discuss its zero-eigenvalues: one eigenvalue is identically zero and, by pinpointing the zero-crossing of a second eigenvalue, we reproduce the bifurcation condition Eq. (8). In Subsection 3 we determine the number of negative eigenvalues of $`𝐉`$ in the low-driving limit $`\sigma 0`$. Likewise, in Subsection 4 we determine the number of positive eigenvalues in the (mathematical) limit $`\sigma \mathrm{}`$. Combining these two results, in Subsection 5, we finally find the number of positive eigenvalues of $`𝐉`$ for general values of $`\sigma `$, and this gives the stability of the branches over the entire bifurcation diagram.

### .1 Eigenvalues of matrix $`𝐌`$

The matrix $`𝐌`$ in Eq. (4) is closely related to the $`N\times N`$ tridiagonal matrix $`\mathrm{tridiag}(1,2,1)`$ associated with the second difference operator known from numerical schemes for solving second order pde’s. Its eigenvalue problem can be solved exactly guardiola82 , and the same is true for $`𝐌`$. The eigenvalues of $`𝐌`$ are given by:

$$\lambda _k=4\mathrm{sin}^2\left(\frac{k\pi }{N}\right)$$
(14)
where $`k`$ runs from $`0`$ to $`N/2`$ for $`N`$ even, and from $`0`$ to $`(N1)/2`$ for $`N`$ odd. The corresponding eigenvectors are:

$$a_i(k)=C_1\mathrm{cos}\left(\frac{(2i+1)k\pi }{N}\right)+C_2\mathrm{sin}\left(\frac{(2i+1)k\pi }{N}\right)$$
(15)
with $`i=1,\mathrm{},N`$ and arbitrary coefficients $`C_1`$ and $`C_2`$.

As we see, the first eigenvalue ($`k=0`$) is zero and the corresponding eigenvector is $`\mathrm{𝟏}=(1,1,\mathrm{},1)`$. Physically, this eigenvector represents simultaneous filling of all $`N`$ compartments, and the eigenvalue $`0`$ expresses the fact that this is prohibited (because the number of particles in our system is conserved).

All non-zero eigenvalues are negative and (except the one for $`k=N/2`$ in the case of even $`N`$) doubly degenerate. This means that the corresponding eigenvectors span a two-dimensional subspace, reflected by the two terms $`C_1`$ and $`C_2`$ in Eq. (15). Since $`𝐌`$ is symmetric, and therefore normal, linear subspaces corresponding to different eigenvalues are orthogonal. Especially, the eigenvectors of all non-zero eigenvalues span a $`N1`$ dimensional subspace perpendicular to $`\mathrm{𝟏}=(1,1,\mathrm{},1)`$.

The matrix $`𝐌^{(nc)}`$ for the non-cyclic case, given by Eq. (13), has a different set of eigenvalues:

$$\lambda _k^{(nc)}=4\mathrm{sin}^2\left(\frac{k\pi }{2N}\right)$$
(16)
Here $`k`$ runs from $`0`$ to $`N1`$. The corresponding eigenvectors are:

$$a_i^{(nc)}(k)=\mathrm{cos}\left(\frac{(2i+1)k\pi }{2N}\right)$$
(17)
Just like in the cyclic case, the first eigenvalue equals zero, and all the others are negative. However, they are non-degenerate and the corresponding eigenspaces are one-dimensional.

### .2 Zero-eigenvalues of matrix $`𝐉`$

Now we turn to the Jacobian matrices. We consider the cyclic version $`𝐉`$, with components as given in Eq. (9), but the results are also valid for the non-cyclic version. This matrix can be written as the product of $`𝐌`$ and a diagonal matrix $`𝐃=\mathrm{diag}(F^{}(n_1),F^{}(n_2),\mathrm{},F^{}(n_N))`$:

$$𝐉=𝐌𝐃=\left(\begin{array}{cccccc}2F^{}\left(n_1\right)& F^{}\left(n_2\right)& 0& \mathrm{}& 0& F^{}\left(n_N\right)\\ F^{}\left(n_1\right)& 2F^{}\left(n_2\right)& F^{}\left(n_3\right)& \mathrm{}& 0& 0\\ 0& F^{}\left(n_2\right)& 2F^{}\left(n_3\right)& \mathrm{}& 0& 0\\ :& :& :& & :& :\\ :& :& :& & :& :\\ 0& 0& 0& \mathrm{}& 2F^{}\left(n_{N1}\right)& F^{}\left(n_N\right)\\ F^{}\left(n_1\right)& 0& 0& \mathrm{}& F^{}\left(n_{N1}\right)& 2F^{}\left(n_N\right)\end{array}\right)$$
(18)
In the context of the bifurcation diagram, the main thing one wants to know is the number of positive eigenvalues of $`𝐉`$ for each branch. This is what we are going to determine now.

First we note that the eigenvalues of $`𝐉`$ are real, even though the matrix is not symmetric. This is a consequence of the following similarity relationship between $`𝐉`$ and $`𝐉^{}`$:

$$𝐉^{}=(𝐌𝐃)^{}=𝐃𝐌=𝐃(𝐌𝐃)𝐃^1=𝐃𝐉𝐃^1$$
(19)
This implies that $`𝐉`$ and $`𝐉^{}`$ have the same eigenvalues, and hence they must be real. Because $`𝐌`$ is singular, $`𝐉`$ must be too (it has a zero eigenvalue) and so its determinant $`det(𝐉)`$ is zero. More explicitly:

$$det(𝐉)=det(𝐌)det(𝐃)=\left(\underset{k}{}F^{}(n_k)\right)det(𝐌)=0$$
(20)
where, for a fixed point with $`m`$ filled compartments, the product term equals $`[F^{}(n_+)]^m[F^{}(n_{})]^{(Nm)}`$.

For the other eigenvalues we have to look at the characteristic equation $`det(𝐉\lambda 𝐈)=0`$. This is a polynomial expression in $`\lambda `$, of which the constant term is zero since it is equal to $`det(𝐉)`$. The coefficient $`L`$ of the linear term is:

$$L=\underset{k}{}det(𝐉^{(k,k)})=\underset{k}{}\left(\underset{lk}{}\lambda _l\right)$$
(21)
where the matrix $`𝐉^{(k,k)}`$ is the $`(N1)\times (N1)`$ matrix obtained from $`𝐉`$ by deleting its $`k`$-th row and its $`k`$-th column. In the right-hand side of this equation, the only product that survives is the one that does not contain the trivial (zero) eigenvalue. So:

$$L=\underset{\text{all non-trivial}\lambda _l}{}\lambda _l$$
(22)
Alternatively, the determinant of $`𝐉^{(k,k)}`$ in Eq. (21) can be written in terms of $`det(𝐌^{(k,k)})`$, by deleting the $`k`$-th factor from the product in Eq. (20):

$$L=\underset{k}{}\left(\underset{lk}{}F^{}(n_l)\right)det(𝐌^{(k,k)})$$
(23)
It can be shown that for all $`k`$ the determinant $`det(𝐌^{(k,k)})`$ is a constant, $`C`$, which equals $`(N1)(1)^{N1}`$ in the cyclic, and $`(1)^{N1}`$ in the non-cyclic case. Thus, Eq. (23) reduces to:

$$L=C\underset{k}{}\left(\underset{lk}{}F^{}(n_l)\right)$$
(24)
For a fixed point with $`m`$ filled compartments, we can write (using that in the above summation each of the products misses either an $`F^{}(n_+)`$ or an $`F^{}(n_{})`$):

$`L`$ $`=C[F^{}(n_+)]^{(m1)}[F^{}(n_{})]^{(Nm1)}\times `$ (25)
$`\left((Nm)F^{}(n_+)mF^{}(n_{})\right)`$
$`=C[F^{}(n_+)]^{(m1)}[F^{}(n_{})]^{(Nm)}\left((Nm)\sigma m\right)`$
From this equation we conclude that $`L`$ becomes zero at $`\sigma =m/(Nm)`$. This is exactly the bifurcation condition already given in the main text \[Eq. (8)\]. Also, with Eq. (22), we see that an eigenvalue crosses zero at this value of $`\sigma `$.

It can be shown, by a similar analysis, that the coefficient of the quadratic term is not equal to zero at $`\sigma =m/(Nm)`$, so not more than one of the eigenvalues changes sign at the bifurcation.

### .3 Number of negative eigenvalues of $`𝐉`$ for $`\sigma 0`$

We now come to the next step in determining the number of positive eigenvalues. We again use the definition of $`\sigma `$ to write: $`𝐉=F^{}(n_{})𝐌\stackrel{~}{𝐃}`$, where $`\stackrel{~}{𝐃}=\mathrm{diag}(1,\mathrm{},1,\sigma ,\mathrm{},\sigma )`$. The factors $`1`$ correspond to the $`Nm`$ nearly empty boxes and the factors $`\sigma `$ to the $`m`$ filled boxes. The precise ordering of the factors is not essential for the following argument, so we may choose the above order for notational convenience.

The factor $`F^{}(n_{})`$ is always positive, so we only have to deal with $`𝐌\stackrel{~}{𝐃}`$. Note that only $`\stackrel{~}{𝐃}`$ depends on $`\sigma `$ and that in the limit $`\sigma 0`$ this matrix becomes<sup>1</sup><sup>1</sup>1Since both the mapping $`\sigma 𝐉_\sigma `$ as $`𝐉_\sigma det(𝐉_\sigma \lambda 𝐈)`$ are $`\text{C}^{\mathrm{}}`$ mappings (from $`^{N\times N}`$ and $`^{N\times N}\text{P}(\lambda )`$ respectively, where $`\text{P}(\lambda )`$ is the space of polynomials of order $`N`$), it is allowed to take the limits $`\sigma 0`$ and $`\sigma \mathrm{}`$, even though the latter does not occur in practice.:

$$\underset{\sigma 0}{lim}\stackrel{~}{𝐃}=\mathrm{diag}(1,\mathrm{},1,0,\mathrm{},0)𝐏$$
(26)
$`𝐏`$ is a projection matrix which projects $`^N`$ to the subspace spanned by the first $`Nm`$ unit vectors. It is obviously non-singular, symmetric, and applying it twice gives the same result as once: $`𝐏^2=𝐏`$.

Instead of taking the matrix $`𝐉_0=𝐌𝐏`$ as input for solving our eigenvalue problem (in the limit $`\sigma 0`$), we will rather look at the matrix $`𝐏𝐌𝐏`$ which is symmetric and has the same eigenvalues as $`𝐉_0`$.

For proof of the last statement, let $`\mu `$ be a (non-zero) eigenvalue of $`𝐉_0`$: $`𝐉_0𝐱=\mu 𝐱`$. Then: $`(𝐏𝐌𝐏)(𝐏𝐱)=𝐏(𝐌𝐏𝐱)=\mu (𝐏𝐱)`$. Note that $`𝐏𝐱0`$, because otherwise also $`𝐉_0𝐱=𝐌𝐏𝐱`$ would be zero, contradicting the assumption that $`\mu `$ is non-zero. This completes the proof.

The matrix $`𝐌`$ is negative semi-definite. This means that $`𝐌`$ has only negative or zero eigenvalues or, equivalently, the inner product $`𝐱,𝐌𝐱0`$ for all $`𝐱`$. This means that also $`𝐏𝐌𝐏`$ is negative semi-definite, because:

$$𝐱,𝐏𝐌𝐏𝐱=𝐏𝐱,𝐌(𝐏𝐱)=𝐲,𝐌𝐲0$$
(27)
In conclusion, $`𝐉_0`$ has negative and zero eigenvalues only.

The remaining task is to identify the number of negative eigenvalues, or otherwise stated, the rank of the matrix $`𝐉_0`$. The statement which we shall prove is that $`\mathrm{rank}(𝐉_\mathrm{𝟎})=\mathrm{rank}(𝐏)=Nm`$.

Proof: Note that the image $`\mathrm{Im}(𝐏)`$ of $`𝐏`$ is spanned by the first $`m`$ unit vectors of $`^N`$. Its kernel $`\mathrm{Ker}(𝐏)`$ is spanned by the remaining $`Nm`$ unit vectors. Since the kernel of $`𝐌`$ is spanned by the vector $`\mathrm{𝟏}`$, the following identities hold:

$$\mathrm{Ker}(𝐏)\mathrm{Ker}(𝐌)=\mathrm{𝟎}$$
(28a)
$$\mathrm{Im}(𝐏)\mathrm{Ker}(𝐌)=\mathrm{𝟎}$$
(28b)
Now, for all $`𝐱\mathrm{Ker}(𝐏)`$ it holds that $`𝐉_0𝐱=𝐌(𝐏𝐱)=\mathrm{𝟎}`$, so $`\mathrm{Ker}(𝐏)\mathrm{Ker}(𝐉_0)`$. On the other hand, for all $`𝐲\mathrm{Ker}(𝐏)`$ one has $`𝐏𝐲𝐳0`$, with $`𝐳\mathrm{Im}(𝐏)`$, and therefore $`𝐉_0𝐲=𝐌𝐳0`$ because of Eq. (28b). This means that $`𝐲\mathrm{Ker}(𝐉_0)`$, and thus $`\mathrm{Ker}(𝐏)\mathrm{Ker}(𝐉_0)`$. Together these two results prove that $`\mathrm{Ker}(𝐏)=\mathrm{Ker}(𝐉_0)`$, so obviously the rank of the two matrices must be equal. Since $`\mathrm{rank}(𝐏)=Nm`$, this is also the rank of $`𝐉_0`$, which completes the proof.

In short, we have shown that in the limit $`\sigma 0`$, the Jacobi-matrix $`𝐉`$ has $`Nm`$ negative eigenvalues.

### .4 Number of positive eigenvalues for $`\sigma \mathrm{}`$

We now turn to the limit $`\sigma \mathrm{}`$. In this limit we rewrite $`𝐉`$ as follows: $`𝐉=F^{}(n_{})\sigma 𝐌\overline{𝐃}`$. Here $`\overline{𝐃}=\mathrm{diag}(\sigma ^1,\mathrm{},\sigma ^1,1,\mathrm{},1)`$, which in the limit $`\sigma \mathrm{}`$ becomes:

$$\underset{\sigma \mathrm{}}{lim}\overline{𝐃}=\mathrm{diag}(0,\mathrm{},0,1,\mathrm{},1)𝐐$$
(29)
Again, $`𝐐`$ is a projection matrix, which now projects $`^N`$ to the subspace spanned by the last $`m`$ unit vectors, so $`𝐐`$ is complementary to $`𝐏`$. Following the same line of reasoning, but keeping in mind that now the constant factor in front of $`𝐉_{\mathrm{}}`$ is negative, we find that in the limit $`\sigma \mathrm{}`$, the matrix $`𝐉`$ has $`m`$ positive eigenvalues.

### .5 Number of positive eigenvalues of $`𝐉`$ for general $`\sigma `$

We are now ready to draw the conclusion. Just below $`\sigma =0`$ the matrix $`𝐉`$ must, by continuity, have at least $`Nm`$ negative eigenvalues. If we now move from $`0`$ towards $`\mathrm{}`$, beyond a certain point there must be at least $`m`$ positive eigenvalues (or equivalently, at most $`Nm1`$ negative eigenvalues). We already know \[cf. Eq. (25)\] that along the way exactly one eigenvalue changes sign, at $`\sigma =m/(Nm)`$. Taken together, this means that $`𝐉`$ has $`m`$ positive and $`Nm1`$ negative eigenvalues for $`\sigma <m/(Nm)`$, and $`m1`$ positive and $`Nm`$ negative eigenvalues for $`\sigma >m/(Nm)`$.

This completes the determination of the number of positive eigenvalues for the various branches in the bifurcation diagram.


text:
# Nonequilibrium transport through a quantum dot weakly coupled to Luttinger liquids

## Figure Captions

* Current, scaled by its maximum at $`\frac{\epsilon }{2T}=5`$, vs bias voltage for $`g=0.9,0.75,\frac{1}{2},\frac{1}{3},\frac{1}{10}`$. $`\epsilon `$ is the energy level of the dot, $`g`$ is the interaction parameters of the leads. Curves for $`g=\frac{1}{2},\frac{1}{3},\frac{1}{10}`$ intersect approximately at a point $`\left(\frac{eV}{2T}=10,\frac{J}{J_{\mathrm{max}}}=\frac{1}{2}\right)`$. The case with FL leads $`\left(g=1\right)`$ is also shown for comparison..
* Differential conductance scaled by its maximum at $`\frac{eV}{2T}=10`$ as a function of $`\frac{\epsilon }{2T}`$ for $`g=1,\frac{1}{2},\frac{1}{10}`$. The location of the maximum is shifted away from $`\epsilon =eV`$ for $`g=1`$ to $`\epsilon =\frac{eV}{2}`$ for $`g0`$.
* Current vs $`\frac{\epsilon }{2T}`$ for $`g=1,\frac{1}{2},\frac{1}{10}.`$ The current is scaled by its maximum at bias voltage $`\frac{eV}{2T}=10,`$ For $`g<1`$(Luttinger liquid), a peak occurs at $`\epsilon =\frac{eV}{2}`$, and the half width of the peak becomes narrower with the decrease of $`g`$. The curve with $`g=1`$ (Fermi liquid) is also shown for comparison.