text:
# Lattice Models of Ionic Systems with Charge Asymmetry

## I Introduction

The nature of critical phenomena in ionic systems has been a subject of numerous theoretical studies in recent years.fisher93 ; stell ; levin96 ; weingartner01  Due to the long-range nature of Coulombic interactions, construction of a full renormalization group treatment, which was so successful in describing critical behavior of non-ionic fluids, meets both conceptual and technical difficulties.fisher93  However, in recent years some progress has been achieved in obtaining physically reasonable, well-based mean-field theories for ionic systems.fisher93 ; stell ; levin96 ; weingartner01 ; levin02  These theoretical studies have been supported and, sometimes in substantial part, initiated by intensive Monte Carlo simulations panagiotopoulos ; orkoulas99 ; camp99 ; lujiten01 ; pana02 ; yan02  of charged systems.

To investigate the thermodynamics of ionic fluids, two main mean-field approaches have been developed. The first onefisher93 ; levin96  extends the pioneering work of Debye and Hückeldebye  (DH) on dilute solutions of strong electrolytes, while the second approachstell  relies on integral equations for correlation functions. Analysis of thermodynamic energy bounds and comparison with the best Monte Carlo estimates for the critical parameters suggests that the DH-based theory gives a better description of the thermodynamics of electrolytes, at least in the critical region.levin96 ; zuckerman97

The simplest model of ionic fluids, the so-called restrictive primitive model (RPM), considers a system of spherical equisized charged particles, half of them carrying a charge $`q`$ and the other half with charge $`q`$. The charge symmetry of this model plays a crucial role in the determination of its universality class and in the ability to obtain analytic solutions. This raises the question of how the breaking of the symmetry will affect the thermodynamics and critical properties of electrolyte systems. An important extension of the RPM is the charge-asymmetric primitive model, where the sizes of negative and positive particles are the same while absolute values of charges for positive and negative ions are different. Recent Monte Carlo simulationspana02 ; camp99  have revealed that, as charge asymmetry increases, the critical temperature $`T_c`$ of the gas-liquid transition decreases, while the critical density $`\rho _c`$ grows. However, most of the current theories give different predictions. Simple DH theory and the mean-spherical approximation (MSA) both predict no dependence on the asymmetry parameter $`z`$.levin96 ; MSA  In symmetric Poisson-Boltzmann and modified Poisson-Boltzmann integral equation methodssabir  the charge asymmetry hardly changes $`T_c`$ , while $`\rho _c`$ increases. However, the absolute values of the critical parameters are very different from Monte Carlo estimates. A field theoretical approach by Netz and Orlandnetz  predicts large increases in critical temperatures and, similarly, a large decrease in critical densities, in sharp contrast with computer simulations trends. The only theory that produces reasonable results for the effect of charge asymmetry on thermodynamics and criticality of ionic systems, as judged by comparison with Monte Carlo computer simulations, is the DH approach augmented by Bjerrum cluster formation and cluster-ion interactions (DHBjCI).banerjee02

Lattice models, such as the Ising model, have played an important role in understanding criticality in non-ionic systems. In recent years, lattice models have also attracted the attention of researchers as a tool for investigating thermodynamics and criticality in Coulomb systems.fisher93 ; dickman99 ; ciach01 ; kobelev02 ; brognara02 ; anisotropy  A systematic study of electrolytes on lattices, which utilizes the Debye-Hückel approach, has been presented recently.kobelev02  In this work the thermodynamics of a $`d`$-dimensional system of equal numbers of positive and negative ions, i.e., a charge-symmetric lattice RPM, has been investigated. Specific calculations for Coulomb systems on three-dimensional lattices, which included ion pairing and ion-dipole interactions, predicted a gas-liquid phase separation at low densities. However, taking into account the lattice symmetry yielded a different scenario — the phase diagrams of electrolytes on simple cubic and body-centered cubic lattices show order-disorder phase transitions with a tricritical point, while gas-liquid phase separation is suppressed. The introduction of charge asymmetry in lattice models of ionic fluids tends to suppress the possibility of order-disorder phase transitions, and gas-liquid phase coexistence reappears, although the position of critical point may change.

In this paper, we present a thermodynamic investigation of charge-asymmetric lattice models of electrolytes. By explicitly including the clustering of oppositely charged particles and ion-cluster interactions, we obtain phase diagrams for 2:1 and 3:1 lattice electrolytes, and we locate the critical point for the 4:1 ionic system. Our results accord well with the trend obtained in recent Monte Carlo simulations and the continuum DHBjCI theory.banerjee02  The paper is organized as follows. In Sec.II we present an overview of our thermodynamic approach to multicomponent charged species mixtures and we outline the pure Debye-Hückel theory. The full theory, which accounts for charged and neutral cluster formation and their interactions with the residual ions, for a $`2:1`$ system is presented in detail in Sec.III. Sec.IV describes the general scheme of thermodynamic calculations for $`3:1`$ and $`4:1`$ lattice electrolytes. Finally, the results are discussed in Sec.V and our conclusions are given in Sec.VI.

## II Lattice Debye-Hückel Theory of Charge-asymmetric Electrolytes

### II.1 Thermodynamic Overview

Consider a system of charged particles on a three-dimensional simple cubic lattice with a unit cell length $`a`$, which initially has $`N_0`$ ions carrying a charge $`zq`$ and $`zN_0`$ ions with a charge $`q`$, i.e., the total number of ions is $`N=(z+1)N_0`$. Because of the electrostatic interactions, ions with opposite charges tend to form clusters. As a result, there will be many species present in the system: dimers, trimers, etc., with respective charges $`(z+1)q,(z+2)q,..,0`$. If the number of particles of type $`i`$ is given by $`N_i`$, then we define $`\rho _i=N_i/V`$ and $`\rho _i^{}=\rho _i\mathrm{v}_0`$ to be the number density and the reduced density of the $`i`$-th species, with $`\mathrm{v}_0=a^3`$ being the unit lattice cell volume.

The Helmholtz free energy is central for the determination of the thermodynamic behavior of charge-asymmetric lattice electrolytes. It can be approximated by summing consecutively the free energies describing the interactions between different species,

$$F=F^{Id}+\underset{i}{}F_i,$$
(1)
where $`F^{Id}`$ is the ideal lattice gas (entropic) term and $`F_i`$ is the electrostatic energy of the $`i`$-th species. Once the reduced free energy density $`\overline{f}F/k_BTV`$ is known, the reduced chemical potentials for every component $`\overline{\mu }_i\mu /k_BT`$ can be computed via

$$\overline{\mu }_i=\overline{f}/\rho _i.$$
(2)
Finally, the reduced pressure is given by

$$\overline{p}p/k_BT=\overline{f}+\underset{i}{}\rho _i\overline{\mu }_i.$$
(3)
Then the possible phase equilibria are defined by matching pressures and chemical potentials for each component in different phases.

In multicomponent systems with charged particles it is the electrochemical potential that must be equal in coexisting vapor ($`v`$) and liquid ($`l`$) phases,banerjee02  namely,

$$\mu _{i,v}+q_i\varphi _v=\mu _{i,l}+q_i\varphi _l,$$
(4)
where $`\varphi _{v(l)}`$ is the electrostatic potential in the corresponding phase, where, in general, there is a nonzero Galvani potential difference $`\mathrm{\Delta }\varphi =\varphi _v\varphi _l`$ between the phases.galvani  However, for calculating phase equilibria in multicomponent systems of charged particles, it is more convenient to use the single-component thermodynamic picture.banerjee02  Since every thermodynamic phase is electroneutral, the multicomponent system with $`N=(z+1)N_0`$ ions can be viewed as a single-component system of $`N_0`$ molecules, each of them consisting of one negative ion and $`z`$ positive ions. Then phase equilibrium between the liquid and vapor at temperature $`T`$ is ensured by

$$P(T,\rho _v)=P(T,\rho _l)\text{and}\mu (T,\rho _v)=\mu (T,\rho _l),$$
(5)
where $`\rho _v`$ and $`\rho _l`$ are the overall particle densities in gas and liquid phases, respectively, while $`\mu =\mu _{}+z\mu _+`$. The pressure in each phase can still be calculated using Eq. (3). This approach accounts for the electroneutrality of each phase, and utilizes only one chemical potential, which significantly simplifies calculations of phase diagrams.

### II.2 Pure DH theory of charge-asymmetric lattice electrolytes

As a first approximation, assume that there is no clustering between oppositely charged particles, and that only free ions are present in the system. The free energy density can be written as $`\overline{f}=\overline{f}^{Id}+\overline{f}^{DH}`$, where the first term is the entropic ideal gas contribution, which is given by

$$\overline{f}^{Id}=\frac{\rho _+^{}}{v_0}\mathrm{ln}\rho _+^{}\frac{\rho _{}^{}}{v_0}\mathrm{ln}\rho _{}^{}\frac{1\rho ^{}}{v_0}\mathrm{ln}(1\rho ^{}).$$
(6)
The subscripts “+” and “–” denote positive and negative ions, respectively. Owing to overall electroneutrality in the system, the densities of free ions are related to each other by

$$\rho _+^{}=z\rho ^{}/(1+z),\rho _{}^{}=\rho ^{}/(1+z).$$
(7)
The second term in the free energy density is the electrostatic contribution $`\overline{f}^{DH}=\overline{f}_+(\rho _+)+\overline{f}_{}(\rho _{})`$, which is the result of ion-ion Coulombic interactions. By solving the lattice version of the linearized Poisson-Boltzmann equation for the electrostatic potential, and subsequently applying the Debye charging procedure, as was done for the charge-symmetric lattice model of electrolytes ,kobelev02  it can be easily shown that

$$\overline{f}^{DH}=\frac{1}{12v_0}_0^{x^2}\left[P(1)P\left(\frac{6}{x^2+6}\right)\right]d(x^2),$$
(8)
where

$$P(\zeta )=\frac{1}{(2\pi )^3}_\pi ^\pi _\pi ^\pi _\pi ^\pi \frac{d𝐤}{1\frac{\zeta }{3}(\mathrm{cos}k_1+\mathrm{cos}k_2+\mathrm{cos}k_3)}$$
(9)
is the integrated lattice Green’s function for the simple cubic lattice,katsura71 ; joyce  and $`x^2=\kappa ^2a^2`$ with the reciprocal squared Debye screening length given by

$$\kappa ^2=\frac{4\pi }{Dk_BT}(\rho _+q_+^2+\rho _{}q_{}^2)=\frac{4\pi zq^2\rho }{Dk_BT}.$$
(10)
Now the chemical potential for each type of ion can be easily calculated through Eq. (2), which yields

$$\overline{\mu }_i=\mathrm{ln}\left(\frac{\rho _i^{}}{1\rho ^{}}\right)+\frac{\pi }{3T^{}}\left[P\left(\frac{6}{x^2+6}\right)P(1)\right],$$
(11)
where we have defined the reduced temperature by

$$T^{}=\frac{Dk_BTa}{zq^2},$$
(12)
which is related to the reduced density as $`\rho ^{}=x^2T^{}/4\pi `$. Then for the chemical potential of the neutral “molecule,” consisting of one negative and $`z`$ positive ions, we obtain

$$\overline{\mu }=(1+z)\mathrm{ln}\left(\frac{\rho ^{}}{1\rho ^{}}\right)+(1+z)\frac{\pi }{3T^{}}\left[P\left(\frac{6}{x^2+6}\right)P(1)\right]+C(z),$$
(13)
where the density-independent constant is $`C(z)=\left[z\mathrm{ln}z(1+z)\mathrm{ln}(1+z)\right]`$. Having the chemical potentials and free energy, we can calculate the pressure using Eq. (3) to obtain

$$\overline{p}v_0=\mathrm{ln}(1\rho ^{})+\frac{1}{12}\left[x^2P\left(\frac{6}{x^2+6}\right)_0^{x^2}P\left(\frac{6}{x^2+6}\right)d(x^2)\right].$$
(14)
Note that the expression for the pressure is independent of the charge asymmetry, and the chemical potential of the “molecule” is $`(1+z)/2`$ times the value for the 1:1 electrolytes (if we neglect the constant $`C(z)`$ which does not affect the thermodynamics of the system). Then, the predicted phase separation and the critical point,

$$T_c^{}=0.1018,\rho _c^{}=0.0996,$$
(15)
are the same as for the charge-symmetric lattice electrolytes. This result fully agrees with the pure DH theory in continuum space,banerjee02  which also predicts no change in the critical parameters for charge-asymmetric electrolytes.

## III Lattice Debye-Hückel Theory with Clustering and Cluster-Ions Interactions for 2:1 Electrolytes

### III.1 Bjerrum clustering

One of the main deficiencies of the pure DH theory, which describes the system of free positive and negative ions, is the total neglect of the ion clustering. Oppositely charged particles attract each other and form clusters in order to reduce the free energy. This process significantly decreases the number of free ions. In charge-symmetric lattice electrolytes in the critical region the number of neutral ion pairs is 2-3 times larger than the number of free ions,kobelev02  and one can expect that for charge-asymmetric ionic systems this ratio should be even larger because of the larger free energy gain for clusters with the increase of the parameter $`z`$.

The importance of ion pairing in charged particles systems was first recognized by Bjerrum.bjerrum  In his original approach a cut-off distance between two oppositely charged ions was introduced to define a bound pair. Later, the definition of a pair became a subject of many discussions.fisher93 ; levin96  However, as was shown in a careful analysis by Levin and Fisher, levin96  the precise value of this cut-off distance has little influence on the critical parameters and coexistence curves (less than 0.5%). Meanwhile, the further problem of calculation of the electrostatic energy of the dipole particle poses another, new technical difficulty, since the ion pair does not possess spherical symmetry and the problem cannot be solved exactly analytically or numerically. By approximating the ideal bispherical exclusion zone by a symmetrically centered sphere, Levin and Fisherlevin96  succeeded in obtaining a precise, numerically tractable solution. For three- or four- particle clusters in charge-asymmetric continuum ionic fluids, the same strategy also yields reasonable results, banerjee02  although the technical complexity of computations increases significantly.

In lattice systems the situation is intrinsically simpler because of the discrete rather than continuous symmetry. In $`1:1`$ electrolytes it allows one to define clearly an ion pair as two oppositely charged particles sitting on neighboring lattice sites.kobelev02  Similarly, for charge-asymmetric lattice electrolytes we can easily define different cluster configurations. These particles may be viewed as independent chemical species, and the processes of clustering can be considered as a set of chemical reactions between them.

To illustrate our approach, consider a $`2:1`$ lattice ionic fluid in which clusters are allowed to form. Then, following Ref.banerjee02 , in addition to the free positive particles with charge $`+q`$ and free negative ions with charge $`2q`$, we suppose it suffices to consider three basic types of clusters: specifically we include dimers carrying a charge $`q`$ (with the number density $`\rho _2^{}`$) and two possible configurations of neutral trimers, linear and angular, with number densities $`\rho _{3a}^{}`$ and $`\rho _{3b}^{}`$, respectively: see Fig.1.

In the spirit of Bjerrum’s approach, we assume first that neutral clusters do not interact with charged particles. Then these neutral particles contribute only to the ideal free energy, and the total free energy density is given by

$$\overline{f}=\overline{f}^{Id}(\rho _+,\rho _{},\rho _2,\rho _{3a},\rho _{3b})+\overline{f}_+^{El}+\overline{f}_{}^{El}+\overline{f}_2^{El}.$$
(16)
To calculate the free energy of the charged dimers $`f_2^{El}`$, we adopt a view of them as a combination of neutral symmetric $`(+,)`$ dipoles, which occupy two neighboring sites, and single-site particles (with the charge $`q`$) which sit on the top of negative part of the dipoles, i.e., $`(+,2)=(+,)+()`$. Then the electrostatic free energy densities for the free ions and charged dimers are given by the corresponding expressions from the pure DH theory, although with a different inverse squared Debye screening length, namely,

$$\kappa ^2=\frac{4\pi }{Dk_BT}(\rho _+q^2+4\rho _{}q^2+\rho _2q^2).$$
(17)
To determine the ideal gas (entropic) contribution to the free energy one needs to know the corresponding densities and chemical potentials. As shown in Fig.1, we identify five different species in the $`2:1`$ lattice Coulombic system. Let us define $`\alpha `$ as a set of free positive and negative ions, $`\beta `$ as a set of negatively charged dimers and free positive ions, $`3a`$ as the linear neutral trimers, and $`3b`$ as the neutral angular trimers. Then the clustering in these system can be described by three chemical reactions

$`\alpha `$ $`\stackrel{K_1}{}\beta ,`$ (18)
$`\beta `$ $`\stackrel{K_{2a}}{}3a,`$ (19)
$`\beta `$ $`\stackrel{K_{2b}}{}3b,`$ (20)
with the corresponding equilibrium constants $`K_1`$,$`K_{2a}`$ and $`K_{2b}`$. The chemical equilibrium in the system is described by the following relations between chemical potentials

$$2\mu _++\mu _{}=\mu _++\mu _2=\mu _{3a}=\mu _{3b}.$$
(21)
The chemical potential for each particle can be computed by first separating the entropic and electrostatic parts, i.e., $`\overline{\mu }_i=\overline{\mu }_i^{Id}+\overline{\mu }_i^{El}`$, where the latter is determined simply by $`\overline{\mu }_i^{El}=\overline{f}^{El}/\rho _i`$. For the ideal-gas parts of the chemical potentials the situation is more complex. For free single-site ions the chemical potentials can be easily found applying the potential distribution theorem widom  or from simple entropic considerations,

$`\overline{\mu }_+^{Id}`$ $`=`$ $`\mathrm{ln}\rho _+^{}\mathrm{ln}(1\rho _+^{}\rho _{}^{}2\rho _2^{}3\rho _{3a}^{}3\rho _{3b}^{}),`$ (22)
$`\overline{\mu }_{}^{Id}`$ $`=`$ $`\mathrm{ln}\rho _{}^{}\mathrm{ln}(1\rho _+^{}\rho _{}^{}2\rho _2^{}3\rho _{3a}^{}3\rho _{3b}^{}).`$ (23)
However, for dimers and trimers there are no similar exact expressions. Nevertheless, these entropic contributions to the chemical potentials can be estimated by applying the Bethe approximation, nagle  which has been successful in the DH theory of $`1:1`$ lattice electrolytes. kobelev02  Note, that in lattice electrolyte systems phase transitions predominantly take place in the low-density regimes, and thus the error in using the Bethe approximation is expected to be small.

Defining the activities of the particles via $`z_i`$, for every species we have

$$z_i=\zeta _i/\mathrm{\Lambda }_i^{3n}e^{\mu _i},$$
(24)
where $`n`$ is the number of sites occupied by the particle, $`\zeta _i`$ is the corresponding internal partition function and $`\mathrm{\Lambda }_i`$ is the de Broglie wavelength (see Ref.levin96 ). For the latter the equality holds $`\mathrm{\Lambda }_+=\mathrm{\Lambda }_{}=\mathrm{\Lambda }_2=\mathrm{\Lambda }_{3a}=\mathrm{\Lambda }_{3b}`$. We refer to the Appendix A for a detailed calculation of the dimers and trimers activity using the Bethe approximation method, and we give here only the final expressions for activities,

$`z_2`$ $`=`$ $`{\displaystyle \frac{\rho _2^{}(1\frac{1}{3}\rho _2^{})}{(1\rho _+^{}\rho _{}^{}2\rho _2^{}3\rho _{3a}^{}3\rho _{3b}^{})^2}}e^{\overline{\mu }_2^{El}},`$ (25)
$`z_{3a}`$ $`=`$ $`{\displaystyle \frac{\rho _{3a}^{}(1\frac{1}{2}\rho _2^{})^2}{(1\rho _+^{}\rho _{}^{}2\rho _2^{}3\rho _{3a}^{}3\rho _{3b}^{})^3}}e^{1/(2\rho _{3a})}e^{\overline{\mu }_{3a}^{El}},`$ (26)
$`z_{3b}`$ $`=`$ $`{\displaystyle \frac{\rho _{3b}^{}(1\frac{1}{2}\rho _2^{})^2}{64(1\rho _+^{}\rho _{}^{}2\rho _2^{}3\rho _{3a}^{}3\rho _{3b}^{})^3}}e^{1/(2\rho _{3b})}e^{\overline{\mu }_{3b}^{El}}.`$ (27)
Once the activities are known, the chemical potentials are derived by utilizing Eq.(24) for each species, and the free energy is obtained by integrating the chemical potentials. Finally, we arrive at the chemical potential of an $`\alpha `$-“molecule” (which consists of two positive and one negative ions) in the form

$$\overline{\mu }_\alpha ^{DHBj}=\mathrm{ln}\rho _+^{}+2\mathrm{ln}\rho _{}^{}3\mathrm{ln}(1\rho _+^{}\rho _{}^{}2\rho _2^{}3\rho _{3a}^{}3\rho _{3b}^{})+\overline{\mu }^{El},$$
(28)
and the reduced pressure

$`\overline{p}^{DHBj}v_0`$ $`=`$ $`\mathrm{ln}(1\rho _+^{}\rho _{}^{}2\rho _2^{}3\rho _{3a}^{}3\rho _{3b}^{})+3\mathrm{ln}(1\rho _2^{}/3)+3\mathrm{ln}(1\rho _{3a}^{}/2)`$ (29)
$`+`$ $`3\mathrm{ln}(1\rho _{3b}^{}/2){\displaystyle \frac{\rho _{3a}^{}}{2\rho _{3a}^{}}}{\displaystyle \frac{\rho _{3b}^{}}{2\rho _{3b}^{}}}+\overline{p}^{El}v_0,`$
The electrostatic part of the chemical potential and the pressure have the same form as in the pure DH theory (14), namely,

$`\overline{\mu }^{El}`$ $`=`$ $`{\displaystyle \frac{\pi }{3T^{}}}\left[P\left({\displaystyle \frac{6}{x^2+6}}\right)P(1)\right],`$ (30)
$`\overline{p}^{El}`$ $`=`$ $`{\displaystyle \frac{1}{12}}\left[x^2P\left({\displaystyle \frac{6}{x^2+6}}\right){\displaystyle _0^{x^2}}P\left({\displaystyle \frac{6}{x^2+6}}\right)d(x^2)\right],`$ (31)
where $`x^2=\kappa ^2a^2`$ with $`\kappa ^2`$ given by (17).

In terms of activities, the chemical equilibrium conditions (21) can be presented in the form of laws of mass action

$`z_+^2z_{}K_1`$ $`=`$ $`z_2z_+,`$ (32)
$`z_2z_+K_{2a}`$ $`=`$ $`z_{3a},`$ (33)
$`z_2z_+K_{2b}`$ $`=`$ $`z_{3b},`$ (34)
As can be seen from the chemical equilibrium conditions (21) and the definition of activities (24), the association constants in (32)-(34) are equal to the internal partition functions of dimers and trimers,

$$K_1=\zeta _1(T),K_1K_{2a}=\zeta _{3a}(T),K_1K_{2b}=\zeta _{3b}(T).$$
(35)
The above definitions of the association constants lead to the expressions

$$K_i=v_0c_i\underset{j=1}{\overset{n}{}}e^{u_{ji}/k_BT}$$
(36)
where $`u_{ji}=q_j\phi _j`$ is the electrostatic energy of $`j`$th ion of a cluster particle $`i`$ in the potential $`\phi _j`$, which is due to the other ions entering the multimer cluster. The coefficient $`c_i`$ in Eq.(36) is an entropic factor which takes into account all allowable configurations of the cluster. The electrostatic potentials $`\phi _j`$ can be determined by solving the lattice version of the linearized Poisson-Boltzmann equation

$$\mathrm{\Delta }\phi _j=\frac{4\pi }{Dv_0}\underset{kj}{}q_k\delta (𝐫_k),$$
(37)
in which $`q_k`$ and $`𝐫_k`$ are the charge and the position of the $`k`$th ion entering the multimer cluster. Since lattice Coulomb potentials can be calculated exactly in numerical terms,kobelev02  we obtain the following expressions for association constants

$`K_1`$ $`=`$ $`6v_0\mathrm{exp}\left[{\displaystyle \frac{1.08152}{T^{}}}\right],`$ (38)
$`K_{2a}`$ $`=`$ $`v_0\mathrm{exp}\left[{\displaystyle \frac{0.812033}{T^{}}}\right],`$ (39)
$`K_{2b}`$ $`=`$ $`4v_0\mathrm{exp}\left[{\displaystyle \frac{0.734737}{T^{}}}\right].`$ (40)
Substituting (38)-(40) and the expressions for activities (25)-(27) into the expressions for chemical equilibrium (32)-(34) yields a set of equations, which define implicitly the dimer and trimer densities $`\rho _2^{}`$, $`\rho _{3a}^{}`$ and $`\rho _{3b}^{}`$ in terms of the monomer densities $`\rho _+^{}`$ and $`\rho _{}^{}`$. The electroneutrality of the system requires that $`\rho _+^{}=2\rho _{}^{}+\rho _2^{}`$. Then the chemical potential and the pressure can be expressed in terms of only one variable, the total reduced density $`\rho ^{}=\rho _+^{}+\rho _{}^{}+2\rho _2^{}+3\rho _{3a}^{}+3\rho _{3b}^{}`$, which allows for the construction of the coexistence curve: see Fig.2.

As in the case of charge-symmetric continuum levin96  and latticekobelev02  electrolytes, the coexistence curve in the DHBj approximation has an unphysical banana-like shape. This is related to the fact that, as the temperature is lowered, the number of neutral clusters quickly grows, and this depletes the number of free charges present in the system. Since it is the density of free charges that plays the role of the order parameter and governs the gas-liquid transition, the phase separation takes place at higher overall densities. At this level of approximation, neutral clusters are electrically inactive, and hence contribute only to the hard-core part of the free energy, which is the same for both phases. Thus, DHBj theory simply superimposes the pressure of ideal gas of clusters on the pure DH pressure, and both sides of the coexistence curve shift to higher densities by equal amounts. The critical density is now substantially higher, $`\rho _c^{}0.0807`$, while the critical temperature slightly decreases, $`T_c^{}0.099.`$

### III.2 Cluster-Ion interactions

The predictions of DHBj theory for lattice ionic systems are thermodynamically unreasonable and should be corrected by taking into account the effects of interactions between multimers in clusters and free ions.levin96 ; kobelev02  As shown earlier for continuum and lattice electrolytes,levin96 ; kobelev02  these solvation effects eliminate the unphysical banana-shape phase coexistence curves. It is reasonable to expect that these cluster-ion interactions are even more important in charge-asymmetric ionic systems due to a larger fraction of neutral clusters in equilibrium with free ions.

The exact calculations of interactions between the multimers and free ions are very complicated. Instead we use a reasonable assumption to obtain closed analytic expressions. We approximate the neutral clusters as a set of overlapping non-interacting symmetric dipoles. For example, the neutral trimers can be viewed as a combination of two dipoles, i.e., $`(+,2,+)=(+,)+(,+)`$. The charged clusters, as we discussed above, are approximated as symmetric dipoles overlapping with free ions, i.e., $`(+,2)=(+,)+()`$. We expect that differences between our approximation and the exact cluster-ion interaction contribution to the free energy to be small since the corrections correspond to higher order dipole-dipole interactions, and thus can be neglected.

According to our approximation the reduced density of dipoles from all cluster particles in the system is equal $`\rho _2^{}+2\rho _{3a}^{}+2\rho _{3b}^{}`$. Since the energy of solvation of a single free dipole in the lattice electrolyte system is known exactly,kobelev02  the cluster-ion interactions contribute to the free energy density as

$$\overline{f}^{CI}=(\rho _2^{}+2\rho _{3a}^{}+2\rho _{3b}^{})\frac{\pi q^2a^2}{21Dk_BTv_0^2}\left[\frac{x^2}{2}+\frac{1}{x^2}_0^{x^2}G(x^2)d(x^2)\right],$$
(41)
where

$$G(x^2)=\frac{x^2(x^2+7)}{x^2+6}P\left(\frac{6}{x^2+6}\right).$$
(42)
Then the corresponding contributions to the chemical potentials and pressure are given by

$$\overline{\mu }_\alpha ^{CI}=\frac{2\pi ^2(\rho _2^{}+2\rho _{3a}^{}+2\rho _{3b}^{})}{21(T^{})^2}\left[\frac{1}{2}\frac{1}{x^2}G(x^2)+\frac{1}{x^4}_0^{x^2}G(x^2)d(x^2)\right]$$
(43)
and

$$\overline{p}^{CI}v_0=\overline{\mu }_\alpha ^{CI}x^2T^{}/4\pi ,$$
(44)
which must be added to the values provided by Eqns. (28) and (29) in order to obtain the complete expressions for the chemical potential and pressure.

The full DHBjCI theory predicts a phase separation as exhibited in Fig.2. As in the charge-symmetric lattice and continuum electrolytes,levin96 ; kobelev02  taking into account the solvation of clusters by the residual free ions eliminates the unphysical banana-like shape of the coexistence curve. The critical parameters are now given by

$$T_c^{}=0.08735,\rho _c^{}=0.05001,Z_cp_c/\rho _ck_BT=0.2431.$$
(45)
Comparison with the critical parameters of 1:1 lattice ionic system indicates that the critical temperature is about 11$`\%`$ lower, while the critical density is 1.67 times higher. These trends, the decrease in the critical temperature and inrease in the critical density, are in agreement with the results of continuum calculations for charge-asymmetric electrolytes.banerjee02

It is interesting to estimate the relative amounts of different species in the critical region. In terms of molar fractions $`y_i=n_i\rho _i^{}/\rho ^{}`$ ($`n_i`$ is the size of the particle), our theoretical approach predicts

$`y_+`$ $`=`$ $`0.1002,y_{}=0.0388,y_2=0.0452,`$ (46)
$`y_{3a}`$ $`=`$ $`0.0075,\mathrm{and}y_{3b}=0.8083.`$
As expected, the fraction of neutral trimers significantly exceeds the fraction of free charges. Also, the number of charged dimers is very low, which can be attributed to their propensity to combine with free ions to form energetically more favorable neutral trimers. What is surprising, at first glance, is that the ratio of linear trimers to angular trimers is much less than unity. Indeed, because electrostatic interactions favor a linear arrangement of the ions in a trimer cluster, these clusters are more energetically stable and should prevail over angular trimers. However, these naive arguments do not take into account the entropic considerations. First of all, the angular trimers have more different possible configurations than the linear trimers. In addition, they are more compact and can be packed more densely. As a consequence, there are more ways to arrange angular trimers on the lattice and they dominate trimer clusters.

## IV Lattice DHBjCI Theory for 3:1 and 4:1 Electrolytes

The thermodynamic calculations for $`z=3`$ asymmetric lattice electrolytes can be performed following the method presented in full detail in Sec.III. As shown in Fig.3, there are seven basic clusters to be considered in this system, namely, single-site positive and negative ions, dimers with the charge $`2q`$, linear and angular trimers (with the charge $`q`$), and two types of neutral tetramers, which are connected by the network of six chemical reactions. Note that we do not consider clusters with bonds between the same charges. For example, we assume that charged trimers have the configurations (+,3$``$,+), but not (+,+,3$``$). Indeed, these particles are less stable and their contributions to free energies can be neglected.

In addition, as in the case of 2:1 lattice Coulombic systems, we view the clusters as combinations of noninteracting (+,$``$) symmetric dipoles and single-site charges. The activities of multimers are calculated by applying the Bethe approximation method outlined in the Appendix. These approximations allow us to calculate the free energy of the system as a sum of the corresponding entropic, electrostatic and cluster solvation terms for each particle. The resulting phase coexistence curve is given in Fig.4. The critical parameters are

$$T_c^{}=0.0688,\rho _c^{}=0.0847.$$
(47)
The molar fractions of different clusters in the critical region are

$`y_+`$ $`=`$ $`0.08872,y_{}=0.00347,y_2=0.0088,`$ (48)
$`y_{3a}`$ $`=`$ $`0.0033,y_{3b}=0.2049,y_{4a}=0.3043,y_{4b}=0.3864.`$
As expected, neutral clusters again dominate.

Analogous thermodynamic calculations can be done for 4:1 lattice electrolytes. However, the number of different types of clusters and the number of chemical reactions between them become fairly large, and the full thermodynamic analysis is difficult to complete. Instead, we focus on the critical region of the system where calculations can be completed. The resulting critical parameters are

$$T_c^{}=0.060,\rho _c^{}=0.148.$$
(49)
The relative density of all neutral clusters here is about 87$`\%`$.

Full phase diagrams for $`z:1`$ lattice ionic systems ($`z=1,2`$ or 3) calculated using DHBjCI approach are exhibited in Fig.4.

## V Discussion

Our analysis of charge-asymmetric lattice electrolytes using Debye-Hückel theory with Bjerrum clustering and cluster-ion interactions indicates that charge asymmetry strongly influences the thermodynamics, especially in the critical region. It is found that critical temperatures decrease, while critical densities increase with charge asymmetry. Our theoretical predictions for critical parameters are compared in Fig.5 with the results of continuum Monte Carlo simulationscamp99 ; pana02 ; yan02  and with the predictions of the continuum DHBjCI theory. banerjee02

Clearly, the overall consistency between the lattice and continuum DHBjCI approaches and computer simulations results indicates that the Debye-Hückel method correctly captures the physics of phase transitions in these ionic systems. In Coulomb systems with charge asymmetry the formation of clusters is strongly favored. Clusters contribute to electrostatic interactions much less than free ions, and the effective electrostatic energy per particle decreases. However, the phase separation is driven by charged particles, and the temperature is normalized by the strongest electrostatic interaction between the positive and negative ions: see Eq.(12). Therefore, the reduced critical temperature falls.

Cluster formation also significantly depletes the number of free ions. At the critical point, the molar fractions of free single-site positive and negative ions are 21.0$`\%`$, 13.8$`\%`$, 9.2$`\%`$ and 4.2$`\%`$ for 1:1, 2:1, 3:1 and 4:1 lattice electrolyte systems, respectively. As a consequence, larger overall densities are required to achieve the phase separation, and the reduced critical density increases.

At the same time, the relative numbers of neutral particles as a function of charge asymmetry shows a non-monotonic behavior. The molar fractions of neutral clusters are equal to 79% in 1:1 ionic system, 81% in 2:1 electrolytes, 69.1% in 3:1 ionic system and about 88% in 4:1 electrolytes. The decrease of the molar fraction of neutral clusters in the critical region in 3:1 system as compared to 2:1 electrolytes has also been found in continuum calculations.banerjee02  The actual number of neutral clusters depends on two factors: the free energy gain of cluster formation and the temperature (through the equilibrium constant). It is possible that in 3:1 electrolytes the lowering of the critical temperature is not enough to overcome the repulsion between the positive ions in neutral tetramers, while in 4:1 ionic systems the decrease in critical temperature is enough to favor more strongly the formation of neutral clusters.

The absolute values of the critical temperature in our charge-asymmetric lattice models is larger than the Monte Carlo simulation datacamp99 ; pana02  by a factor of $`1.681.94`$, with lesser factors corresponding to higher $`z`$. This agrees with the general feature of lattice models to over-estimate the values of critical temperatures in comparison with continuum models. Note, however, that the simulations are performed with the continuum potential, while our theory employs the lattice Coulombic potentials. In addition, the reduced temperature (12) is normalized by the energy of two ions interaction via continuum or $`1/r`$ Coulombic potential. However, the lattice potential differs from $`1/r`$ at short distances,kobelev02  and normalizing the temperature by the energy of lattice Coulombic interactions of nearest neighbors would lower the temperatures by an additional 8$`\%`$.

The absolute values of the critical densities in our theoretical approach are only 38-69$`\%`$ of the current Monte Carlo estimates.camp99 ; pana02  This fact reflects the discrete nature of the lattices and that not all possible clusters have been taken into account.

## VI Conclusions

We have investigated the effects of charge asymmetry on thermodynamics and critical properties of the lattice ionic systems using Debye-Hückel theory with Bjerrum clustering and cluster-ion interactions. Phase diagrams for $`z:1`$ lattice electrolytes have been obtained for $`z=2`$ and 3, while for $`z=4`$ the location of the critical point was determined. Our results agree well with the Monte Carlo simulations and with the continuum Debye-Hückel theory,banerjee02  predicting that the increase in the charge asymmetry lowers the critical temperature and increases the critical density.

Our theoretical approach may be extended in several directions. In this paper we investigated charge asymmetric electrolyte systems on the simple cubic lattices. It is interesting to consider the lattices with different symmetry such as body-centered cubic and face-centered cubic lattices, for which thermodynamic calculations for charge symmetric systems are already performed.kobelev02  Another interesting question is how the charge asymmetry will affect the ionic systems in two dimensions, where interesting phase transitions may appear.

## Acknowledgments

Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund (Grant No. 37867-G6) for support of this research. The authors also acknowledge the support of the Camille and Henry Dreyfus New Faculty Awards Program (Grant No. NF-00-056). The authors would like to thank Professor M.E. Fisher for critical comments and useful suggestions, and for providing us with a draft manuscript of his work.

## Appendix A The Bethe approximation for multimers

Consider a system of $`N_m`$ multimers (each particle occupies $`m`$ lattice sites), which interact with each other only via on-site exclusion, on a simple cubic lattice with the total number of sites $`N`$. We can calculate the entropy of the multimer distribution following a direct method of counting probabilities in the Bethe approximation employed by Naglenagle . (For other methods of calculating the partition function for lattice systems see Refs.kikuchi ; burley .) The Bethe approximation does not take into account any correlations around cycles and thus it is exact for lattices with only tree-like paths (Bethe lattices). However for $`3D`$ cubic lattices this approximation gives fairly good results.

Consider, first, a system consisting only of dimers. Then the number of ways to arrange $`N_2=N\rho _2^{}`$ dimers on the simple cubic lattice can be approximated bynagle

$$W_2(\rho _2^{})=\left[\left(\begin{array}{c}N\\ 2N\rho _2^{}\end{array}\right)6^{2N\rho _2^{}}\right]\left[\left(\frac{2\rho _2^{}}{6}\right)^{1/2}\left(1\frac{2\rho _2^{}}{6}\right)^{5/2}\right]^{2N\rho _2^{}}\left[\left(1\frac{2\rho _2^{}}{6}\right)^{6/2}\right]^{N(12\rho _2^{})}.$$
(50)
Here the first square bracket represents the total number of ways to arrange $`N\rho _2^{}`$ dimer vertices on the lattice, taking into account six possible orientations for each dimer, provided one its vertex is fixed. The second bracket gives the probability that a dimer vertex configuration is compatible with all its nearest-neighbors’ vertex configurations. Here $`2\rho _2^{}/6`$ is the probability of finding dimer’s edge along certain lattice direction starting from a given point, and $`(12\rho _2^{}/6)^5`$ ensures there is only one dimer at this point. The third square bracket factor is the probability that no dimers occupy $`N(12\rho _2^{})`$ empty lattice sites. The square roots are taken in order not to compute the probabilities twice. The Bethe approximation for dimer’s activity on the sc lattice then yields

$$z=\mathrm{exp}\left(\frac{1}{N}\frac{d}{d\rho _2^{}}\mathrm{ln}W_2(\rho _2^{})\right)=\frac{\rho _2^{}/3(1\rho _2^{}/3)}{(12\rho _2^{})^2}.$$
(51)
Similar calculations can be performed for the systems containing trimers, tetramers or higher order multimers. However, since now not all multimer vertices are equivalent, we have to account for each type of vertex separately.

In particular, consider the system with $`N_3=N\rho _3^{}`$ trimers. Then there are $`3N\rho _3^{}`$ trimer vertices, with $`N\rho _3^{}`$ “center” and $`2N\rho _3^{}`$ “end” vertices. For linear trimers the number of possible orientations is 3 if we fix a “center” vertex and 6 if we fix an “end” vertex. Then the number of ways to arrange linear trimers on the lattice is given by

$`W_3(\rho _3^{})`$ $`=`$ $`\left[\left(\begin{array}{c}N\\ 3N\rho _3^{}\end{array}\right)6^{2N\rho _3^{}}3^{N\rho _3^{}}\right]\left[\left({\displaystyle \frac{3\rho _3^{}}{6}}\right)^{1/2}\left(1{\displaystyle \frac{3\rho _3^{}}{6}}\right)^{5/2}\right]^{2N\rho _3^{}}`$ (54)
$`\times `$ $`\left[\left({\displaystyle \frac{3\rho _3^{}}{6}}\right)^{2/2}\left(1{\displaystyle \frac{3\rho _3^{}}{6}}\right)^{4/2}\right]^{N\rho _3^{}}\left[\left(1{\displaystyle \frac{3\rho _3^{}}{6}}\right)^{6/2}\right]^{N(13\rho _3^{})}.`$ (55)
Again the first square bracket gives the maximum possible number of distinct trimer vertices configurations, the second square bracket accounts for the probability that all “end” vertex configurations are compatible with their nearest-neighbors’ configurations, the third bracket corresponds to the “center” vertices’ compatibility, and the last factor ensures that there are $`N(13\rho _3^{})`$ lattice sites free of trimers. Then the activity of linear trimers is given by

$$z_{3a}=\frac{\rho _{3a}^{}(13\rho _{3a}^{}/6)^2}{(13\rho _{3a}^{})^3}e^{1/(2\rho _{3a}^{})}.$$
(56)
For angular trimers, the only difference is that now there are 24 distinct trimer orientations if we start from the “end” vertices, and 12 orientations if we start from the “center” vertex. This yields

$$z_{3b}=\frac{1}{64}\frac{\rho _{3b}^{}(13\rho _{3b}^{}/6)^2}{(13\rho _{3b}^{})^3}e^{1/(2\rho _{3b}^{})}.$$
(57)
Similarly, the activity of tetramers is given by the following expression

$$z_4=C\frac{\rho _4^{}(12\rho _4^{}/3)^2}{(14\rho _4^{})^4}e^{1/(12\rho _4^{}/3)},$$
(58)
where the coefficient $`C`$ equals to $`1/96`$ for planar tetramers and $`C=1/128`$ for non-planar tetramers.

Figure captions.

Fig.1. Different charged and neutral particles and chemical reactions in 2:1 lattice electrolyte systems.

Fig.2. Phase diagrams for 2:1 lattice electrolyte: (a) pure DH theory, (b) DHBj approximation, (c) DHBjCI theory.

Fig.3. Different charged and neutral particles and chemical reactions in 3:1 lattice electrolyte.

Fig.4. Phase diagrams of (a) 1:1, (b) 2:1, (c) 3:1 lattice electrolytes in DHBjCI theory.

Fig.5. A) Critical temperatures $`T_c^{}`$ and B) critical densities $`\rho _c^{}`$ as functions of charge asymmetry: (a) our lattice DHBjCI predictions, (b) continuum DHBjCI predictions,banerjee02  (c) Monte Carlo simulations.camp99 ; pana02 ; yan02


format: text

text:
# Critical behavior of gravitating sphalerons

## I Introduction

Over the last two decades, substantial effort has been brought to bear on the study of gravitating Yang–Mills fields. This research has resulted in the discovery of numerous solutions to the coupled Einstein–Yang–Mills(–Higgs and/or –dilaton) equations. These solutions include both black hole and regular or particle–like solutions. In addition to confirming the richness of these nonlinear systems, this work has also been helpful in clarifying the standing of the black hole uniqueness theorems and various “no-hair” ideas.

While many of these solutions have been found by solving the appropriate static equations, it was realized early on that understanding the stability of these solutions is important in order to ascribe relative significance to these solutions within the context of some of the no–hair conjectures. The primary means for evaluating stability have been linear perturbation analyses of various static solutions. As a result, many of these static, gravitating Yang–Mills(–scalar) configurations have been found to be unstable to small time–dependent perturbations.

That many of these solutions appear unstable in linear perturbation theory does not necessarily mean that such solutions are without significance. Indeed, it is now widely accepted within the context of gravitational critical phenomena that some of these solutions will have relevance as attractors in the critical collapse of gravitating fields at the threshold of black hole formation . As an example, the Bartnik–McKinnon solutions of the spherically symmetric, static Einstein–Yang–Mills equations are a countably infinite family of regular solutions characterized by the integer number, $`n`$, of zero–crossings of the gauge potential. These solutions are unstable in linear perturbation theory with the $`n^{\mathrm{th}}`$ member of the family having $`n`$ unstable modes.Strictly, this is true if only the radial, gravitational perturbations are excited. If additional components of the gauge potential are also perturbed (i.e. the sphaleron sector), the $`n^{\mathrm{th}}`$ solution will have $`2n`$ unstable modes, the sphaleron sector contributing an additional $`n`$ unstable modes to those from the gravitational sector. Thus, the $`n=1`$ member of this family has a single unstable mode and is a candidate for being a critical solution in the gravitational collapse of a pure Yang–Mills field. Indeed, this is exactly the result found in  where the full evolution equations for the model were solved. The $`n=1`$ Bartnik–McKinnon solution is the critical solution that sits on the threshold between the complete dispersal of the collapsing field and the formation of a finite size black hole (Type I collapse). The value of the unstable mode for this solution then correctly predicts the scaling relation for the lifetime of near–critical solutions.

In addition to these regular solutions, the Einstein–Yang–Mills equations also admit another countably infinite family of solutions, but with a nonzero horizon. These black hole solutions have non-trivial hair outside their horizons and are again characterized by the number, $`n_{\mathrm{BH}}`$, of zero-crossings of the gauge potential. They too are unstable in linear perturbation theory with $`n_{\mathrm{BH}}`$ unstable modes from the gravitational sector. In agreement with expectations, it has been shown that the $`n_{\mathrm{BH}}=1`$ solution is also a critical solution . But this non-abelian or colored black hole solution, rather than separating dispersion and black hole formation as does the $`n=1`$ Bartnik–McKinnon solution, sits on the threshold between two different kinds of dynamical collapse.

Given this behavior, it is natural to conjecture that other configurations of gravitating Yang–Mills fields should likewise exhibit critical phenomena. With that in mind, we consider here the nonlinear evolution of gravitating $`SU(2)`$ sphalerons in Einstein–Yang–Mills–Higgs theory. We provide evidence that there can be three critical transitions in the initial data space. These include the now standard Type I and Type II transitions as well as the transition mentioned above between different kinds of dynamical collapse on which sits a colored or “hairy” black hole as the intermediate attractor. In the process of examining critical gravitational collapse within this system and the formation of such “hairy” black holes as attractors in the black hole regime, we have also confirmed the stability of two additional endstates of collapse. One is a regular, gravitationally bound configuration of Yang–Mills–Higgs field forming a stable “sphaleron star.” The second is a family of stable, hairy black holes different from those that serve as the critical solutions in the black hole regime. The existence of such stable solutions appears to have first been predicted by Maison .

The outline of the remainder of the paper is then as follows. Section 2 summarizes the equations which constitute the full evolution problem and our numerical approach to their solution. Section 3 describes the numerical results, including results of our parameter space searches, the critical solutions and the nature of the stable solutions. Section 4 offers some conclusions and thoughts for future directions.

## II The model

Our starting point in studying the gravitational collapse of configurations of Yang–Mills–Higgs fields is the action

$`S`$ $`=`$ $`{\displaystyle d^4x\sqrt{g}\left[\frac{1}{16\pi G}R\frac{1}{4g^2}F_{}^{a}{}_{\mu \nu }{}^{}F^{a\mu \nu }\frac{1}{2}(D_\mu \mathrm{\Phi })^{}D^\mu \mathrm{\Phi }V(|\mathrm{\Phi }|^2)\right]}`$ (1)
where $`F_{\mu \nu }^a`$ is the Yang-Mills field strength tensor given by

$`F_{\mu \nu }^a`$ $`=`$ $`_\mu A_\nu ^a_\nu A_\mu ^a+ϵ^{abc}A_\mu ^bA_\nu ^c,`$ (2)
$`D_\mu `$ is the gauge covariant derivative whose action on the Higgs doublet, $`\mathrm{\Phi }`$, is

$`D_\mu \mathrm{\Phi }=_\mu \mathrm{\Phi }+A_\mu \mathrm{\Phi },`$ (3)
and the potential $`V`$ is taken to be

$`V(|\mathrm{\Phi }|^2)`$ $`=`$ $`{\displaystyle \frac{\lambda }{4}}(\mathrm{\Phi }^{}\mathrm{\Phi }\eta ^2)^2.`$ (4)
Varying the action with respect to the metric, the field strength and the Higgs field result in the Einstein equations and the curved space Yang-Mills and Higgs equations, respectively. These are

$`{\displaystyle \frac{1}{8\pi G}}G_{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{1}{g^2}}\left[F_{\mu \lambda }^aF^a{}_{\nu }{}^{}{}_{}{}^{\lambda }{\displaystyle \frac{1}{4}}g_{\mu \nu }F_{\alpha \beta }^aF^{a\alpha \beta }\right]`$ (6)
$`+\left(D_\mu \mathrm{\Phi }\right)^{}D_\nu \mathrm{\Phi }{\displaystyle \frac{1}{2}}g_{\mu \nu }\left(D_\lambda \mathrm{\Phi }\right)^{}D^\lambda \mathrm{\Phi }g_{\mu \nu }V(|\mathrm{\Phi }|^2)`$
$`D_\mu F^{a\mu \nu }`$ $``$ $`_\mu F^{a\mu \nu }+ϵ^{abc}A_\mu ^bF^{c\mu \nu }=0.`$ (7)
$`D^\mu D_\mu \mathrm{\Phi }`$ $`=`$ $`\lambda \mathrm{\Phi }(\mathrm{\Phi }^{}\mathrm{\Phi }\eta ^2)`$ (8)
With these general forms for the equations of motion, we make some simplifying assumptions. In particular, we will restrict ourselves to spherically symmetric gravitational collapse and work exclusively with an $`SU(2)`$ gauge group. We also make the assumption that the Higgs field lives in the fundamental representation of $`SU(2)`$. The corresponding flat space version of this theory includes the so–called sphaleron solutions .

Our intent is to solve the full set of nonlinear, evolution equations representing gravitational collapse. In order to do this numerically we must fix both the coordinate freedom and the gauge freedom in our model. There are, of course, numerous possibilities, but we will try to hew fairly closely to related work of others in the field. For the coordinate system, we will work in maximal, areal coordinates. If the general form of the spherically symmetric metric is written as

$`ds^2`$ $`=`$ $`(\alpha ^2+a^2\beta ^2)dt^2+2a^2\beta dtdr+a^2dr^2+b^2r^2d\theta ^2+b^2r^2\mathrm{sin}^2\theta d\varphi ^2`$ (9)
where the metric components depend only on $`t`$ and $`r`$, the choice of areal coordinates amounts to $`b=1`$ while choosing maximal time slices corresponds to the vanishing of the trace of the extrinsic curvature, $`K=0`$.

For the Yang-Mills field, the most general form for a spherically symmetric gauge potential is the Witten ansatz :

$`A`$ $`=`$ $`u\tau _rdt+v\tau _rdr+\left(w\tau _\theta +(\stackrel{~}{w}1)\tau _\varphi \right)d\theta +\left((1\stackrel{~}{w})\tau _\theta +w\tau _\varphi \right)\mathrm{sin}\theta d\varphi `$ (10)
where the $`\tau _i`$ ($`i\{r,\theta ,\varphi \}`$) are the spherical projection of the Pauli spin matrices and form an anti-Hermitian basis for the group $`SU(2)`$ satisfying $`[\tau _i,\tau _j]=ϵ_{ijk}\tau _k`$. With this ansatz for the gauge potential there is some gauge freedom that allows us to simplify its form, namely the potential is invariant under a transformation of the form $`U=e^{\psi (t,r)\tau _r}`$. We can fix some of that gauge freedom by choosing $`v0`$. This choice effectively eliminates the $`r`$ dependance in the above gauge transformation. If we choose to work within the so-called “magnetic ansatz” we can fix the remaining freedom in the following way. It can be shown that in this ansatz, the component $`u`$ is a function only of $`t`$, i.e. it is now pure gauge and can be set to zero as part of our gauge fixing. The remaining fields, $`w`$ and $`\stackrel{~}{w}`$, under the remaining constant gauge transformations are merely sent into linear combinations of each other and hence we can fix the last bit of gauge freedom by setting $`w=0`$. This leaves $`\stackrel{~}{w}`$ as the sole non-zero component of the gauge potential.

Our form for the Higgs field, taken from , is

$`\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(\gamma \mathrm{𝟏}2\psi \tau _r)\left(\begin{array}{c}0\\ 1\end{array}\right)`$ (11)
and though not strictly spherically symmetric, results in a spherically symmetric energy density . We will consider in this work only the case in which $`\psi =0`$. This is not a gauge choice but an additional assumption made merely to simplify the resulting equations and dynamics. A similar thing is done, for instance, in .

With these assumptions, the evolution equations for the Yang-Mills field become

$`\dot{\stackrel{~}{w}}`$ $`=`$ $`{\displaystyle \frac{\alpha }{a}}P+\beta Q`$ (12)
$`\dot{Q}`$ $`=`$ $`\left({\displaystyle \frac{\alpha }{a}}P+\beta Q\right)^{}`$ (13)
$`\dot{P}`$ $`=`$ $`\left(\beta P+{\displaystyle \frac{\alpha }{a}}Q\right)^{}+{\displaystyle \frac{\alpha a\stackrel{~}{w}}{r^2}}\left(1\stackrel{~}{w}^2\right){\displaystyle \frac{\alpha a}{4}}\gamma ^2\left(\stackrel{~}{w}1\right)`$ (14)
while the evolution equations for the Higgs field are given as

$`\dot{\gamma }`$ $`=`$ $`{\displaystyle \frac{\alpha }{a}}G+\beta E`$ (15)
$`\dot{E}`$ $`=`$ $`\left({\displaystyle \frac{\alpha }{a}}G+\beta E\right)^{}`$ (16)
$`\dot{G}`$ $`=`$ $`{\displaystyle \frac{1}{r^2}}\left[r^2(\beta G+{\displaystyle \frac{\alpha }{a}}E)\right]^{}{\displaystyle \frac{\alpha a\gamma }{2r^2}}\left(\stackrel{~}{w}1\right)^2{\displaystyle \frac{\lambda }{2}}\alpha a\gamma \left(\gamma ^22\eta ^2\right),`$ (17)
where, as usual, overdots and primes denote differentiation with respect to $`t`$ and $`r`$ respectively. Both of these sets of evolution equations are supplemented with the first order definitions $`Q=\stackrel{~}{w}^{}`$ and $`E=\gamma ^{}`$ as well as the constraints on the metric components coming from the Einstein equations

$`\alpha ^{\prime \prime }`$ $`=`$ $`\alpha ^{}\left({\displaystyle \frac{a^{}}{a}}{\displaystyle \frac{2}{r}}\right)+{\displaystyle \frac{2\alpha }{r^2}}\left(a^21+{\displaystyle \frac{2ra^{}}{a}}\right)+4\pi G\alpha (S3\rho )`$ (18)
$`a^{}`$ $`=`$ $`a{\displaystyle \frac{1a^2}{2r}}+{\displaystyle \frac{3}{8}}ra^3K_r{}_{}{}^{r}{}_{}{}^{2}+4\pi Gra^3\rho .`$ (19)
$`K_r{}_{}{}^{r}^{}`$ $`=`$ $`{\displaystyle \frac{3}{r}}K_r{}_{}{}^{r}+8\pi G{\displaystyle \frac{1}{a}}[{\displaystyle \frac{PQ}{r^2}}+{\displaystyle \frac{1}{4}}EG]`$ (20)
$`\beta `$ $`=`$ $`\alpha rK_{\theta }^{}{}_{}{}^{\theta }`$ (21)
The matter stress-energy terms in these equations are given by

$`\rho +S_r^r`$ $`=`$ $`{\displaystyle \frac{P^2+Q^2}{a^2r^2}}+{\displaystyle \frac{E^2+G^2}{4a^2}}`$ (22)
$`\rho S_{r}^{}{}_{}{}^{r}`$ $`=`$ $`{\displaystyle \frac{(\stackrel{~}{w}^21)^2}{2r^4}}+{\displaystyle \frac{\gamma ^2(\stackrel{~}{w}1)^2}{8r^2}}+{\displaystyle \frac{\lambda }{16}}(\gamma ^22\eta ^2)^2`$ (23)
$`S_{\theta }^{}{}_{}{}^{\theta }`$ $`=`$ $`{\displaystyle \frac{1}{8a^2}}(G^2E^2)+{\displaystyle \frac{(\stackrel{~}{w}^21)^2}{4r^4}}{\displaystyle \frac{\lambda }{32}}(\gamma ^22\eta ^2)^2`$ (24)
$`j_r`$ $`=`$ $`{\displaystyle \frac{PQ}{ar^2}}{\displaystyle \frac{EG}{4a}}`$ (25)
Boundary conditions are implemented by demanding regularity at the origin and requiring the presence of only outgoing radiation at large distances (see ). The resulting constraints on the metric components require $`\alpha ^{}(t,0)=a^{}(t,0)=\beta (t,0)=K_r^r(t,0)=0`$. The matter fields may satisfy one of two possible regular configurations at the origin. Either, $`\gamma (t,0)=0`$ and $`\stackrel{~}{w}(t,0)=1`$, or $`\gamma ^{}(t,0)=0`$ and $`\stackrel{~}{w}(t,0)=1`$. These two choices correspond to the odd and even node solutions, respectively . In order to find the critical solution we choose the former and look for the solution with a single unstable mode. The outgoing conditions require that at the edge of our grid,

$`P`$ $`=`$ $`Q`$ (26)
$`\dot{Q}`$ $`=`$ $`[({\displaystyle \frac{\alpha }{a}}\beta )Q]^{}`$ (27)
$`G`$ $`=`$ $`{\displaystyle \frac{1}{r}}(\beta {\displaystyle \frac{a}{\alpha }}1)(\gamma \sqrt{2}\eta )`$ (28)
$`\dot{E}`$ $`=`$ $`[(\beta {\displaystyle \frac{\alpha }{a}})(E+{\displaystyle \frac{1}{r}}[\gamma \sqrt{2}\eta ])]^{}.`$ (29)
For the initial pulse, we use a “time-symmetric kink” as in  for the gauge potential, namely

$`w(0,r)`$ $`=`$ $`\left[1+a\left(1+{\displaystyle \frac{br}{s}}\right)e^{2\left(r/s\right)^2}\right]\mathrm{tanh}\left({\displaystyle \frac{xr}{s}}\right)`$ (30)
$`\dot{w}(0,r)`$ $`=`$ $`0`$ (31)
where the parameters $`a`$ and $`b`$ are chosen so $`w(0,0)=1`$ and $`w^{}(0,0)=0`$. The parameters $`x`$ and $`s`$ are the center and width of the kink, respectively.

The Higgs field is initialized as

$`\gamma (0,r)`$ $`=`$ $`\sqrt{2}\eta \mathrm{tanh}{\displaystyle \frac{r}{\delta }}+A_He^{(rr_0)^2/d^2}`$ (32)
$`\dot{\gamma }(0,r)`$ $`=`$ $`0`$ (33)
where the parameter $`\delta `$ is usually set to $`\delta =10`$. This is primarily due to the fact that varying $`\delta `$ does not significantly change the final result of the collapse. As a consequence, we perturb the Higgs field via a gaussian pulse. Similar to the initialization for the Yang–Mills field, the parameters $`A_H`$, $`r_0`$, and $`d`$ which describe the initialization of the Higgs field represent the amplitude, center, and width of the gaussian pulse, respectively. These initial data parameters for the Yang–Mills and Higgs fields will constitute our initial data set and will be used when tuning our evolutions to the critical solutions.

Our numerical approach closely follows that of . We use a uniform grid recognizing that we will not have sufficient resolution to investigate Type II collapse in a completely satisfactory way. Nonetheless, we have indications verifying the existence of Type II behavior in our model. For this paper, therefore, we focus our primary interest on the black hole transition and the dynamics occuring within the black hole regime.

We use an iterative Crank–Nicholson scheme for the evolution equations while for the constraints, we simply integrate outward from the origin. As we want to consider evolutions that extend to the future of black hole formation, our use of maximal slicing is crucial. In our coordinates, the apparent horizon equation is an algebraic relation

$`arK_{}^{\theta }{}_{\theta }{}^{}`$ $`=`$ $`1`$ (34)
We use the same black hole excision technique developed in  and used in . As discussed there, we set a threshold value slightly larger than one such that if $`arK_{}^{\theta }{}_{\theta }{}^{}`$ exceeds that value for certain grid points, we discard those points at future time steps considering them inside the apparent horizon. At the boundary of this region, we need no new boundary conditions for the evolution equations. For those variables solved via constraint equations we either switch to solving an evolution equation subsequent to the formation of a horizon or we “freeze” the variable (e.g. $`\alpha `$) such that it retains the value it had when the horizon formed . As a result, though we can observe matter falling into the horizon, we cannot comment on any dynamics within the apparent horizon as the evolution is effectively frozen for values of the radial coordinate less than the horizon radius. This procedure thus allows us to evolve past the formation of the black hole and thereby investigate such things as the final endstates as well as the critical dynamics in the vicinity of transition regions.

We have tested the resulting code and shown it to be second–order convergent and to conserve mass. It also reproduces the results of  in the limit that the Higgs field and its coupling vanish. Finally, we note that we made extensive use of RNPL (Rapid Numerical Prototyping Language) , a language written expressly to aid the differencing and solution of PDE’s.

## III Results

In attempting to evolve these equations, it quickly became clear that the size of the initial data sets that could be varied is somewhat unwieldy and we had to make choices in order to restrict the possible sets of initial data parameters. Although we have performed numerous evolutions by varying the elements of different sets of initial data parameters, we will focus on the evolution of two sets of parameters to highlight our results. Other sets would appear to give qualitatively similar conclusions. When we evolve these equations, we confirm many of the same aspects that have come to be expected in similar models. However, there are, at the same time, a number of unexpected surprises.

To begin our examination of the dynamics of this model, we consider varying two parameters describing the initialization of the gauge potential, namely $`x`$ and $`s`$, the center and width of the kink, respectively. The amplitude of the perturbing, gaussian pulse for the Higgs field is set to zero, $`A_H=0`$, and the width of the $`\mathrm{tanh}`$ function describing the Higgs is $`\delta =10`$. We call this configuration the “bland” Higgs field. Note that in this section, all the pictured results are for values of the Higgs coupling parameters, $`\eta =0.1`$ and $`\lambda =0.5`$.

On varying the center and width of the Yang-Mills potential, $`\stackrel{~}{w}`$, one finds three distinct regions of the initial data space. These correspond, as in , to generalized Type I collapse, generalized Type II collapse and a “dispersive” region in which no black hole forms. On the boundaries between these regions sit appropriate critical solutions. It is worth noting, however, that in the region in which no black hole forms, we no longer observe the complete dispersal of all the matter fields. Instead, while a majority of the fields do escape to infinity, a nontrivial portion of the fields forms a bound state, or “sphaleron star.” Shortly after formation, this solution oscillates rapidly, but settles down to what appears to be a static solution. Long evolutions with $`t30,000`$ confirm the stability of this solution. The mass of this stable star is, to within a few percent, independent of any of the initial field parameters. The solution and its mass do appear to depend on the coupling parameters $`\eta `$ and $`\lambda `$ . Snapshots of a typical evolution in the “dispersive” regime are shown in Figure 1.

In the region with generalized Type I collapse, we note no significant change in the dynamics from those similarly exhibited in . The collapsing matter forms a black hole with finite mass which, after the residual fields have dispersed to infinity, settles down to the Schwarzschild solution. On the Type I critical line separating black hole formation and the sphaleron star configuration, we find a regular sphaleron as the critical solution analogous to the Bartnik–McKinnon $`n=1`$ solution. An example of the critical behavior at this Type I transition is given in Figure 2 in which a sub–critical and a super–critical evolution are shown.

Within the generalized Type II region, there are some new features. As in , this region is again characterized by the dynamical formation of a black hole with finite mass. As the critical line which separates dispersion (or strictly, sphaleron star formation) from the black hole region is approached, the mass of these black holes begins to decrease such that we interpret the critical transition as Type II. However, the black holes that form away from the critical line after the transient hair has dispersed to infinity do not settle down to Schwarzschild black holes. Instead, the final endstate would appear to be a stable, colored black hole with non–trivial Yang–Mills and Higgs fields outside the event horizon. This, of course, is analogous to the sphaleron star that forms in the no–black–hole region of this system rather than the complete dispersal of the fields seen in . The basic dynamics in this case are illustrated in Figure 3. The collapsing configuration forms a finite mass black hole with a significant portion of the remaining field escaping to infinity. Nevertheless, some of the hair remains behind and within the vicinity of the event horizon. This hair oscillates for some time and eventually settles down to a stable configuration. Evolutions on the order of $`t30,000`$ show no appreciable diminution or instability in the fields.

The mass of the hair in these black hole solutions also seems to be independent of the initial data parameters. Though the radius of the black hole will vary with the initial parameters, the exterior mass remains unchanged. This observation is consistent with and similar to that for the sphaleron stars in which a single stable, regular solution is found throughout the no–black–hole region. In addition, like their regular counterparts, the black hole solutions will depend on the parameters $`\eta `$ and $`\lambda `$. Curiously, the mass of this exterior hair is very nearly the same value as the mass of the sphaleron star. Thus, in one sense, these hairy black hole solutions can be thought of as sphaleron star solutions within which the central density increases to the point that a horizon forms. This is similar in turn to gravitating t’Hooft-Polyakov monopoles. For certain values of the coupling, the monopoles can have a black hole form at their center.

For the solutions near the threshold separating black hole formation and dispersion, we find hints that these are indeed Type II critical solutions and that the black hole mass scales as expected. However, we stress again that our unigrid code is not able to settle this issue definitively and that it awaits additional study.

We are able, though, to consider the transition between the different types of dynamical collapse in the black hole regime. Again, we find a family of critical solutions separating generalized Type I collapse from generalized Type II collapse. These critical solutions are themselves sphaleron black holes parameterized by their horizon radius such that as one moves away from the “triple point” in Figure 5, that radius increases. On the Type II side of this line, near–critical evolutions have dynamics described above with the collapsing configuration forming a black hole of finite mass with non–trivial hair outside. However, as the transition between Type II and Type I is approached, there is an intermediate solution – a hairy black hole – which forms and to which the evolving solution is attracted. This intermediate solution is unstable and eventually collapses. For initial configurations on the Type II side, the collapse is distinctive in that very little of the exterior fields falls into the black hole. Rather, some of it is dispersed to infinity while the remainder reconstitutes in a new and different colored configuration outside the black hole already present. This is shown in Figure 3 for a generic collapse in the Type II regime as well as in the last frames of Figure 4 for a near–critical evolution.

On the Type I side of the critical line, similar near–critical evolutions exhibit the same early time dynamics with the formation of a finite mass black hole and the approach to the intermediate hairy black hole. However, as the critical line is approached, this unstable black hole now collapses and loses most of its hair into the black hole causing it to grow in size. A picture of both collapse dynamics is shown in Figure 4.

So far, our entire description has been within the context of varying two of the initial data parameters that describe the Yang–Mills field. We have chosen to vary the center and width of the Yang–Mills kink, but the same results hold for other Yang–Mills parameters. However, there does appear to be a non–trivial difference in our description of the initial data space if we include a parameter of the Higgs field in the set of two parameters which we vary in order to search for critical behavior.

In the original, so-called “bland” Higgs case, the structure of the phase space is very similar to that found in , namely Type I and II critical collapse separated by the previously discussed colored, critical solution. This is shown in Figure 5.

However, in the event that we turn on the exponential, perturbing pulse for the Higgs field, and look at the phase space including one Yang–Mills parameter and one Higgs parameter, we see a different picture, shown in Figure 6. In this case and for the region of the parameter space explored, there appears to be no critical boundary between Type I collapse and the formation of a sphaleron star. Instead, Type II collapse would seem to border the region in which the regular solution is formed. It would therefore seem that in this region of the initial data space, there is no “triple point.” One might have imagined that as one varied the amplitude of the exponential portion of the Higgs pulse that there would be such a point in each $`xs`$ plane for a given $`A_H`$. However, this does not seem to be the case and is an issue deserving of more investigation.

Finally, we note that the critical solutions separating the two types of collapse and Type I collapse from sphaleron star formation exhibit time scaling as would be expected. As the single, unstable mode characteristic of each critical solution is tuned out, near–critical solutions spend increasing amounts of time as measured by an asymptotic observer on the critical solution. These scaling relations are given by $`T\lambda \mathrm{ln}|pp^{}|`$ where $`\lambda `$ is the characteristic time scale for the collapse of the unstable critical solution. It corresponds to the inverse Lyapounov exponent of the unstable mode. Such scaling relations specific to points on the relevant critical lines are shown in Figures 8 and 7.

## IV Discussion

We have presented evidence for critical phenomena in the gravitational collapse of $`SU(2)`$ sphaleron configurations of Yang–Mills–Higgs fields. In many respects, this collapse is qualitatively similar to that in the Einstein–Yang–Mills system but does have some notable suprises. The critical behavior is seen again in three possible transitions. On each of these transitions sit critical solutions which serve as intermediate attractors for nearby evolutions in the initial data space. Near the critical line separating Type I collapse from the regular solution as well as for the critical line separating the two types of dynamical black hole formation, there are time scaling relations as the near–critical solutions approach the critical solutions. In addition, the mass of the black holes formed in the appropriate region will exhibit a mass gap in crossing these critical lines. Near the critical line separating Type II collapse from the regular endstate region, we have indications that the mass of the black hole scales without a mass gap, but again, due to our unigrid code, we can not settle this conclusively although expectations and indications would bear this out.

Among the surprises in this model are that in certain regions of the initial data space, we find that regular, stable, sphaleron solutions are produced rather than the purely dispersive regime seen in . However, we also have stable, hairy or colored black holes produced in the supercritical, black hole regime. This contrasts again with earlier results in which the endstate in spherical symmetry was always a Schwarzschild black hole. Thus, the transition here between generalized Type I and Type II collapse is a transition between different types of black holes. Within the Type I region, the endstate is always a Schwarzschild black hole with the exterior gauge and Higgs fields either falling into the existing black hole or dispersing to infinity. In the Type II region, the final black holes are stable, colored black holes. It is worth reemphasizing that the colored black hole solutions sitting on the critical line separating types of black hole collapse are not the same as the stable, colored black holes that are the final endstates in the Type II region. This can be seen most easily in Figure 4.

It is also noteworthy that the existence of all the solutions which we find is contingent on the magnetic ansatz within which we have chosen to work. In general, both the regular and colored black hole solutions which we find to be stable endstates of collapse are expected to be unstable based on a linear perturbation analysis . However, such an analysis assumes that both the gravitational and sphaleron sectors in the theory are perturbed. Our evolutions perturb only the gravitational sector. It is reasonable to assume that the stable solutions which we find will become unstable on perturbation of the Yang–Mills gauge field away from the magnetic ansatz. We hope to address this issue in future work.

Another issue for future consideration is the structure of the initial data space. A curiosity of our current results is that the “triple point” would seem to disappear as the Higgs field is varied. As a result, the boundary between regular endstates (sphaleron star formation here) and black hole formation is taken up entirely by a Type II transition.

Nonetheless, given the structure of the initial data space, one can draw an analogy with the gravitating monopole case in which a small black hole can form within a t’Hooft-Polyakov monopole coupled to gravity. This stable object can be rendered unstable above a maximum value of the horizon radius at which point the exterior Yang-Mills-Higgs hair will either fall into the black hole or disperse leaving a final Schwarzschild black hole. A similar thing happens in the current sphaleron case. For example, following a line of constant $`x`$ in Figure 5 that intersects each region we see that as $`s`$ decreases, one can interpret the process in a similar way. A regular solution develops a small, stable black hole at the center which (with decreasing width of the initial Yang-Mills potential, $`s`$) increases in size until the combined sphaleron and black hole system becomes unstable and is replaced with a larger Schwarzschild black hole. As a result, it would be interesting to consider the full dynamical evolution of the gravitating monopole and compare with the sphaleron case reported here.

## Acknowledgements

This research has been supported in part by NSF grants PHY-9900644 and PHY-0139782. EWH would like to thank A. Wang for useful discussions as well as M. Choptuik and R. Marsa for their work on an earlier code which served as the precursor to the one described here.