text:
# The stiffness exponent of two-dimensional Ising spin glasses for non-periodic boundary conditions using aspect-ratio scaling

## I Introduction

Although spin glassesreviewSG  have been studied for more than two decades, finite-dimensional spin glasses are still far from being well understood. This is true also in $`d=2`$ dimensions, where no stable spin-glass phasemcmillan1984 ; bray84 ; rieger1996 ; stiff2d  exists at $`T>0`$. The fact that the transition temperature, $`T_c`$, is zero can be seen by studying domain walls between ground states induced by changing the boundary conditions. For a system of size $`L\times M`$, if the domain wall is forced to run in the $`y`$-direction, the mean domain-wall energy scales like $`M^\theta `$ in the thermodynamic limit, where $`\theta `$ is the stiffness exponent. If $`\theta <0`$, as in two-dimensional spin glasses, no stable phase can exist at finite temperature.

Recentlystiff2d  slightly different values of the stiffness exponent have been found for Gaussian distribution of the bonds for different choices of the boundary conditions. Here, by applying the aspect-ratio methodcarter2002 , we show that this difference is due to corrections to scaling which persisted, at least for one set of boundary conditions, even for the quite large system sizes used in Ref. stiff2d, . The results presented here are consistent with $`\theta =0.287(4)`$ for all boundary conditions. This agrees with the result found in Ref. carter2002,  where periodic boundary conditions were applied in all directions for smaller systems than investigated here, indicating a high degree of universality for the exponent $`\theta `$. Apart from using larger systems and different boundary conditions, we go in another respect beyond the previous work, because aspect ratios, $`RL/M`$, less than unity are studied here. This enables us to verify the proposed scaling behavior in the region $`R<1`$, where the finite-size scaling depends on the boundary conditions.

The model we study consists of $`N=L\times M`$ Ising spins $`S_i=\pm 1`$ on a square lattice with the Hamiltonian

$$=\underset{i,j}{}J_{ij}S_iS_j,$$
(1)
where the sum runs over all pairs of nearest neighbors $`i,j`$ and the $`J_{ij}`$ are quenched random variables. Here, we consider two kinds of disorder distributions:

1. Gaussian with zero mean and variance unity,
2. bimodal, i.e. $`J_{ij}=\pm 1`$ with equal probability.

To study whether an ordered phase is stable at finite temperatures, the following procedure is usually applied mcmillan1984 ; bray84 ; rieger1996 ; kawashima1997 . First a ground state of the system is calculated, having energy $`E_0`$. Then the system is perturbed to introduce a domain wall and the new ground state energy, $`E_0^{\mathrm{pertb}}`$ is evaluated. For each sample, the domain-wall energy is given by $`E_{\mathrm{DW}}=|E_0^{\mathrm{pertb}}E_0|`$. Here, we apply free boundary conditions (bc) in the $`y`$-direction. To study the scaling behavior of the domain-wall energy, we consider two means to introduce domain walls:

1. “P-AP”: First a ground state with periodic bc in the $`x`$-direction is calculated. Then the system is perturbed by introducing antiperiodic bc in that directions, e.g. by reversing one line of bonds parallel to the $`y`$-direction.
2. “F-DW” (called “F-AF” in Ref. carter2002, ): First a ground state with free bc in the $`x`$-direction is calculated. For the perturbed system, we then add extra bonds which wrap around the system in the $`x`$-direction, and have a sign and strength such that they force the spins they connect to have the opposite relative orientation to that which they had in the original ground state.

In both cases, $`M`$ denotes the size of the edge along which the boundary conditions are changed to induce a domain wall (i.e. the size in the $`y`$-direction). The scaling behavior

$$E_{\mathrm{DW}}=M^\theta F(L/M)$$
(2)
has been predictedcarter2002 , with the following limiting forms:

$$F(R)\{\begin{array}{ccc}R^1\hfill & & \text{for }R\mathrm{}\hfill \\ R^{\theta (d1)/2}\hfill & \text{(P-AP)}\hfill & \text{for }R0\hfill \\ R^{\theta (d1)}\hfill & \text{(F-DW)}\hfill & \text{for }R0.\hfill \end{array}$$
(3)
We shall verfiy that the data does follow this scaling form very well.

## II The Algorithm

In greater than two dimensions, or in the presence of a magnetic field, the exact calculation of spin-glass ground states belongs to the class of NP-hard problemsbarahona1982 ; opt-phys2001 . This means that only algorithms with exponentially increasing running times are known. However, for the special case of a planar system without magnetic field, e.g. a square lattice with periodic boundary conditions in at most one direction, there are efficient polynomial-time “matching” algorithmsbieche1980 . The basic idea is to represent each realization of the disorder by its frustrated plaquettestoulouse1977 . Pairs of frustrated plaquettes are connected by paths in the lattice and the weight of a path is defined by the sum of the absolute values of the coupling constants which are crossed by the path. A ground state corresponds to the set of paths with minimum total weight, such that each frustrated plaquette is connected to exactly one other frustrated plaquette. This is called a minimum-weight perfect matching. The bonds which are crossed by paths connecting the frustrated plaquettes are unsatisfied in the ground state, and all other bonds are satisfied.

For the calculation of the minimum-weight perfect matching, efficient polynomial-time algorithms are availablebarahona1982b ; derigs1991 . Recently, an implementation has been presentedpalmer1999 , where ground-state energies of large systems of size $`N1800^2`$ were calculated. Here, an algorithm from the LEDA libraryleda1999  has been applied, which limits the system sizes due to the restricted size of the main memory of the computers which we used. Furthermore, each system border with a free boundary decreases the running time relative to a boundary with pbc. The reason is that frustrated plaquettes near a free boundary may be connected even when they are far apart from each other, increasing the number of possible connections. On the other hand, full pcb cannot be realized with the matching algorithm. Hence, we have limited our system sizes to $`L1024,M16`$ for $`R1`$, and $`L12,M768`$ for $`R<1`$.

## III Results

For Gaussian interactions we have considered systems with P-AP and F-DW boundary conditions with aspect ratios $`R=1/8`$ to $`R=64`$. We took system sizes $`L=1`$ to $`L=12`$ for $`R=1/8`$, $`1/4`$, $`1/2`$, and $`M=1`$ to $`M=16`$ for $`R=1,2,4,\mathrm{},64`$. Typically $`10^5`$ samples were treated per system size $`(L,M)`$, except for few combinations, notably $`R=64`$ where we studied $`10^4`$ samples.

For bimodal interactions we have studied P-AP boundary conditions for $`M=1`$ to $`M=16`$ ($`R=1,2,4,\mathrm{},64`$), with $`10^4`$ samples per system size.

We start with the results for the Gaussian distribution. In Figs. 1 and 2 the scaled average domain-wall energies $`LE_{\mathrm{DW}}`$ are displayed for P-AP and F-DW boundary conditions respectively. From the scaling theorycarter2002  a $`M^{1+\theta }`$ behavior is expected for $`R1`$. Indeed for larger values of $`R`$ and $`M`$ straight lines are visible in the double logarithmic plot. From fitting to algebraic functions, we have obtained for both cases the values of $`\theta `$ as a function of $`R`$; for small aspect-ratio the small system sizes were omitted from the fits. The results are displayed in Table 1.

For the P-AP case, the exponent $`\theta `$ is more or less independent of the system size. By contrast, for the F-DW case a significant increase with system size is observed, with the value seeming to converge near the value of the P-AP case for large $`R`$. This $`R`$-dependence explains why in previous workstiff2d  ($`R=1`$) a difference between P-AP and F-DW has been found. The result found here is compatible with the exponents being equal for both cases. We quote $`\theta =0.287(4)`$, which is a conservative estimate including all values of $`\theta `$ found for the P-AP case. This is compatible with the result of $`\theta =0.282(3)`$ obtainedcarter2002  for systems with full periodic bc.

Part of the data in Ref. carter2002,  corresponds to the F-DW case with pbc in the $`y`$-direction, and we show this here in the final column in Table 1, where $`M12`$ was used. For this case an even stronger dependence of the effective exponent on $`R`$ is observed, indicating that corrections to scaling are even larger in this case. We suggest that this is because the pbc forces the domain wall to take at the top and bottom edges the same positions along the $`x`$-axis, thereby raising its typical energy. In Ref. carter2002,  it was shown that $`\theta _{\mathrm{F}\mathrm{DW}}^\mathrm{P}`$ (called $`\theta _{\mathrm{F}\mathrm{AF}}`$ in Ref. carter2002, ) nevertheless extrapolates, for $`R\mathrm{}`$, to an asymptotic value consistent with that obtained from the present work, where we employed free boundary conditions in the $`y`$-direction.

We have also checked the scaling predictionscarter2002  for $`E_{\mathrm{DW}}`$ explicitly. In Fig. 3 the scaled domain-wall energy $`M^\theta E_{\mathrm{DW}}`$ is plotted for all considered values of $`M`$ as a function of the aspect ratio $`R`$. A very good data collapse can be observed (see e.g. right inset of Fig. 3 for P-AP). The functional form of the scaling function agrees very well with the predictions shown in Eq. (3).

In Fig. 4 the normalized standard deviation of the scaling plot is shown as a function of $`R`$. Again, the high quality of the scaling can be seen. Note that the statistics for both P-AP and F-DW cases are comparable, while the statistics for $`R=64`$ are worse than for the other values due to the smaller number of samples. From Fig. 4 we see that, for small values of $`R`$, the F-DW case scales less well than the P-AP case. This may explain, why, for F-DW, the stiffness exponent obtained for small aspect-ratios differs significantly from the $`R\mathrm{}`$ limit. The fact that F-DW scales less well than P-AP is probably due to the larger influence of the system boundaries induced by the free bc. As a result, the scaling near $`R=1`$ is more complicated than a simple $`M^{1+\theta }`$ behavior.

Finally, we turn to the system with bimodal distribution of the bonds. In Ref. stiff2d,  a saturation of the domain-wall energy has been observed for $`R=1`$ when $`L=M\mathrm{}`$, i.e. a fundamentally different result from the Gaussian case. This raises the question of whether the aspect-ratio approach is also able to make those results consistent.

As an example, in Fig. 5 the domain-wall energy is plotted for $`M=4`$ as a function of $`L`$. Clearly an exponential decrease with $`L`$ can be observed. This is due to the discrete nature of the interactions and can be explained as follows. For $`LM`$, the system can be decomposed into almost independent subsystems of size $`M\times M`$. Let $`p`$ be the probability to find a non-zero domain-wall energy in a subsystem. Then the probability to find a non-zero $`E_{\mathrm{DW}}`$ in the full system is $`p^{L/M}`$ which decreases exponentially with $`L`$.

As a result, for large $`R`$ (which implies a large value of $`L`$ for fixed $`M`$), the probability to find a non-zero domain-wall energy is so small that for a reasonable number of samples $`E_{\mathrm{DW}}`$ will be exactly zero. This is what we have indeed found for $`R=32`$ and $`R=64`$. It means that for the bimodal model the aspect-ratio approach cannot be applied. This is in contrast to the Gaussian case, where aspect ratio scaling works better and better with increasing $`R`$.

## IV Summary and Discussion

We have studied the scaling behavior of the energy of domain walls in two-dimensional spin glasses obtained by changing the boundary conditions at $`T=0`$. For the numerical calculations, we have employed a very efficient matching algorithm, enabling us to calculate exact ground states of systems up to size $`N=1024\times 16`$. To eliminate finite-size effects, we have applied the aspect-ratio scaling approach. Our main results are:

1. The value of the stiffness exponent for Gaussian distribution of the bonds in $`d=2`$ seems to be independent of the choice of the boundary conditions and the way the domain walls are introduced. We give a final value of $`\theta =0.287(4)`$.
2. The scaling behavior predicted by the aspect-ratio approach applies very well for the Gaussian distribution of the bonds, again indicating that the value of $`\theta `$ is independent of the boundary conditions.
3. The aspect-ratio technique does not allow us to gain further insight for systems with the bimodal distribution of the bonds.

## V Acknowledgments

The simulations were performed on a Beowulf Cluster at the Institut für Theoretische Physik of the Universität Magdeburg and at the Paderborn Center for Parallel Computing, both in Germany. The main part of the work was done during a visit of AKH to Manchester. Financial support for the visit was provided by the ESF programme SPHINX. AKH also obtained financial support from the DFG (Deutsche Forschungsgemeinschaft) under grants Ha 3169/1-1 and Zi 209/6-1. ACC acknowledges support from EPSRC (UK). APY acknowledges support from the NSF through grant DMR 0086287 and the EPSRC under grant GR/R37869/01. He would also like to thank David Sherrington for hospitality during his stay at Oxford.


format: text

text:
# Heat Conduction and Magnetic Phase Behavior in Electron-Doped Ca1-xLaxMnO3 (0≤𝑥≤0.2)

## Abstract

Measurements of thermal conductivity ($`\kappa `$) vs temperature are reported for a series of Ca1-xLaxMnO3($`0x0.2`$) specimens. For the undoped ($`x=0`$), G-type antiferromagnetic compound a large enhancement of $`\kappa `$ below the Neél temperature ($`T_N125K`$) indicates a strong coupling of heat-carrying phonons to the spin system. This enhancement exhibits a nonmonotonic behavior with increasing $`x`$ and correlates remarkably well with the small ferromagnetic component of the magnetization reported previously \[Neumeier and Cohn, Phys. Rev. B 61 14319 (2000).\] Magnetoelastic polaron formation appears to underly the behavior of $`\kappa `$ and the magnetization at $`x0.02`$.

Electronic phase separation has emerged as a paradigm for describing the ground state of strongly correlated electron systems.Dagotto  It may underly the phenomenon of colossal magnetoresistance (CMR) in hole-doped (Mn3+-rich) perovskite manganites studied extensively in recent years. A few studies of electron-doped (Mn4+-rich) manganites Earlye-doped  revealed anomalous magnetic properties and drew attention to these compounds. Subsequent investigationsNeumeierCohn ; Martin ; Savosta ; Mahendiran ; Santhosh ; Respaud ; Aliaga  suggest electron-doped manganites to be novel systems for studies of phase separation and polaron physics.

Undoped CaMnO3 exhibits G-type antiferromagnetic (AFM) order below $`T_N`$125 K. With $`20\%`$ trivalent substitution for divalent Ca, the system adopts a C-type, orbitally-ordered AFM ground state with $`T_N150200`$K, depending on the dopant ions. For Mn3+ concentrations between these end points a small ferromagnetic (FM) moment is observed with $`T_C=T_N(G)`$ and a nonmonotonic doping behavior; data from Ref. NeumeierCohn,  on Ca1-xLaxMnO3 ($`0x0.2`$) are shown in the inset of Fig. 1.

A fundamental issue is whether this FM component reflects a homogeneous canting of AFM moments for all $`x`$ in this regime, as originally proposed by deGennes,deGennes  or whether some portion of the phase behavior can be attributed to FM droplets or polarons as found for the hole-doped compounds.Hennion  A very recent neutron scattering studyLingPreprint  of Ca1-xLaxMnO3 constrains models for the regime $`0.06x0.16`$: (1) the FM moment is perpendicular to the G-AFM moment, consistent with homogeneous canting, and (2) FM clusters with size smaller than $`800`$Å do not exist. Similar experiments at low doping $`x0.05`$ have not been reported, but this regime has been investigated theoretically.Batista ; ChenAllen

The present study of thermal conductivity ($`\kappa `$) vs temperature on Ca1-xLaxMnO3 ($`0x0.2`$) was motivated by the novel phase behavior of electron-doped manganites and by prior workCohnKmang ; CohnRev  demonstrating that the lattice thermal resistivity of manganites is a sensitive measure of bond disorder arising from distorted MnO6 octahedra. Electron hopping via double exchange couples the spins to these octahedral distortions. For the present system the thermal conductivity clearly reflects the lattice response to the FM moment throughout the doping range. At low doping ($`x0.02`$) the evolution of $`\kappa `$ and the magnetization suggests a competition between long-range antiferromagnetism and magnetoelastic polaron formation.

Ca1-xLaxMnO3 polycrystals were prepared by standard solid-state reaction; the preparation methods along with magnetization and resistivity measurements are reported elsewhere.NeumeierCohn  Iodometric titration indicated the oxygen content of all specimens fell within the range 3.00$`\pm 0.01`$. The thermal conductivity was measured in a radiation-shielded vacuum probe using a differential chromel/constantan thermocouple and steady-state technique. Typical specimen dimensions were $`1\times 1\times 3`$mm3. Heat losses via radiation and conduction through leads were measured in separate experiments and the data corrected accordingly. This correction was typically 10-15 % near room temperature and $`2`$% for $`T150`$ K. The specimens have a density of 78$`\pm 4\%`$ that of fully dense material, with no systematic dependence on doping; no porosity corrections have been applied.

Fig. 1 shows $`\kappa (T)`$ for a series of specimens. The data for CaMnO3 ($`x=0`$) indicate a large enhancement of $`\kappa `$ at temperatures below $`T_N(G)=130`$K (for the remainder of the paper we write $`T_N`$ to mean $`T_N(G)`$). The other compounds also show an enhancement, but with a diminished magnitude. For all values of $`x`$ the electronic contribution to $`\kappa `$ at $`T<150`$ K, as inferred from the electrical resistivityNeumeierCohn  and Wiedemann-Franz relation, is less than 5% of the measured value. Furthermore, no substantial changes were observed in resistivities through $`T_N`$. Thus the enhancement is associated with either a magnon or phonon contribution to heat conduction.

To characterize the enhancement, we define the dimensionless change in slope, evaluated just below

$`T_N`$, as $`\mathrm{\Gamma }d(\kappa /\kappa _p)/dt|_{t1}`$, where $`\kappa _p`$ is the $`T`$-dependent thermal conductivity in the paramagnetic state ($`T>T_N`$) and $`t=T/T_N`$. The behavior of $`\kappa _p`$ at $`T<T_N`$ is taken as the extrapolation of polynomial fits to data at $`T>T_N`$ (solid curves in Fig. 1). The doping behavior $`\mathrm{\Gamma }(x)`$ is shown in Fig. 2; of particular interest is the nonmonotonic behavior. $`\mathrm{\Gamma }(x)`$ appears to be composed of two contributions: a term strongly decreasing with $`x`$ and operative for $`x0.02`$, and a term proportional to the FM saturation moment (open triangles and right ordinate, Fig. 2) operative for $`x0.03`$.

Increases in $`\kappa `$ at AF transitions have been observed previously in MnOSlackNewman  and LaMnO3ZhouGoodenough  crystals. The former material undergoes a substantial crystallographic distortion below $`T_N`$ (magnetostriction), and this suggests changes in the lattice heat conduction as a likely mechanism for its $`\kappa `$ enhancement. Lattice anomalies associated with magneto or exchange striction are quite small forRitterNeutron  LaMnO3 andGranadoRaman ; UndopedUltrasound ; NeumeierCornelius  CaMnO3. Nevertheless, prior workCohnKmang ; CohnRev  demonstrated that the lattice thermal resistivity of manganites at $`T300`$K is controlled principally by distortions of the MnO6 octahedra through their influence on phonon scattering rates. The latter can be substantially more sensitive to internal structural modifications than are lattice or elastic constants. Heat conduction by magnons could contribute to the enhancement, but this seems less likely given that magnons contribute negligibly to $`\kappa `$ near $`T_C`$ for ferromagnetic compositions.CohnKmang  The subsequent analysis supports a lattice response to magnetic order as the mechanism for $`\mathrm{\Gamma }`$ in these compounds.

The observation $`\mathrm{\Gamma }M_{sat}`$ for $`x0.03`$ is reminiscent of the behavior found for CMR compounds.CohnKmang  The lattice contribution to $`\kappa `$ in CMR materials is enhanced below the zero-field FM transition temperature and in applied magnetic field at fixed $`T`$ near $`T_C`$. This lattice response correlates with a reduced distortion of the MnO6 octahedra that accompanies double-exchange mediated charge delocalization. Underlying the field and $`T`$-dependent thermal resistivity ($`W=1/\kappa `$) of La0.83Sr0.17MnO3 is a simple magnetization dependence,CohnKmang  $`W(M)W(0)M^2(H,T)`$.

We now show that the nonmonotonic behavior of $`\mathrm{\Gamma }(x)`$ for Ca1-xLaxMnO3 is remarkably well reproduced by a similar phenomenological assumption for $`T<T_N`$: $`W(M,T)W(0,T)M(T)M_p(T)`$. $`M_p(T)`$ is the magnetization of the paramagnetic phase. For zero-field measurements, $`M(T)`$ is the spontaneous magnetization, $`W(M,T)=\kappa ^1`$ is the measured thermal resistivity, and $`W(0,T)\kappa _p^1`$ (the hypothetical thermal resistivity in the absence of magnetic order). This assumption implies that $`\mathrm{\Gamma }`$ should be proportional to the normalized change in slope of the magnetization, evaluated just below $`T_N`$,

$$\mathrm{\Gamma }=\beta \frac{T_N}{M(T_N)}\left(\frac{dM}{dT}|_{TT_N}\frac{dM_p}{dT}|_{TT_N}\right)\beta \mathrm{\Gamma }^M$$
(1)
Using the available $`M(T)`$ curvesNeumeierCohn  measured at $`H=2`$kOe, excellent agreementnoteonEq1  of $`\mathrm{\Gamma }(x)`$ with Eq. (1) was found for all $`x`$, with proportionality constant $`\beta =2.52\times 10^4`$ (open circles and solid curve, Fig. 2).

This result is most easily interpreted in the regime $`x0.03`$ where $`\mathrm{\Gamma }M_{sat}`$. Though no insulator-metal transition takes place at $`T_N`$ akin to that at $`T_C`$ in CMR compounds, there is clear evidence in electrical resistivity measurements near $`x=0.10`$ that some electron delocal-

ization takes place below $`T_N`$; the slopes, $`(1/\rho )d\rho /dT`$, exhibit an abrupt change at the transition (Fig. 3). In analogy with the case of CMR materials, the correlation between $`\mathrm{\Gamma }`$ and magnetization in this regime is plausibly attributed to enhanced electron hopping mediated by double-exchange between aligned (FM droplet scenario) or partially aligned (canting scenario) spins. Enhanced electron transfer reduces the average distortion of the MnO6 octahedra and associated phonon scattering within the FM regions of the specimen. That $`\mathrm{\Gamma }`$ follows both $`M_{sat}`$ and $`\mathrm{\Gamma }^M`$ is consistent with a conventional magnetization of the form $`M=M_{sat}f(T)`$ where $`f(T)`$ reflects the order parameter of the FM phase.

The regime $`x0.02`$ is more complicated since $`M_{sat}`$ has a very different $`x`$ dependence from that of $`\mathrm{\Gamma }^M`$ (and $`\mathrm{\Gamma }`$). These different doping behaviors entail a crossing of the $`M(T)`$ curves for different $`x`$ at $`T<T_N`$ (Fig. 4). The data suggest that two independent components contribute to the magnetization in this regime, one with $`T_C=125`$K, characterizing the undoped specimen, and the other with $`T_C115`$K, characterizing $`x0.04`$. A smooth evolution between the two is reflected in the $`M(T)`$ data for the intervening compositions. This coexistence is most evident in the curve for $`x=0.02`$ as an inflection near $`T=110`$ K.

$`\mathrm{\Gamma }`$ follows the diminution of the higher-$`T`$ transition and presumably has its origin in the same coupling between spins and octahedral distortions underlying the response at $`x0.03`$. Supporting this hypothesis is a recent study of Raman scattering in similar compounds.GranadoRaman  Sharply enhanced Raman intensities at $`T<T_N`$ for low-frequency, rotational and bending modes of the oxygen octahedra were observed for CaMnO3. With increasing La doping this enhancement below $`T_N`$ was diminished, and was absent for $`x`$=0.03, very similar to the trend observed here for $`\mathrm{\Gamma }`$. It seems likely that the two phenomena are related.

Another important experimental result relevant

to a description of the data at low $`x`$ comes from recent magnetic neutron scattering studies of an $`x=0.02`$ compound in magnetic field.GranadoPrivate  The field dependence of the AF scattering intensity was inconsistent with a FM component arising from either uniform spin canting or ferrimagnetism; the FM component is decoupled from the AF background in applied field.

These observations implicate magnetoelastic polarons in the FM and lattice response at low $`x`$. In this regard the behavior of the slope $`dM_{sat}/dx`$ is of interest (inset, Fig. 1). It increases from $`1\mu _B`$/Mn-ion-electron for $`x0.02`$ to $`7\mu _B`$/Mn-ion-electron for $`0.03x0.08`$. The latter is the value expected if each La dopant adds a symmetric 7-site FM polaron,NeumeierCohn  determined in calculationsChenAllen  to be the stable ground state for this system. The value $`x0.03`$ appears to mark a crossover between regimes. The mean spacing between dopant ions is estimated as, $`r_{La}=(3V/16\pi x)^{1/3}`$ ($`V=207\mathrm{\AA }^3`$ is the unit cell volume,LingPreprint  containing 4 f.u.’s). For $`x=0.03`$, $`r_{La}7.4\mathrm{\AA }`$, equal to the 3rd-nearest neighbor Mn distance. This would be the distance between the centers of 7-site polarons sharing a single Mn site (a polaron “cluster”), but other polaron configurations are close in ground-state energy and might be favored if interactions not yet considered in calculationsChenAllen  (e.g. defects, next-nearest neighbor exchange) play a role. For example, larger polarons with canted-spin arrangements are interesting candidates because the large-scale clusters anticipated at higher doping could evolve smoothly into the long-range, spin-canted state proposed on the basis of neutron scatteringLingPreprint  for $`x0.06`$.

A model that distinguishes between isolated and clustered polarons describes the qualitative features of the experimental results. We propose that the $`T_C=125`$ K transition in Fig. 4 is associated with isolated polarons (i.e. those without adjacent polarons) and the lower-$`T_C`$ transition to clusters of two or more polarons. This is plausible since the balance between kinetic, lattice and spin energies, possibly different for a cluster, could yield a polaronic state less robust against thermal fluctuations. Presumably a small (principally isolated) polaron density at $`x=0`$ is associated with native defects (e.g. oxygen or Mn vacancies). With increasing La doping, polaron clusters are produced at the expense of isolated polarons. For a random distribution, the probability for two or more polarons to share a Mn site grows rapidly for 7-site or larger polarons given the large number of Mn sites over which electrons may be distributed to define a cluster. This provides a natural explanation for the rapid decrease of $`\mathrm{\Gamma }`$ and $`\mathrm{\Gamma }^M`$ with $`x`$ at $`x0.02`$, and for their minima at $`x0.03`$ where the FM transition is maximally nonuniform.

What is the origin of the small value for $`dM_{sat}/dx`$ at $`x0.02`$ and its increase near $`x0.03`$? One possibility is that these features reflect differing magnetic moments for isolated and clustered polarons (e.g. different canting angles for spin-canted polarons). Alternatively, these features might reflect a change in the number of polarons induced per dopant, e.g. if perturbations associated with defects (vacancies, La ions) suppress polaron formation at $`x0.02`$. Experimental investigations that better define the polaron characteristics are required to refine these ideas.

In summary, systematic changes in the behavior of the thermal conductivity at the AF transition in La-doped CaMnO3 and their correlation with magnetization measurements indicate that the lattice thermal resistivity is a sensitive probe of FM interactions through the coupling of spins to local distortions of the MnO6 octahedra. This extends similar conclusions obtained previously for hole-doped, CMR compoundsCohnKmang ; CohnRev  to the present system where the ground state appears to consist of magneto-elastic polarons ($`x0.02`$) and a spin-canted phase ($`x0.06`$). The crossover between these two regimes, clearly manifested in both the thermal conductivity and magnetization data, appears to reflect novel polaron physics and is of particular interest for further study.

The authors acknowledge experimental assistance from Dr. B. Zawilski and Dr. R. Maier. The work at the University of Miami was supported, in part, by NSF Grant No.’s DMR-9631236 and DMR-0072276, and at Florida Atlantic University by DMR-9982834.


text:
# 1 Introduction

## 1 Introduction

In  we have shown how it is possible to build a physically self-consistent model for an accretion disc corona: assuming that the turbulent magnetic stresses generated by MRI are responsible for angular momentum transport in the disc, and that the field saturates mainly due to buoyancy (as should be the case if strong coronae are to be generated), then the fraction of power released into the corona as a function of radius, $`f(R)`$, can be uniquely determined by solving the algebraic equations for accretion disc structure. For uniform discs (which may not be the case in the radiation pressure dominated parts, see ), we have

$$f(R)\sqrt{\alpha }\left(1+\frac{P_{\mathrm{rad}}(R)}{P_{\mathrm{gas}}(R)}\right)^{1/4},$$
(1)
where $`\alpha `$ is a constant of the order of unity and $`P_{\mathrm{rad}}`$ and $`P_{\mathrm{gas}}`$ represent radiation and gas pressure in the disc, respectively.

A major consequence of such a disc-corona coupling, which is of interest in the context of microquasars, is that coronae are more powerful at low accretion rates, in particular below the critical rate at which the radiation pressure dominated part of the disc disappears altogether. Indeed, if the $`\alpha `$ viscosity parameter is high enough, $`f`$ can approach unity, and we are therefore left with a flow in which most of the energy goes into the magnetically dominated corona. The cold, geometrically thin disc, although very dim, may manifest itself as a reprocessor of X-ray coronal radiation. This makes our solution different from a simple magnetically dominated non radiative accretion flow (NRAF): the two may be distinguished, for example, by the presence of relativistically smeared reflection features.

## 2 MHD outflows from powerful coronae

A magnetically dominated corona, with magnetic flux tubes corotating with the underlying disc at nearly Keplerian speed, $`v_{\mathrm{Kep}}`$, is prone to generate powerful MHD outflows. Let us examine in some detail the disc-corona-outflow energetics. At a distance $`R`$ from the center, the (magnetic) energy flux emerging from the disc is

$$Q_\mathrm{c}(R)=f(R)Q(R)=\frac{3GM\dot{M}f(R)}{8\pi R^3}\mathrm{ergs}\mathrm{cm}^2\mathrm{s}^1,$$
(2)
Part of this energy will be dissipated in reconnection events, heating the coronal plasma that is then cooled by inverse Compton scattering soft photons. Such energy flux will emerge as hard X-rays, and can be parametrized as $`F_\mathrm{H}(R)=(1\eta )Q_\mathrm{c}(R)`$, where $`\eta `$ is in general a function of the distance. The remaining fraction of the coronal magnetic energy flux is carried off by the jet/outflow: $`Q_\mathrm{j}(R)=\eta Q_\mathrm{c}(R)`$. Most models of MHD generated outflows from rotating systems  agree on the conclusion that the final power carried by such jets depends on the poloidal component of the magnetic field and on the angular velocity of the field lines. We have therefore $`dL_\mathrm{j}=\xi 2\pi RQ_\mathrm{j}(R)dR=(B_\mathrm{p}^2/8\pi )d𝒜_c\mathrm{\Omega }R`$, where $`d𝒜_c\xi 2\pi RdR`$ is the element of disc area covered by the corona (which has a covering fraction $`\xi `$) and the poloidal component of the coronal magnetic field is a function of the corona scaleheight, $`H_c`$ (the height of a reconnection site, from which the outflow is launched): $`B_\mathrm{p}B(H_c/R)`$. The intensity of the coronal magnetic field, in turn, depends on how quickly the field dissipates $`B^2=\frac{8\pi Q_\mathrm{c}(1\eta )}{\xi v_{\mathrm{diss}}}`$.

Combining the expressions for $`dL_\mathrm{j}`$ and $`B`$, we obtain

$$\eta =\left(1+\xi \left(\frac{R}{H_c}\right)^2\frac{v_{\mathrm{Kep}}(R)}{v_{\mathrm{diss}}}\right)^1.$$
(3)
The fraction of coronal power that is channeled in a MHD jet/outflow, depends crucially on the height of a reconnection site (or, more generally, on the scaleheight of the magnetically dominated corona).

## 3 Radial profiles

Most of the variability of black hole binaries (both the broad band one and the QPOs) can be spectrally associated to the hard component . Moreover, QPO frequencies tend to correlate with various X-ray spectral parameters; in particular with those sensitive to the system size. In general, higher frequencies are associated to softer spectra. In principle, if a magnetic corona is responsible for generating (or amplifying) the observed modulation of the hard X-ray flux, and the corona is vertically extended (large $`H_c`$), both radial and vertical oscillation may be observed. This implies that both the radial and vertical size of the corona play an important role in shaping the Power Density Spectra of these sources . For example, if jet and corona are coupled systems, as we envisage, QPO properties (in particular their amplitude and narrowness) are bound to be influenced by the jet. Indeed,  shows that when the microquasar GRS 1915+105 is in a hard state characterized by strong radio emission (states $`\chi _1`$ and $`\chi _3`$ according to the classification of ), the high frequency noise component in the PDS is strongly suppressed as compared to the hard states that do not show radio emission (states $`\chi _2`$).

Radial profiles of coronal emissivity and of the fraction of the total gravitational power channeled either into coronal heating or into an outflow, should be derived in order to build more physical models for the QPO mechanism. As a final goal, time dependent versions of such models should be implemented. Here we present a simple, illustrative example of a possible calculation. We have fixed the free parameters $`\xi =0.5`$ and $`c/v_{\mathrm{diss}}=30`$. The coronal optical depth is fixed at $`\tau =0.8`$ (variation of this parameters affect the slope of the X-ray spectrum but not the overall energetics), and the temperature is calculated self-consistently from the heating-cooling balance (including both coronal synchrotron and disc radiation fields as sources of soft photons for Comptonization) as in , at every annulus of width $`dR`$. The most critical assumption is that regarding the radial profile of the coronal scaleheight $`H_c`$, for the influence it has in determining the fraction of power that goes into the jet. In general, the coronal scaleheight will be larger for a larger value of $`f`$ (because more strongly buoyant flux tubes should rise more) and proportional to the number of twists a coronal loop, placed at a distance $`R`$ and with footpoints separated by a distance $`dR`$, experience due to the disc differential rotation that torques the flux tube. Also, it should have a power-law dependence on the distance due to the natural disc flaring. In summary, we may write $`H_c/RR^pf(R)d\mathrm{ln}\mathrm{\Omega }/dR=Af(R)R^{p1}`$, with $`0<p<1`$. In our calculations, we chose $`p=1/2`$ and $`A=1`$.

Once the accretion rate is fixed, we can integrate the disc structure equations with the prescription of Eq. (1) for the coronal fraction, and solve for $`f(R)`$ and $`\eta (R)`$. With the solution in hand, we can calculate the intensity of the magnetic field in the corona as a function of $`R`$ and then compute the amount of energy flux that is converted locally into hard X-rays, into soft quasi thermal emission (including both intrinsic and reprocessed disc emission), into self-absorbed synchrotron radiation and into jet-power. In Figure 1 we show the solutions for a 10 $`M_{}`$ black hole accreting at $`\dot{m}=0.005`$ (typical low/hard state accretion rate, thick lines) and at $`\dot{m}=0.25`$ (thin lines).

At low accretion rates, the hard X-rays from the corona dominate over the quasi-thermal disc radiation. In the innermost region of the corona the energy density of the synchrotron radiation field dominates over the external disc one, and coronal heating produces strong, highly variable IR/Optical non-thermal emission. The kinetic power channeled into the jet is the main repository of the gravitational energy of the accretion flow in the inner 10 Schwarzschild radii. At high accretion rates the inner part of the disc is radiation pressure dominated, the corona is strongly suppressed, and no strong outflow is generated, as expected.

## 4 Conclusions

If a standard, geometrically thin and optically thick accretion disc is coupled to a corona through buoyancy of magnetic flux generated by MRI in the disc itself, then a powerful corona is generated whenever gas pressure dominates over radiation pressure, i.e. at low accretion rates. This implies that microquasar in the low/hard state may be characterized by strong coronae on top of very dim, cold discs, that may be seen only as reprocessors. Depending on the detailed coronal geometry, and in particular on its vertical extent, low accretion rate systems may produce strong MHD driven jets/outflows whose total kinetic energy flux may exceed the radiated one. Beside their radio emission, such outflows should manifest themselves through their influence on the variability properties (noise and QPOs), for example reducing the intensity and narrowness of such timing features. In order to understand the complex interplay between coronal physics and dynamics and the variability properties of microquasars, full time- and radial-dependent models are needed. As a first illustrative step in this direction, we have shown here a stationary model for the radial profiles of the different spectral components emerging from a disc-corona system at various accretion rates.


text:
# Recent results on spectra and yields from NA49

## 1 Introduction

In the last few years NA49 has taken data at 40, 80 and 158 $`A`$GeV beam energy to study the energy dependence of particle production in nucleus-nucleus collisions and to search for indications of the onset of deconfinement in these collisions. In this paper recent results on particle spectra and yields are reported. Results by NA49 on other observables such as HBT, fluctuations and system size dependences from p+p, p+A, and A+A (peripheral Pb+Pb, central Si+Si and C+C) are presented elsewhere in these proceedings .

In addition to the energy scan programme, a large dataset of central Pb+Pb data has been collected at 158 $`A`$GeV, allowing to study rare particles such as $`\mathrm{\Omega }`$, $`D`$, and $`\varphi \mathrm{e}^+\mathrm{e}^{}`$. Results on $`\mathrm{\Omega }`$ and open charm production are presented here.

## 2 Experiment and centrality selection

The NA49 experimental setup  consists of four TPCs, two of which are in a magnetic field. The Main TPCs (MTPC) are outside the magnetic field and are used to identify particles by measuring the energy loss $`dE/dx`$. The particle identification capabilities of the MTPCs are augmented by two time-of-flight (TOF) detector arrays which have acceptance around midrapidity for kaons and close to midrapidity for protons and deuterons. Centrality selection of the events is done using a calorimeter which detects the projectile spectators. At 40 and 80 $`A`$GeV the 7.2% most central events were selected. For the 158 $`A`$GeV data the on-line centrality cut is at 10%, but for the analysis of pions, kaons and $`\varphi `$ an off-line cut was made at 5%. The high statistics data sample used for the $`\mathrm{\Omega }`$ analysis at 158 $`A`$GeV was taken with a centrality trigger at 20% of the inelastic cross section.

## 3 Transverse mass spectra and thermal freeze-out

In Fig. 1 all available results on $`m_T`$-spectra as measured in central Pb+Pb collisions by NA49 are collected.222The $`\overline{p}`$, $`\mathrm{\Omega }^{}`$ and $`\overline{\mathrm{\Omega }}^+`$ $`m_T`$-spectra at 158 $`A`$GeV appearing in the printed version of these proceedings (Nucl. Phys. A715 161c) are unfortunately incorrect and therefore differ slightly from those presented here. The spectra are measured at or close to midrapidity except for the $`\varphi `$, $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ which are integrated over approximately 1 unit of rapidity, to increase statistics.

The $`\mathrm{\Xi }`$ and $`\varphi `$ at 158 $`A`$GeV have already been published  while the kaon and pion spectra are submitted for publication . All other presented results are still preliminary.

In a hydro-dynamical picture, the $`m_T`$-spectra of particles are sensitive to transverse flow. To characterise the flow, the spectra were fitted with

$$\frac{dN}{m_Tdm_Tdy}m_TK_1\left(\frac{m_T\mathrm{cosh}\rho }{T}\right)I_0\left(\frac{p_T\mathrm{sinh}\rho }{T}\right).$$
(1)
A combined fit of several spectra with this function allows to determine the thermal freeze-out temperature $`T`$ and the mean transverse flow velocity $`\beta _T`$ ($`\rho =\mathrm{atanh}\beta _T`$). The $`\pi ^{}`$ and deuteron spectra were excluded from the fit (dashed lines in Fig. 1) because the pions are expected to have a significant contribution from resonance decays and the deuterons may be formed by coalescence. The particles (baryons and $`K^+`$) and anti-particles (anti-baryons and $`K^{}`$) were fitted separately. The $`\varphi `$ was included in the anti-particle fit, because it is not expected to be sensitive to the baryon density. The fit describes the data reasonably well, even including the pions at higher $`m_T`$ and the deuterons. The relatively large $`\chi ^2`$ per degree of freedom can partly be explained by the fact that only statistical errors were taken into account in the fit (systematic errors are not yet available for all spectra). Since the slopes of the $`m_T`$-spectra of particles and anti-particles are similar, the differences between the freeze-out parameters for both groups of particles are small.

The differences between the three different energies are not larger than the difference between the particles and anti-particles at each energy, implying that the energy dependence of the thermal freeze-out conditions is small at these energies.

The $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ spectra, which are only available at 158 $`A`$GeV, do not show large deviations from the fit. Thus, based on these data, there is no indication of an early freeze-out as has been suggested by several authors . The slope parameters of the $`\mathrm{\Omega }`$ spectra are compatible with earlier measurements by WA97 .

## 4 Rapidity spectra and baryon density

The $`m_T`$-spectra of kaons as shown in Fig. 1 were obtained using combined $`dE/dx`$ and TOF information at midrapidity. The $`dE/dx`$ measurement in the TPC alone allows to identify kaons at forward rapidities. The resulting spectra are shown in Fig. 2, together with $`\pi ^{}`$ distributions which were obtained using unidentified negatively charged particles, corrected for contributions of $`K^{},\overline{p}`$ and non-vertex tracks.

The $`K^{}`$ distributions are narrower than the $`K^+`$ distributions at all energies. This could indicate that the kaons are sensitive to the local baryon density, which is reflected in the difference between the $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ spectra as well (for details on the $`\mathrm{\Lambda }`$ analysis, see  in these proceedings).

The $`\varphi `$ spectra are obtained from invariant mass distributions of $`K^+K^{}`$ pairs. The resulting rapidity spectra are shown in Fig. 3, where the new results at 40 and 80 $`A`$GeV are compared to the published result at 158 $`A`$GeV . Clearly, the width of the distribution, as well as the total yield, increases with beam energy.

The large statistics data sample of the 20% most central Pb+Pb collisions at 158 $`A`$GeV beam energy, which was taken in the year 2000, has been used for the $`\mathrm{\Omega }`$ analysis. The resulting $`m_T`$-spectra are shown in Fig. 1 and the corresponding rapidity spectra in Fig. 5. Integrating over rapidity gives mean multiplicities per event of $`\mathrm{\Omega }^{}=0.47\pm 0.07`$ and $`\overline{\mathrm{\Omega }}^+=0.15\pm 0.02`$. The observation that the ratio $`\overline{\mathrm{\Omega }}^+/\mathrm{\Omega }^{}`$ is below unity cannot be explained in a pure string fragmentation model but fits to expectations in an equilibrium hadron gas model . Also note that the $`\mathrm{\Omega }^{}`$ has a slightly wider rapidity distribution than the $`\overline{\mathrm{\Omega }}^+`$. This indicates that even the $`\mathrm{\Omega }`$ is sensitive to the effect of baryon density.

## 5 Total yields and chemical freeze-out

The measured rapidity spectra have been integrated to obtain the mean multiplicities in full phase space. To determine the chemical freeze-out parameters at the three different energies, the total yields were fitted with a hadron gas model with partial strangeness saturation . The obtained temperatures for chemical freeze-out (145–160 MeV) are larger than the thermal freeze-out temperatures (115–140 MeV). The chemical freeze-out temperature $`T`$ and baryon chemical potential $`\mu _B`$ are compared to results at AGS and RHIC in Fig. 5. The freeze-out points follow a smooth curve of increasing temperature and decreasing baryon density with energy. The observed freeze-out parameters are close to the calculated phase boundary  at SPS energies and above.

## 6 Energy dependence of strangeness production

The mean multiplicities of the most abundant carriers of strangeness, kaons and lambdas, are plotted as a function of the centre of mass energy in Fig. 6. Only the charged kaons are shown because the neutral kaons are not measured at all energies. All yields are divided by the total pion yield $`\pi =1.5(\pi ^++\pi ^{})`$. It is seen that the $`K^+/\pi `$ ratio peaks at 40 $`A`$GeV, whereas $`K^{}/\pi `$ increases monotonically. On the other hand, the $`\mathrm{\Lambda }`$ production also peaks around 40 $`A`$GeV, whereas the $`\overline{\mathrm{\Lambda }}`$ yield increases monotonically with energy. Apparently, both the total $`s`$ production and the $`\overline{s}`$ production have a maximum close to 40 $`A`$GeV. The maximum in the $`\overline{s}`$ production is directly seen $`K^+/\pi `$ ratio because the contribution of $`\overline{\mathrm{\Lambda }}`$ is small. Although the $`K^{}`$ represent a large fraction of the $`s`$ production, the maximum in the $`s`$ production is reflected in the $`\mathrm{\Lambda }/\pi `$ ratio.

In Fig. 8 the energy dependence of the $`K^+/\pi ^+`$ ratio is compared to several models which do not assume the formation of a QGP in the collision. The UrQMD model  fails to describe the data because the pion yields are over-estimated by about 30%. The RQMD model  gives a much better descriptions, but over-estimates the ratios at 80 and 158 $`A`$GeV. The same is true for the hadron gas model . In this model, the $`K^+/\pi ^+`$ ratio peaks around 40 $`A`$GeV due to the interplay of the decreasing baryon density and the increasing temperature with energy.

In Fig. 8 the energy dependence ($`F(\sqrt{s}2m_N)^{3/4}/\sqrt{s}^{1/4}`$) of the strangeness to entropy ratio $`E_S=(\mathrm{\Lambda }+K+\overline{K})/\pi `$ is compared to a model  which assumes that above a certain energy a QGP is formed in the collision. The peak in the model curve is caused by this transition. In  it is shown that the model also describes the energy dependence of pion production.

NA49 intends to complete the energy scan by taking data at 20 and 30 $`A`$GeV in 2002.

## 7 Open charm upper limit

All available central Pb+Pb data at 158 $`A`$GeV (900k events at 10% and 3M events at 20% centrality) were used in an invariant mass analysis to measure open charm. The invariant mass was calculated for pairs of pion and kaon candidates, to detect the decays

$$D^0K^{}+\pi ^+\mathrm{and}\overline{D}^0K^++\pi ^{},$$

which have a branching fraction of 3.83%. For tracks which reach the MTPC, the kaon candidates were selected by $`dE/dx`$ measured in the MTPC. For tracks which do not reach the MTPC, all candidates were used. The invariant mass histogram for $`D^0`$ before background subtraction is shown in the left panel of Fig. 9. The result after adding the $`\overline{D}^0`$, subtracting the background and correcting for acceptance and efficiency losses is shown in the right panel of Fig. 9. No signal of the decays is visible. Detector simulations have shown that the shape of the expected signal can be parametrised as a Breit-Wigner with a width $`\mathrm{\Gamma }=11`$ MeV. Using this parametrisation, predictions by different models are shown in the figure. It is seen that thermal production of charm quarks in a Quark Gluon Plasma at $`T=265`$ MeV, as described in , is excluded (dotted curve in Fig. 9). The expectation from the quark coalescence model ALCOR  is at the limit of the sensitivity of the present analysis (dashed curve). Lower values of the charm yield, such as the expectation from perturbative QCD, or a factor 3 enhancement  (full curve) cannot be excluded. In order to quantify the result, the data have been fitted with the parametrised peak shape, leaving the normalisation as the only free parameter. The result is $`D^0+\overline{D}^0=1.07\pm 0.84`$.

## 8 Summary and outlook

Spectra of $`\pi ,K,p,d,\varphi `$, and $`\mathrm{\Lambda }`$ in central Pb+Pb collisions at 40, 80 and 158 $`A`$GeV are presented, together with new results on $`\mathrm{\Omega }`$ rapidity and $`m_T`$ spectra at 158 $`A`$GeV.

The $`m_T`$-spectra are compatible with radial flow, with similar temperature $`T130`$ MeV and transverse flow velocity $`\beta _T0.45`$ at the different energies. The $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ spectra, which are only measured at 158 $`A`$GeV, do not show indications of early freeze-out.

The chemical freeze-out parameters as determined from the measured total yields approach the phase boundary in the SPS energy range.

The energy dependence of strangeness production shows a maximum in the relative strangeness yield at 40 $`A`$GeV. This maximum is observed in the $`K^+/\pi `$ ratio for the $`\overline{s}`$ and in the $`\mathrm{\Lambda }/\pi `$ ratio for the $`s`$ quarks. The maximum can be explained by the onset of deconfinement in the energy range between AGS and SPS. Upcoming runs at 20 and 30 $`A`$GeV will provide a more accurate localisation of the maximum.

To investigate open charm production in heavy ion collisions, an invariant mass analysis has been performed using a large sample of 158 $`A`$GeV Pb+Pb events. No signal has been observed. The result excludes equilibrium charm production in a QGP of $`T=265`$ MeV.

### Acknowledgements

We would like to thank F. Becattini for providing the results of his hadron gas fits.

This work was supported by the Director, Office of Energy Research, Division of Nuclear Physics of the Office of High Energy and Nuclear Physics of the US Department of Energy (DE-ACO3-76SFOOO98 and DE-FG02-91ER40609), the US National Science Foundation, the Bundesministerium fur Bildung und Forschung, Germany, the Alexander von Humboldt Foundation, the UK Engineering and Physical Sciences Research Council, the Polish State Committee for Scientific Research (5 P03B 130 23 and 2 P03B 02418), the Hungarian Scientific Research Foundation (T14920 and T32293), Hungarian National Science Foundation, OTKA, (F034707), the EC Marie Curie Foundation, and the Polish-German Foundation.


text:
# 1 Results

## 1 Results

Greiner et al.  identified a K-M III giant as the mass-donating component of GRS 1915+105 in a 33.5 day orbit around a 14 $`M_{}`$ black hole (see Fig. 1). We assume that the donor fills its Roche lobe and that the mass loss is determined by evolutionary expansion along the giant branch, conserving total mass and angular momentum. The growth of the He-core, resulting in an increase of the radius, is determined by the luminosity due to hydrogen shell burning which depends completely on the core mass .

Assuming the donor mass is around 1.2 $`M_{}`$ it is rather easy to estimate the mass transfer from the donor (1 – 2) $`\times `$ 10-8 M yr-1 (Fig. 2 and Table 1, for details see ).

The mean X-ray luminosity over the past 10 years (1.2$`\times `$1039 erg s-1 using a distance of 12.5 kpc, see Fig. 3) gives an accretion rate of (0.2/$`\eta `$) 10-7 M yr-1 onto the primary black hole where $`\eta `$ is the conversion factor of mass infall into radiation (L = $`\eta `$ $`\dot{M}`$ c2). Finally, the ratio (mass transfer/accretion) gives the duty cycle (0.05 - 0.1)($`\eta `$/0.1), the active time as a fraction of the total life time.

We suggest that the quiescent OFF-time is the viscosity time scale at the circularization radius (15 solar radii from the BH) where the local Keplerian angular momentum equals to that in the L1-point. For small $`\alpha `$ (= 0.001) the time scale equals to 300 – 400 years. In this particular case the present active phase may last 10 – 20 years more.

The whole microquasar phase may have started when the secondary made its first contact with its Roche lobe just after finishing the main sequence phase or while crossing the Hertzsprung gap, in a 10-day binary. The final end product after 107 years will probably be a long-period detached white dwarf – black hole system.

## Acknowledgments

.

This research is part of the consortium ‘High Energy Astrophysics and Space Astronomy’ financed by the Academy of Finland and the Finnish Technology Agency TEKES within the framework of the Space Research program ANTARES.


text:
# Dipolar origin of the gas-liquid coexistence of the hard-core 1:1 electrolyte model

## I Introduction

In recent years there has been an increasing interest in phase transitions between fluid phases of different electrolyte concentrations in ionic solutions. Two different regimes have been identified experimentally Fisher-review . The “solvophobic” regime occurs for large solvent dielectric constants that effectively turn off Coulombic interactions. Consequently, solvophobic phase transitions are primarily driven by unfavourable interactions between solute and solvent. This behavior is well described by the usual theory of nonelectrolyte solutions with short range interactions, and clearly leads to Ising-like critical behavior. By contrast, in the “Coulombic” regime the solvent has a low dielectric constant and electrostatic interactions between the solute ions drive the phase separation. In this case, the phase diagrams are quite asymmetric, and apparently mean-field critical behavior has been claimed Singh ; Weingatner ; Zhang , although Ising-mean field crossover has also been seen within a narrower range of temperatures around the critical than in non-ionic fluids Narayanan ; Jacob ; Anisimov ; Gutkowski . Some of these experimental studies suggest that there exists a new characteristic length in these systems that competes with the correlation length for density fluctuations Anisimov ; Gutkowski .

Electrolyte systems in the Coulomb regime are often modelled as charged hard spheres embedded in a uniform dielectric continuum (primitive models). Most studied is the “restricted primitive model” (RPM), in which the ions are of equal size. A vapor-liquid phase transition at very low temperatures and densities was predicted theoretically 25 years ago Stell1  and by early computer simulation studies Vorontsov . Improvement of computer simulation techniques has allowed an increasingly precise determination of the coexistence parameters Orkoulas ; Caillol ; DePablo1 ; Thanos ; Thanos4 ; Luijten . There have also been recent studies of primitive models with asymmetry in size E2PM ; DePablo2 ; Thanos2  and charge Thanos2 ; DePablo3 . Very recent results suggest that the critical behavior of the RPM belongs to the Ising universality class Luijten .

From a theoretical point of view, different approaches have been used in order to explain the vapor-liquid transition of the primitive models. Integral equations such as the mean spherical approximation (MSA) Levin ; Tovar , as well as the Debye-Hückel theory Levin  and Poisson-Boltzmann approaches Outhwaite  have succesfully been applied to it. For the RPM, the most succesful theories are the pairing theories, that consider the ionic fluid as a mixture of bound pairs and free ions in chemical equilibrium Levin ; Raineri , in the spirit of the Bjerrum’s ideas Bjerrum . However, in all these theories the transition is driven by the free ions, even when in same cases, as in the Debye-Hückel-Bjerrum approach, the associated pairs are the dominant species. Analytical Gillan  and computer simulation Bresme  results, on the other hand, show that the structure of the vapor phase is dominated by neutral clusters, mostly dimers and tetramers. Computer simulations have also demonstrated that the phase envelope of the RPM resembles that of the charged hard dumbbell model Shelley ; E2PM . This result has been confirmed by theoretical studies Kalyuzhnyi ; Jiang . The question that the present work examines in detail is the influence of ionic association on the vapor-liquid transition of primitive model electrolytes. In contrast to an earlier study E2PM , which only considered tightly bound dimers, here we examine a range of models with varying values of $`R_c`$, the maximum separation between ions in a dimer. $`R_c\mathrm{}`$ corresponds to the Stillinger-Lovett description Stillinger  of the free ion system.

The structure of this paper is as follows. In Section II, we examine the microscopic structure of the coexisting phases of the ionic fluid and compare the correlation functions to those of a tightly tethered dimer model. In Section III, an exact chemical representation of the ionic fluid as a mixture of associated pairs and free ions is introduced. The role of the associated pairs on the phase coexistence is studied in detail in Section IV and the paper closes with discussion and conclusions.

## II Microscopic structure of the ionic and tethered dimer systems at coexistence

In this Section we analyse the role of pairing on gas-liquid coexistence of primitive model ionic fluids. We have studied by computer simulation an 1:1 size-asymmetric primitive model, in which the ions are modelled as hard spheres of diameters $`\sigma _+`$ and $`\sigma _{}`$, and carrying charges $`+q`$ and $`q`$, respectively, embedded on a dielectric continuum of dielectric constant $`D`$ ($`D=1`$ for the vacuum). The interaction potential between two ions separated by a distance $`r_{ij}`$ is given by:

$$u_{ij}(r_{ij})=u_{hs}(r_{ij})+\frac{q_iq_j}{Dr_{ij}}$$
(1)
where $`u_{hs}(r_{ij})`$ is the hard-core potential that takes the value $`+\mathrm{}`$ if $`r_{ij}<\frac{1}{2}\left(\sigma _i+\sigma _j\right)`$ and 0 otherwise.

The size asymmetry is characterized by the parameter $`\lambda `$, defined as:

$$\lambda =\frac{\sigma _{}}{\sigma _+}$$
(2)
Monte Carlo simulations in the neutral grand-canonical ensemble were performed, characterized by a temperature $`T`$ and the *configurational* chemical potential for a pair of unlike ions $`\mu =2\mu _{real}3k_BT\mathrm{ln}\left(\mathrm{\Lambda }_+\mathrm{\Lambda }_{}\right)`$, with $`\mu _{real}(\mu _++\mu _{})/2`$ and $`\mathrm{\Lambda }_\pm `$ are the thermal de Broglie wavelengths of the ionic species:

$$\mathrm{\Lambda }_\pm =\sqrt{\frac{h^2}{2\pi m_\pm k_BT}}$$
(3)
As usual, cubic boxes of length $`L`$ under periodic boundary conditions are used. The long-ranged character of the Coulombic interactions is handled by the use the Ewald summation tecnique with conducting boundary conditions, with 518 Fourier-space wavevectors and real-space damping parameter $`\kappa =5`$. The relative error due to the infinite sums truncation in the electrostatic energy is less than $`10^5`$ for random configurations in small systems Hummer , and this choice has been also validated by direct simulations of the RPM Thanos4 . In order to speed up the simulations, the basic steps of the simulations (insertion and deletions of pairs of unlike ions) are biased following Orkoulas . Moreover, a fine-discretization approximation is used: the positions available to each ion are the sites of a simple cubic lattice of spacing $`a`$. This methodology has succesfully been applied to the RPM Thanos , to 1:1 size-asymetric primitive models E2PM  and to $`z`$:1 size-asymetric primitive models Thanos , and allows a speedup relative to the continuum calculations of a factor of 100 for small systems. The results are almost indistinguishable from the continuum ones for a discretization parameter $`\zeta \sigma _\pm /a=10`$ Thanos ; E2PM ; Thanos2 , where $`\sigma _\pm =\frac{1}{2}\left(\sigma _++\sigma _{}\right)`$ is the unlike-ion collision diameter. This value of the discretization parameter ($`\zeta =10`$) was used in the present study. Histogram reweighting techniques Ferrenberg  and mixed-field finite size scaling methods Bruce  were used to obtain the vapor-liquid envelopes and the effective critical points, respectively. For the sake of comparison, we have also studied tightly tethered dipolar dimer systems Shelley ; E2PM , consisting of $`N`$ pairs of a positive and negative ion restricted to remain at separations $`\sigma _d`$ satisfying $`\sigma _\pm \sigma _d1.02\sigma _\pm `$. Simulation details and some preliminary results for a range of values of $`\lambda `$ are presented in Ref. E2PM .

The unlike-ion collision diameter $`\sigma _\pm `$ provides the basic length scale appropriate for defining both the reduced temperature and reduced density via

$$T^{}=k_BTD\sigma _\pm /q^2\mathrm{and}\rho ^{}=2N\sigma _\pm ^3/L^3$$
(4)
where $`N_+=N_{}N`$ is the particle number of each ionic species Tovar ; Orkoulas ; E2PM . The reduced simulation box length is defined similarly via $`L^{}=L/\sigma _\pm `$, and the reduced energies and chemical potential as $`U^{}=UD\sigma _\pm /q^2`$ and $`\mu ^{}=\mu D\sigma _\pm /q^2`$.

The value of the asymmetry parameter considered in this paper is $`\lambda =3`$. For that case, the gas-liquid coexistence of the ionic fluid shows a shift in both temperature and density respect to the tethered dimer fluid (see Fig. 1), which is a general feature when comparing ionic and tethered dimer systems E2PM . On the other hand the asymmetry is not high so as to favor large chain-like neutral clusters, as occurs for bigger values of $`\lambda `$ E2PM ; DePablo2 . These features qualify this case to be a typical example for moderately asymmetric 1:1 electrolytes, including the RPM.

Fig. 1 makes the similarity between the phase diagram of the ionic and the tethered dimer fluids clear. It also suggests in an indirect way that the pairing plays a decisive role on the gas-liquid coexistence in the ionic fluid. In order to clarify such a role, we have studied the ion-ion radial distribution functions $`g_{ij}(r)`$ corresponding to the ionic systems, and the corresponding ones for the tethered dimer systems. For this purpose, we have considered the gas and liquid states at coexistence for a temperature $`T0.95T_c`$, with $`T_c`$ the corresponding critical temperature. The ion-ion radial distribution functions of the ionic systems, for $`T^{}=0.0425`$ and the coexisting densities $`\rho ^{}=0.0076(15)`$ and $`\rho ^{}=0.162(3)`$ are plotted in Figs. 2 and 3, respectively. In both cases, the unlike-ion radial distribution function becomes very large close to contact, indicating the association in bound pairs of unlike ions. Moreover, the like-ion radial distribution functions show maxima around $`r^{}=1.5`$ (in the case of $`g_{}`$, this peak coincides with the contact value), and there is a secondary maximum in $`g_+`$ for $`r^{}2.5`$. These observations allow us to conclude that there are high correlation between pairs of associated unlike ions. We recall that the range of densities in which the gas-liquid coexistence occurs prevents packing effects, so the structure is completely given by the Coulombic interactions.

The ion-ion radial distribution functions for the tethered dimer fluid are plotted on the Fig. 4 for the gas branch ($`T^{}=0.0405`$, $`\rho ^{}=0.0119(15)`$) and Fig. 5 for the liquid branch ($`T^{}=0.0405`$, $`\rho ^{}=0.206(2)`$). The comparison between the ionic fluid and tethered dimer fluid microscopic structures confirms the qualitative similarity between both systems. A further test of this similarity is found in the comparison of the neutral cluster populations in the gas branch. We use Gillan’s definition of a cluster Gillan . Two ions $`i`$ and $`j`$ are directly bound when the distance between them is less than $`R_{ij}^c`$. This condition defines mathematically an equivalence relationship, and the equivalences classes in which the ions group are the clusters. As the interaction between like ions is repulsive, the dependence of the cluster definition on $`R_{++}^c`$ or $`R_{}^c`$ should be very weak (if they are taken smaller than the mean distance between two like ions). On the other hand, the cluster definition is not very sensitive to the value of $`R_+^c`$ if that value lies between the first minimum of $`g_+`$ and the mean distance between aggregates. In this work we have used $`R_{++}^c=R_{}^c=R_+^c=1.5\sigma _\pm `$. As expected, the microscopic structure of the gas phase is dominated by the $`N`$-ion neutral clusters Bresme ; E2PM . Their fractions $`f_N`$ for the ionic and tethered ion systems at the gas phase in coexistence for $`T=0.95T_c`$ are quite similar (see Fig. 6), although the tethered dimer systems have slightly higher fractions of neutral clusters than do ionic systems, specially for large $`N`$, in agreement with the results previously reported E2PM .

Despite the similarities found between the ionic fluid and the tethered dimer fluid, there are some differences that can rationalize the quantitative differences between their phase diagrams. The unlike-ion radial distribution functions of the tethered dimer fluid differ qualitatively from the ionic counterparts close to the contact value, since in the former there is a jump at the maximum allowed separation of two ions of a dimer, while the ionic $`g_+`$ takes smoothly higher values on that range of values of $`r^{}`$. Consequently, the condition on the confinement of the ions that compose a tethered dimer must be relaxed in order to refine the pairing concept in the ionic fluids. In the next Section we will introduce an appropriate framework for such a goal.

## III Chemical picture of ionic fluids

We consider an ionic fluid ($`q_+=q_{}q`$) be contained in a volume $`V`$ and in equilibrium with a reservoir at a temperature $`T`$ and chemical potentials $`\mu _+`$ and $`\mu _{}`$. For simplicity, the neutral grand-canonical ensemble, in which only neutral configurations are allowed, will be considered. This ensemble has been shown to be equivalent to the usual grand-canonical ensemble in the thermodynamic limit Lieb . The neutral grand-canonical partition function can be written as:

$$\mathrm{\Xi }=\underset{N_i=0}{\overset{\mathrm{}}{}}\left(\frac{\text{e}^{\beta \mu }}{\mathrm{\Lambda }_+^3}\right)^{N_i}\left(\frac{\text{e}^{\beta \mu }}{\mathrm{\Lambda }_{}^3}\right)^{N_i}Z_U(N_i,N_i)$$
(5)
where $`N_i`$ is the number of ions of each species, $`\beta =1/(k_BT)`$, $`\mu (\mu _++\mu _{})/2`$, $`\mathrm{\Lambda }_+`$ and $`\mathrm{\Lambda }_{}`$ are the thermal de Broglie wavelengths corresponding to each species, and $`Z_U(N_i,N_i)`$ is the canonical configurational partition function corresponding to a fixed number of ions at the temperature $`T`$ and enclosed in the volumen $`V`$:

$$Z_U(N_i,N_i)=\frac{1}{(N_i!)^2}_{V^{2N_i}}𝑑𝐫_1^+𝑑𝐫_1^{}\mathrm{}𝑑𝐫_{N_i}^+𝑑𝐫_{N_i}^{}\text{e}^{\beta U}$$
(6)
with $`UU^{phys}=U_{hc}+_{i<j}q_iq_j/(Dr_{ij})`$ the total potential energy and $`U_{hc}`$ the hard-core contribution, equal to $`+\mathrm{}`$ if two particles overlap, and $`0`$ otherwise (other tempered potentials can be used, but the main results of this Section will remain unchanged). If the fluid particles are strongly associated into bound $`(+,)`$ pairs, as occurs in the low-temperature and low-density region in which the vapor-liquid transition occurs for the ionic fluids, the “physical” representation described above will be inadequate, and it can replaced by a “chemical” picture, in which the fluid is composed by associated pairs and free ions. First, a rule that unequivocally identifies bound pairs for each configuration of the ionic fluid (up to a set of null measure in the configurational space) is needed. Once such a rule is defined, we can write the canonical configurational partition function as:

$$Z_U(N_i,N_i)=\underset{N_p=0}{\overset{N_i}{}}Z_U^{}(N_f,N_f,N_p)$$
(7)
where $`N_p`$ is the number of bound pairs, $`N_f=N_iN_p`$ is the number of free ions of each species, and $`Z_U^{}(N_f,N_f,N_p)`$ is the canonical configurational partition function corresponding to a system of $`N_p`$ associated pairs and $`N_f`$ free ions of each species:

$`Z_U^{}(N_f,N_f,N_p)`$ $`=`$ $`{\displaystyle \frac{1}{N_p!}}{\displaystyle \frac{1}{(N_f!)^2}}{\displaystyle _{V^{2(N_p+N_f)}}}d(1)\mathrm{}d(N_p)`$ (8)
$`\times `$ $`d𝐫_1^+d𝐫_1^{}\mathrm{}d𝐫_{N_f}^+d𝐫_{N_f}^{}\text{e}^{\beta U^{}}`$
where $`(i)\{𝐫^+,𝐫^{}\}_i`$ are the positions of the ions that compose the associated pair $`i`$, and $`𝐫_i^{+()}`$ corresponds to the coordinates of a $`+()`$ free ion $`i`$, respectively. The potential energy $`U^{}U^{chem}`$ *does not* coincide, in general, with the physical potential energy $`U^{phys}`$, since different configurations of the “chemical” mixture can be compatible with a given ionic configuration.

Substituting equation (7) into the equation (5) and rearranging the resulting expression, the neutral grand-canonical partition function can be written as:

$`\mathrm{\Xi }`$ $`=`$ $`{\displaystyle \underset{N_f=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{N_p=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\text{e}^{\beta \mu }}{\mathrm{\Lambda }_+^3}}\right)^{N_f}\left({\displaystyle \frac{\text{e}^{\beta \mu }}{\mathrm{\Lambda }_{}^3}}\right)^{N_f}\left({\displaystyle \frac{\text{e}^{\beta \mu _p}}{\mathrm{\Lambda }_p^6}}\right)^{N_p}`$ (9)
$`\times `$ $`Z_U^{}(N_f,N_f,N_p)`$
with $`\mu _p2\mu `$ and $`\mathrm{\Lambda }_p=\sqrt{\mathrm{\Lambda }_+\mathrm{\Lambda }_{}}`$. These last definitions correspond to the chemical equilibrium conditions between the free ions and the associated pairs Chandler ; Levin . We remark that this derivation is *independent* of the criterion chosen to define the pairs. However, the choice must be such that matches the microscopic structure of the fluid. As we have seen in the previous Section, the vapor structure is dominated by neutral aggregates, mostly dimers and tetramers. This fact is consistent with the Bjerrum ideas of pairing Bjerrum . Consequently, a distance-based criterion seems to be the most appropriate: two unlike ions that are closer than $`R_c`$ ($`R_c`$ being a suitable cutoff distance) are considered as an associated pair. It is easy to see that such a criterion *does not* define uniquely the associated pairs when the population of tetramers and higher order neutral clusters is not negligible, as it can be seen from Fig. 7, since there are different arrangements of associated pairs that are compatible with the same ionic configuration. Hence, a more systematic criterion is needed, in order to be able to define associated pairs for *almost every* ionic configuration, while keeping the intuitive definition of an associated pair as encompassing unlike ions that are at the closest distances. We use a suitable modification of the Stillinger-Lovett pair definition on a given ionic configuration Stillinger . In this prescription, all the distances between two unlike ions are computed, and then the first pair is *defined* as the two unlike ions at closest distance. This step is repeated, taking into account only ions that remain unpaired from previous steps for the evaluation of $`(+,)`$ distances. We stop when the minimum distance between two unlike ions is greater than $`R_c`$, considering the remaining ions as free ions (no free ions were considered in Stillinger , and consequently $`R_c\mathrm{}`$ in that case). This protocol provides an unique configuration of associated pairs and free ions for almost each ionic configuration, since the method is ambiguous only in a subset of the ionic configurations in which at least two $`(+,)`$ distances are equal, and the measure of such a set is null in the configurational space.

Obviously, different pairing prescriptions can be used. However, this protocol provides an explicit expression for the “chemical” potential energy $`U^{chem}`$:

$`U^{chem}`$ $`=`$ $`U^{phys}+{\displaystyle \underset{i=1}{\overset{N_p}{}}}v^p(\{𝐫^+,𝐫^{}\}_i)+{\displaystyle \underset{1i<jN_p}{}}v^{pp}(\{𝐫^+,𝐫^{}\}_i,\{𝐫^+,𝐫^{}\}_j)`$ (10)
$`+`$ $`{\displaystyle \underset{i=1}{\overset{N_p}{}}}{\displaystyle \underset{j=1}{\overset{N_f}{}}}v^{p+}(\{𝐫^+,𝐫^{}\}_i,𝐫_j^+)+{\displaystyle \underset{i=1}{\overset{N_p}{}}}{\displaystyle \underset{j=1}{\overset{N_f}{}}}v^p(\{𝐫^+,𝐫^{}\}_i,𝐫_j^{})+{\displaystyle \underset{i=1}{\overset{N_f}{}}}{\displaystyle \underset{j=1}{\overset{N_f}{}}}v^+(𝐫_i^+,𝐫_j^{})`$
where $`U^{phys}`$ is the physical potential energy due to the hard-core and eletrostatic interactions, and the other terms correspond to effective interactions of entropic origin needed to reduce the configurational space of the “chemical” system. The term $`v^p(\{𝐫^+,𝐫^{}\}_i)`$ depends only on the positions of the ions that compose an associated pair and confines them to be at a distance shorter than $`R_c`$:

$$v^p(\{𝐫^+,𝐫^{}\}_i)=\{\begin{array}{cc}0\hfill & |𝐫_i^+𝐫_i^{}|R_c\hfill \\ & \\ +\mathrm{}\hfill & |𝐫_i^+𝐫_i^{}|>R_c\hfill \end{array}$$
(11)
The pairwise potential energy $`v^{pp}(\{𝐫^+,𝐫^{}\}_i,\{𝐫^+,𝐫^{}\}_j)`$ corresponds to an steric hindrance condition between two associated pairs Stillinger , since the distances between two unlike ions corresponding to different pairs cannot be shorter than the minimum distance between the ions that compose each pair:

$$v^{pp}(\{𝐫^+,𝐫^{}\}_i,\{𝐫^+,𝐫^{}\}_j)=\{\begin{array}{cc}+\mathrm{}\hfill & |𝐫_i^{}𝐫_j^\pm |d_{ij}\hfill \\ & \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(12)
with $`d_{ij}min(|𝐫_i^+𝐫_i^{}|,|𝐫_j^+𝐫_j^{}|)`$.

The interaction between the associated pairs and the free ions is modified by the term $`v^{p\pm }(\{𝐫^+,𝐫^{}\}_i,𝐫_j^\pm )`$, that prevents the free ion to be closer to the unlike ion of the associated pair than its partner:

$$v^{p\pm }(\{𝐫^+,𝐫^{}\}_i,𝐫_j^\pm )=\{\begin{array}{cc}+\mathrm{}\hfill & |𝐫_i^{}𝐫_j^\pm ||𝐫_i^+𝐫_i^{}|\hfill \\ & \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(13)
Finally, two unlike free ions cannot be at a shorter distance than $`R_c`$ due to the term $`v^+(𝐫_i^+,𝐫_j^{})`$:

$$v^+(𝐫_i^+,𝐫_j^{})=\{\begin{array}{cc}+\mathrm{}\hfill & |𝐫_i^+𝐫_j^{}|R_c\hfill \\ & \\ 0\hfill & |𝐫_i^+𝐫_j^{}|>R_c\hfill \end{array}$$
(14)
It is not hard to see that a mixture of associated pairs and free ions with a potential energy given by the Eq. (10) is completely equivalent to the original ionic system. It is interesting to note that, except the term on $`v_p`$ that ties the ions that compose an associated pair, the other terms are pairwise, short-ranged modifications of the hard-core conditions, so the potential energy of an allowed mixture configuration is the same as in the corresponding ionic configuration. Moreover, every configuration of associated pairs and free ions obtained from an ionic configuration by the pairing procedure described above fulfil the constrains induced by the added potentials. Conversely, the mixture configuration is the same as the one obtained from the corresponding ionic configuration by the Stillinger-Lovett protocol.

This analysis provides an exact chemical representation of the ionic system as a mixture of bound pairs and free ions in chemical equilibrium. The equivalence between representations of the system is not only at the level of thermodynamic properties, but also in the microscopic structure, unlike previous studies based on the matching of the “chemical” and the “physical” free energies, e.g. via the virial coefficients Hill ; Zuckerman . In those studies the effective potentials given by Eqs. (11) and (14) could be guessed, but the other terms (that arise from matching third and higher order virial coefficients) were not clearly identified. Furthermore, they have not being used at all in pairing theories of electrolytes. We must stress the importance of the effective new term given by Eq. (12) in the potential energy, that reduces the configuration space available to the associated pairs. This effect is specially important for high values of $`R_c`$. An extreme case case that illustrates the effect of missing the steric hindrance term is the loosely tethered dimer fluid, in which the ions that compose the pair can be at any relative distance greater than $`\sigma _\pm `$. In this case, the canonical partition function of the dimer fluid is $`N_p!\times Z_U(N_p,N_p)`$, where $`N_p`$ is the number of dimers, $`Z_U(N_p,N_p)`$ is given by Eq. (6) and the $`N_p!`$ factor is the number of different ways of pairing the unlike ions in each ionic configuration and lead to different dimer configurations. As a consequence, in the thermodynamic limit the Helmholtz free energy per particle $`f_d`$ of the loosely tethered dimer fluid diverges as:

$$\frac{f_d}{k_BT}\mathrm{ln}N_p+\left(\frac{2f_i}{k_BT}+1\right)$$
(15)
where $`f_i`$ is the ionic Helmholtz free energy per ion, that it is well-behaved in the case of neutral systems Lieb .

The association degree of the ionic fluid in the coexistence region has been studied by computer simulation in the framework of the previous analysis. In the usual grand-canonical simulations, we have identified the associated pairs by the Stillinger-Lovett rule with $`R_c=\sqrt{3}L/2`$, i.e. all the ions are associated, and the minimum image convention is used in order to calculate the distances between unlike ions. The probability distributions $`p(r)`$ of having a pair an internal separation distance $`r`$ between ions is plotted in the Fig. 8 for the gas and liquid phases at $`T^{}=0.0443`$. Surprisingly, both distributions are practically identical, despite the fact that the corresponding densities are quite different. The distributions show a very pronounced peak for $`r=\sigma _\pm `$, and a local maximum around $`r=2.5\sigma _\pm `$, approximately where $`g_+`$ presents the second local maximum. For large values of $`r`$, the distributions $`p(r)`$ decay to zero. These results confirm the strong association in the ionic fluid at low temperatures, and that the structure is weakly affected by variations in density, at least in the range in which the gas-liquid coexistence occurs. For values of $`r`$ close to the contact value, the distribution is well described by the non-interacting pair fluid probability distribution, given by the following expression:

$`p^{ideal}(r)`$ $``$ $`{\displaystyle \frac{r^2\mathrm{exp}\left(\frac{q^2}{Dk_BTr}\right)}{_{\sigma _\pm }^{R_c}y^2\mathrm{exp}\left(\frac{q^2}{Dk_BTy}\right)𝑑y}}`$
$`{\displaystyle \frac{q^2r^2}{Dk_BT\sigma _\pm ^4}}`$ $`\times `$ $`\mathrm{exp}\left({\displaystyle \frac{q^2}{Dk_BT\sigma _\pm }}\left[{\displaystyle \frac{\sigma _\pm }{r}}1\right]\right)`$ (16)
valid for $`T^{}1`$ and an internal cutoff $`R_c\sigma _\pm `$. However, for $`r1.5\sigma _\pm `$, $`p(r)`$ deviates from the ideal expression as a consequence of the interaction between associated pairs. The local maximum showed by the distribution function indicates that the most bound pairs are solvated by less bound pairs, to form stable neutral tetramers (and higher order clusters, but their inclusion does not seem to affect the conclusions of this Section). This fact is in agreement with the features observed in the pair correlation functions in the previous Section. We must stress that the structure that $`p(r)`$ presents is only due to the Coulombic interactions.

It is instructive to compare our results in the coexistence region to the higher temperature ones. We have computed $`p(r)`$ at $`T^{}=1`$ and the same range of densities (see Fig. 9). The qualitative behavior of $`p(r)`$ at high temperatures is similar to that predicted by Stillinger and Lovett Stillinger2 . First, it shows a maximum localized in $`r\rho ^{\frac{1}{3}}`$. On the other hand, $`p(r)`$ takes significant values in a wider range than the corresponding functions at lower temperatures, decaying for large values of $`r`$ as $`r^n`$, with $`n3`$. The latter prediction differs from the value $`n=4`$ given in Ref. Stillinger2  and requires further study to elucidate it.

The probability distribution $`p^{R_c}(r)`$ of having a pair an internal separation distance between ions when the cutoff distance to define a pair is equal to $`R_c`$ can be obtained from $`p(r)`$ as:

$$p^{R_c}(r)=\frac{p(r)}{\mathrm{}(R_c)}=\frac{p(r)}{_{\sigma _\pm }^{R_c}p(t)𝑑t}$$
(17)
where $`\mathrm{}(r)_{\sigma _\pm }^rp(t)𝑑t`$ is the cumulant distribution corresponding to $`p(r)`$. The density of associated pairs is equal to $`\rho \mathrm{}(R_c)/2`$, and the total density of free ions is $`\rho [1\mathrm{}(R_c)]`$, where $`\rho `$ is the total density of ions in the physical picture. It can be seen from Fig. 8 that for $`R_c3\sigma _\pm `$, more than 95% of ions are associated into pairs in both vapor and liquid phases. So, if the associated pair fluid phase separate in the same $`(T,\rho )`$ region as the ionic fluid, it is expected than the free ions do play a mere perturbative role in the phase coexistence. Actually, the free ionic subsystem is a non-additive binary charged hard-sphere mixture, where $`R_c`$ plays the role of the unlike-ion collision diameter. However, the interactions between like ions are purely repulsive, and it is not expected differences between the behavior of this mixture and the RPM as soon as $`\sigma _+`$ and $`\sigma _{}`$ are less than $`R_c`$. Furthermore, this system is in a polar environment given by the associated pair subsystem, and consequently the effective dielectric constant $`D_{eff}`$ of the background will be increased. Then the free ion subsystem is expected to have the same behavior as the RPM at the effective reduced temperature $`\stackrel{~}{T}`$:

$$\stackrel{~}{T}=\frac{D_{eff}k_BTR_c}{q^2}=\frac{D_{eff}}{D}\times \frac{R_c}{\sigma _\pm }\times T^{}$$
(18)
where $`T^{}`$ is the reduced temperature. For $`R_c3\sigma _\pm `$, the critical temperature reduces at least to one third of the RPM reduced critical temperature. Consequently, any possible free ion-driven vapor-liquid transition should happen at much lower temperatures, and thus the free ion subsystem should play no role in the vapor-liquid transition. In order to check this hypothesis, the phase diagram of the associated pair fluid has to be obtained. This issue will be addressed in the next Section.

## IV Phase Behavior of the Associated Pairs Fluid

The results obtained in the previous Section suggest that the “chemical” picture introduced above is a very convenient description of the ionic fluid structure. Furthermore, we can analyse the role played by the associated pairs in the phase equilibrium by eliminating the free ions. We have performed grand-canonical Monte Carlo simulations of the associated dimer system for $`\lambda =3`$ and different values of $`R_c`$. No free ions are allowed ($`N_f0`$) and the potential energy of one configuration is given by the Eq. (10). The grand-canonical free energies corresponding to these systems will consequently provide an upper bound to the ionic fluid grand-canonical free energy at the same temperature and pair chemical potential. As for the ionic system, a fine-discretization approach with a refinement parameter $`\zeta =10`$, and Ewald summation technique with conducting boundary conditions is used to take into account the long range character of the Coulombic interations. The basic steps are either insertion or deletion of associated pairs (chosen randomly with the same probability), biased with a Boltzmann distribution that depends on the separation between the ions that compose the dimer. An associated pair is inserted in the following way: the negative ion is placed at a random place in the box, and its counterion is placed at a relative position $`𝐫_\pm `$ ($`\sigma _\pm |𝐫_\pm |<R_c`$) following a probability distribution proportional to $`\mathrm{exp}\left(\beta _0q^2\varphi _0(𝐫_\pm )\right)`$, where $`\varphi _0(𝐫)`$ is the Ewald potential energy between two unlike ions of unit charge. On the other hand, a pair is deleted with a probability proportional to $`\mathrm{exp}\left(\beta _0q^2\varphi _0(𝐫_\pm )\right)`$, being $`𝐫_\pm `$ the relative position of the positive ion of the pair with respect to the negative one. The value of $`\beta _0`$ has been adjusted to improve the sampling during the simulation. In this work we have set $`\beta _0=8`$.

The acceptance probabilities of a pair insertion $`W_{ij}^i`$ and a pair deletion, $`W_{ij}^d`$ are the following:

$`W_{ij}^i`$ $`=`$ $`\text{min}(1,(L^{})^3\mathrm{exp}(\beta \mu \beta \mathrm{\Delta }U)\left[{\displaystyle \underset{\sigma _\pm |𝐫_\pm |<R_c}{}}{\displaystyle \frac{\mathrm{exp}(\beta _0q^2\varphi _0(𝐫_\pm ))}{\zeta ^3}}\right]\left[{\displaystyle \underset{k=1}{\overset{N_p+1}{}}}\mathrm{exp}(\beta _0q^2\varphi _0(𝐫_\pm ^k))\right]^1)`$ (19)
$`W_{ij}^d`$ $`=`$ $`\text{min}(1,{\displaystyle \frac{\mathrm{exp}(\beta \mu \beta \mathrm{\Delta }U)}{(L^{})^3}}\left[{\displaystyle \underset{\sigma _\pm |𝐫_\pm |<R_c}{}}{\displaystyle \frac{\mathrm{exp}(\beta _0q^2\varphi _0(𝐫_\pm ))}{\zeta ^3}}\right]^1\left[{\displaystyle \underset{k=1}{\overset{N_p}{}}}\mathrm{exp}(\beta _0q^2\varphi _0(𝐫_\pm ^k))\right])`$ (20)
where $`i`$ and $`j`$ are the initial and final configurations, respectively, $`N_p`$ is the number of pairs at the configuration $`i`$, $`\beta =1/(k_BT)`$, $`\mu `$ is the configurational chemical potential, and $`\mathrm{\Delta }U`$ is the “chemical” potential energy variation in the movement, including the steric hindrance terms given by Eq. (12). As it can be seen, the acceptance probabilities reduce to the usual ones as $`\beta _00`$.

Effective critical points for different values of $`L^{}`$ were estimated by using mixed-field finite-size scaling methods Bruce  and assuming Ising-like criticality. Although recent results Fisher2  indicate that the pressure should also enter the field mixing, our approach should be satisfactory to discern the dependence of the critical parameters on $`R_c`$. Moreover, a systematic study for the RPM case has shown that the assumed universality class is very likely to be the right one Luijten . We use histogram reweighting techniques Ferrenberg  to combine the histograms from different runs (typically three) and estimate the critical parameters and their standard errors. Very long runs are needed in order to overcome the critical slowing down and the low acceptance ratios due to the low temperatures involved. Defining a step as a try of a pair insertion/deletion, we have performed simulations of about $`2\times 10^74\times 10^7`$ equilibration steps and $`2\times 10^84\times 10^8`$ sampling steps for $`L^{}=12`$; $`5\times 10^72\times 10^8`$ equilibration steps ans $`6\times 10^81.2\times 10^9`$ sampling steps for $`L^{}=15`$; and $`2\times 10^8`$ equilibration steps and $`1.2\times 10^91.8\times 10^9`$ sampling steps for $`L^{}=18`$, in order to get good statistics for the histograms. The (effective) critical parameters are obtained by minimizing the deviation of the appropriately scaled mixed-field $`\rho su`$ marginal probability distribution ($`u`$ is the potential energy density and $`s`$ is the mixing parameter) with respect to the corresponding critical 3-dimensional Ising universal function Bruce ; Thanos3 . As also occurs in the ionic and tethered dimer systems, the matching improves as $`L^{}`$ increases (see Fig. 10). Our estimations of the critical parameters for different values of $`L^{}`$ and $`R_c^{}R_c/\sigma _\pm `$ are listed in Table 1 and plotted in Figs. 11 and 12. The smallest value of $`R_c^{}`$ is the same as the maximum allowed separation in the tethered dimer system, and as expected the results obtained for the associated pairs fluid are practically indistinguishable from the tethered dimer system. As $`R_c^{}`$ increases, it is observed a sharp increase in the critical temperature but the critical density remains almost unaffected. However, for $`R_c^{}34`$ the critical temperature reaches a maximum and then decreases, converging towards the free ion critical temperature. On the other hand, the critical density also decreases towards the free ion critical density, although the statistical uncertainties are bigger for this parameter and the exact path of convergence is less defined. We must also point out that the statistical uncertainties are big enough not to make possible an estimation of correction-to-scaling effects.

For an heuristic explanation for such a behavior, we have to consider the energetically favoured configurations. It is important to consider the interaction between two different associated pairs, that has its minimum energy configuration in a square ring conformation for low values of $`\lambda `$, with a close-energy secondary minimum conformation for the linear $`(+,)`$ chain E2PM . If $`R_c`$ is too small, by deforming the square it is possible to find energetically relevant configurations (i.e. the potential energy difference per ion of such configurations with respect to the minimum one is of the same order as the thermal energy) in which one of the associated pairs has an internal ionic separation less than $`R_c`$, but the other one does not fulfil such a condition. Consequently, a greater value of $`R_c`$ also increases the number of energetically relevant configurations between two associated pairs. Then the effective interaction between associated pairs is enhanced, and thus an increase of the critical temperature should be expected until all the relevant configurations are allowed (that occurs for $`R_c3\sigma _\pm `$). On the other hand, the average size of the associated pairs and their aggregates is increased, leading to a decrease on the critical density (in order to keep the reduced density in terms of the effective particle size constant). These arguments are not longer valid for $`R_c^{}3`$. Our results show the inclusion of associated pairs with larger internal ionic separations are not crucial to the gas-liquid phase transition but also that they bar the phase transition, as it can inferred from the critical parameters decrease.

As for the ionic and tethered dimer models, histogram reweighting techniques Ferrenberg ; Thanos3  allow us to obtain the coexistence curve up to $`T0.98T_c`$. As the finite-size effects are not important far from the critical point, we have considered $`L^{}=15`$. Taking advantage of the wide range of densities covered by near-critical histograms, we combine them with liquid subcritical state histograms in order to extend the density range. The extra simulations involve shorter runs (typically $`2\times 10^8`$ steps after equilibration). The gas-liquid coexistence curves for different values of $`R_c^{}`$ are plotted in Fig. 13, showing the same tendency as the observed one in the critical parameters. It is interesting to note that the vapor branches coincide (except close to the critical point) for $`R_c>3\sigma _\pm `$. We conclude from this observation that is in the liquid branch where the ionic fluid differs mostly from the associated pair system.

## V Discussion and Conclusions

In summary, an exact “chemical” representation of the ionic system as a mixture of $`(+,)`$ associated ion pairs and free ions has been introduced. This representation, closely related to the Stillinger-Lovett pairing procedure for the RPM, has the advantage that not only provides an exact matching of between the physical and chemical representation thermodynamics, but also from a microscopic point of view. It also avoids an entropy catastrophe that occurs for the tethered dimer model studied at large values of the tether length. The addition to the physical hamiltonian of some new pairwise hard-core interactions between the “chemical” components provides the suitable hamiltonian for the “chemical” representation. In the low temperature and low density regime in which the ionic vapor-liquid transition occurs, such a representation provides a faithful characterization of the microscopic structure of the ionic fluid, and it can be the basis of new theoretical approaches.

The analysis of the phase behavior of the system only composed by associated pairs indicates strongly that the ionic fluid vapor-liquid transition is driven by them. For $`R_c3\sigma _\pm `$, the critical temperature increases with $`R_c`$, as more energetically favoured configurations are allowed. For $`R_c3\sigma _\pm `$, the critical parameters converge smoothly *from above* towards the ionic fluid critical parameters as $`R_c`$ increases. This fact indicates not only that the free ions do not drive the phase transition, but also have the opposite effect in the transition. The value of $`R_c3\sigma _\pm `$ corresponding to the maximum on the critical temperature, can be regarded as the optimal size of the associated pair.

Some remarks on the limitations of the present work are appropriate at this point. We have focussed only on the effect of the associated pairs in the thermodynamical properties. However, the conducting character of the ionic fluid is completely driven by the free ions. On the other hand, the remarkable success of pairing theories (even when the “chemical” representation they implicitly use is not completely correct) remains unexplained. It is possible that the associated pair solvation, in addition to the ion-ion correlations, can mimic the pair-pair interactions. Further studies are needed to completely solve the origin of the ionic fluid vapor-liquid phase transition.

###### Acknowledgements.

J.M.R.-E. and L.F.R. gratefully acknowledge financial support for this research by Grant No. PB97-0712 from DGICyT (Spain) and No. FQM-205 from PAI (Junta de Andalucía). J.M.R.-E. also wishes to thank an FPI scholarship from Ministerio de Educación y Cultura (Spain). AZP acknowledges financial support from the Department of Energy (grant DE-FG02-01ER1512) and ACS-PRF (grant 38165-AC9).