text:
# Liquid-liquid equilibrium for monodisperse spherical particles

## I Introduction

Traditionally, liquid-liquid separation for single component systems was not considered as a real possibility, because of the rapidly changing liquid structure, which smoothly varies on temperature and pressure. This view is changing. Experimentally, a whole family of network forming fluids, including some tetrahedrally coordinated materials, is suspected to have transitions between liquid states. Usually this transition occurs well inside the supercooled region, where the liquid is already in a glassy state, and then it is related to polyamorphism. The numerical modeling of these materials shows that interaction potentials can be constructed for which there are transitions between different liquid states. For instance, the potentials used to mimic the properties of water generate two liquid phases in the supercooled region. The main characteristic of these materials (and of the interaction potentials used to model them) which is responsible for their anomalous properties is the possibility of generating open, low coordinated network structures at low pressure, while collapsing to more compact structures at higher pressures.

This key ingredient is captured by a class of model spherical potentials, the so-called soft core Stell-Hemmer potentials. They include a strict hard core at some distance $`r_0`$ and a soft core at a larger distance $`r_1`$. In this way, two typical distances between particles occur, and there is a collapse from the largest to the smallest distance on applying pressure. The properties of these systems give a remarkably consistent description of most of the anomalous behavior of tetrahedrally coordinated materials.

The possibility of liquid-liquid separation has been extended in the last time by the finding that it can occur even at thermodynamic equilibrium. Katayama et al. report that phosphorus (which is a network forming fluid) shows liquid-liquid coexistence at equilibrium, providing the first experimental evidence of this phenomenon. In addition, Glosli and Ree find numerical evidence of liquid-liquid coexistence at equilibrium in carbon. This was the first direct numerical evidence of the phenomenon for a single component system. The interparticle potential used to model carbon-carbon interaction is non-isotropic, and incorporates information of the internal (electronic) degrees of freedom of the atoms. It is shown that one of the two coexisting liquids is mainly $`sp`$ coordinated, whereas the other shows a dominant $`sp^3`$ character.

On the basis of these results, it is natural to ask whether some kind of soft-core potentials (with no internal degrees of freedom for the particles) can produce an equilibrium liquid-liquid transition. We will show, by providing concrete examples, that soft core Stell-Hemmer potentials are able to do that, then showing that the phenomenon may occur also for systems that, in addition to being monodisperse, have particles which are spherical and ‘rigid’, i.e., in cases where internal degrees of freedom play no role. This is of basic interest since this model is much simpler, and then more transparent to analyze than real examples of material having liquid-liquid coexistence at equilibrium. It also provides a benchmark for the study of tetrahedrally coordinated materials with liquid-liquid critical points in the supercooled region. The relation between the existence of the liquid-liquid critical point and the anomalous thermodynamic properties of these materials can be studied here in a context of thermodynamic equilibrium. On a different context, our results make plausible the existence of liquid-liquid separation at equilibrium in particular cases of colloidal dispersions, which is of potential technological interest.

## II The model

We study a model of particles interacting through a potential which is the sum of attractive and repulsive parts $`V(r)=V_A(r)+V_R(r)`$. The repulsive part is given as

$`V_R(r)=\mathrm{}\mathrm{for}r<r_0`$ (1)
$`V_R(r)=\epsilon _0\left(r_1r\right)/\left(r_1r_0\right)\mathrm{for}r_0<r<r_1`$ (2)
$`V_R(r)=0\mathrm{for}r>r_1,`$ (3)
namely, a hard core at distance $`r_0`$, and a linear (repulsive) ramp between $`r_0`$ and $`r_1`$. The attractive part is taken also to be linear in $`r`$, of the form

$`V_A(r)=\gamma \left(r_2r\right)/\left(r_2r_0\right)\mathrm{for}r<r_2`$ (4)
$`V_A(r)=0\mathrm{for}r>r_2`$ (5)
($`r_2>r_1`$, $`\gamma >0`$). The repulsive part of the potential favors one of the two distances $`r_0`$ or $`r_1`$ between neighbor particles depending on the value of external pressure and the attraction intensity $`\gamma `$. Interparticle distance $`r_1`$($`r_0`$) is favored for pressures lower (higher) than some crossover value $`P_0.`$ The crossover pressure (for $`\gamma =0`$) is given roughly as $`P_0=\epsilon _0/(v_1v_0)`$, where $`v_0`$ and $`v_1`$ are the specific volumes of the states with nearest particles at distances $`r_0`$ or $`r_1`$. For $`P`$ near $`P_0`$ the system has an anomalously large compressibility. The attractive part of the potential may turn this anomaly into a first order transition through a van der Waals mechanism.

We will show results for two-dimensional (2D) and three-dimensional (3D) systems. The parameters we use are: $`r_1=1.72r_0`$, $`r_2=4.8r_0`$, and $`\gamma =0.27\epsilon _0`$ for 2D and $`r_1=1.72r_0`$, $`r_2=3.0r_0`$, and $`\gamma =0.31\epsilon _0`$ for 3D. Temperature will be measured in units of $`\epsilon _0/k_B`$, and pressure in units of $`\epsilon _0/r_0^2`$ in 2D, and $`\epsilon _0/r_0^3`$ in 3D. Calculations were done by standard Monte Carlo techniques, in systems with periodic boundary contitions. Some of the results to be presented where done in the NVT ensemble, where the volume of the system is kept fixed at each simulated value, and pressure is calculated through standard formulae. Other results correspond to the NPT ensemble, where external pressure is fixed, and the volume of the system is taken as an additional variable during the simulation.

## III Results

### A Two dimensions

Fig. 1 shows the isotherms we get from the simulations for a system of 200 particles in the 2D case, near the position in which we find a liquid-liquid critical point for our parameters, namely $`T_{LL}0.10`$, $`P_{LL}0.07`$, and $`v_{LL}1.85r_0^2`$. We see in fact that when reducing $`T`$ below $`0.1`$, the isotherms get a loop, typical of a first order transition. This transition separates two liquid phases that we call (following the notation in water) low density liquid (LDL) and high density liquid (HDL). We have checked that the position and characteristics of this transition do not change when going to a system of 1000 particles. As an example, we show in Fig. 2 the $`T=0.105`$ and $`T=0.0875`$ isotherms of systems of 200 and 1000 particles. Isotherms are seen to coincide with the 200 particles case, except within the coexistence region. Here the isotherm is more flat in the case of larger system. This is in fact what expected, since in an infinitely large system isotherms are strictly flat in the coexistence region. Fig. 3 shows snapshots of the system at $`T=0.0875`$, within the LDL and HDL phases and in the coexistence region. The HDL (LDL) phase is characterized by a larger amount of particles at distance $`r_0`$ ($`r_1`$) from their neighbors. The snapshot at coexistence \[Fig. 3(b)\] shows indeed that different neighborhoods of high density and low density coexist in the system. The existence of a loop in the isotherms below the critical temperature, instead of a flat region, even when the system clearly separates in different regions corresponding to the two different phases \[Fig. 3(b)\] is due to the non-vanishing contribution to the total energy of the surface energy between the two phases. This contribution only vanishes in infinite systems.

From the isotherms of Fig. 1 we see that the LDL has a density anomaly, since there is a temperature range in which volume increases as temperature is reduced, at constant pressure. This density anomaly (which is reminiscent of the similar phenomenon in water) is known to occur in core-softened models even in the absence of an attractive part of the interparticle potential, namely, when $`U_A=0`$.

We will now present the equilibrium $`P`$-$`T`$ phase diagram of the system, close to the liquid-liquid critical point, to show clearly that the liquid-liquid coexistence we are observing occurs in thermodynamic equilibrium, and it is not preempted by the appearance of any other phase (gaseous or crystalline).

The $`P`$-$`T`$ phase diagram is presented in Fig. 4, and the different borders between phases were obtained as explained in the following. The position of the equilibrium first order line between LDL and HDL phases in the $`P`$-$`T`$ phase diagram was determined from the isotherms of Fig. 1 and from the values of the enthalpy $`he+Pv`$ (not shown) by thermodynamic integration, calculating differences in Gibbs free energy $`G`$ from the relation

$$\frac{G_2}{T_2}\frac{G_1}{T_1}=_1^2\frac{v}{T}𝑑P\frac{h}{T^2}dT$$
(6)
where 1 and 2 stand for two set of values $`P_1`$, $`T_1`$ and $`P_2`$, $`T_2`$ and the integration is through an arbitrary (reversible) path in th $`P`$-$`T`$ plane. In this way we can compare $`G`$ on both sides of the first order line, performing an integration surrounding the critical point, and the first order line can be determined with good precision. This procedure is more reliable than the standard Maxwell construction, although the results we have obtained in this case are almost equivalent.

Concerning the crystalline structures, we know in fact that the present model shows a variety of different crystalline structures, which are stable in different regions of the $`P`$-$`T`$ plane. When testing for the stability of these structures, the fact that some of them take a very long time to appear spontaneously in simulations that start in the fluid phase is usually a problem, and then the equilibrium coexistence line between these phases and the fluid has to be determined in some other way. We then applied the procedure already used in  consisting in the inclusion of an additional external potential, with an spatial periodicity corresponding to the crystalline lattice to be studied. The existence of this term (with a strength given by some parameter $`W`$) allows to construct a path in the $`P`$-$`T`$-$`W`$ space that takes the system smoothly (i.e., without crossing any phase transition) from the fluid to the crystalline phase, and the difference in free energy $`\mathrm{\Delta }G`$ between the two phases at some fixed values $`P=P_0`$ and $`T=T_0`$ can be calculated by a generalization of expression (6), namely

$$\mathrm{\Delta }(G/T)|_{P_0,T_0}=_1^2\frac{v}{T}𝑑P+\frac{e_W}{TW}dW\frac{h+e_W}{T^2}dT,$$
(7)
where $`e_W`$ is the potential energy per particle in the external potential. In order to move the system from the crystalline to the fluid phase, the integration consists typically in the following close path: increasing $`W`$ from zero to some large value; increasing $`T`$ from $`T_0`$ to a large value; decreasing $`W`$ down to 0; decreasing $`T`$ down to $`T_0`$. In this way, the difference in free energy between the fluid phase and the possible crystalline phases was determined at fixed points in the $`P`$-$`T`$ phase diagram. Then $`\mathrm{\Delta }G`$ can be calculated for the whole plane by thermodynamic integration of the pure phases. Using this procedure we have verified that the only crystalline phase that is thermodinamically stable within the $`P`$-$`T`$ sector shown in Fig. 4 is the triangular structure, with a lattice parameter $`r_0`$. Other crystalline structures appear only for lower temperatures of higher pressures.

Note in Fig. 4 that the melting of the crystal is anomalous in the region near liquid-liquid coexistence, i.e, both liquid phases are denser than the solid around this point, which is indicated by the negative slope of the liquid-crystal coexistence line. This is probably not a general rule, but is the case we found in all simulations we performed. In this figure the line of maximum density is also indicated. Note that this line does not end exactly at the critical point, but somewhere on the first order line, within the LDL region.

At very low pressure the condensed (liquid or solid) phases are unstable against the formation of a gaseous phase, as in a standard fluid. The condensed-gas coexistence line can be numerically determined in a reliable way by the standard simulation methods we use only close to the (liquid-gas) critical point, where the metastability range is small. Typical isotherms around this point are shown in Fig. 5, from which the location of the liquid-gas critical point is estimated to be at $`P_{LG}0.02`$, $`T_{LG}0.37`$. Using the position of this critical point, we estimate the liquid-gas coexistence line by fitting to a van der Waals equation, chosen to give the correct values of $`T_{LG}`$ and $`P_{LG}`$. This is the line that is seen in the right part of Fig. 4 (note the change in the horizontal scale in this sector). At temperatures ($`0.1`$) where we observe liquid-liquid coexistence, the transition pressure to the gas phase is of the order of $`10^6`$. Then we are sure that the liquid-liquid transition occurs in a region in which liquid is thermodynamically stable with respect to gas.

The crystal-liquid-liquid triple point we obtained turns out into a gas-liquid-liquid triple point (still at equilibrium) by increasing the intensity of the attraction between particles. In fact, this moves all coexistence lines to lower pressures compared to the condensed-gas line, and at a certain point the crystal-liquid-liquid triple point is forced to be in the metastable region with respect to the gas phase. At this point, the gas-condensed line, together with the liquid-liquid line determine a liquid-liquid-gas triple point. For instance, increasing $`\gamma `$ from 0.27 to 0.285 produces a liquid-liquid-gas triple point located at $`P10^6`$, $`T0.10`$. The liquid-liquid critical point is in this case located at $`P0.03`$, $`T0.12`$. At still higher values of the attraction the whole liquid-liquid line moves into the metastable $`P<0`$ region, and the equilibrium liquid-liquid transition is lost.

### B Three dimensions

For the 3D case, we show in Fig. 6 the isotherms obtained for a system of 300 particles. They again show the existence of the liquid-liquid critical point. The curves are smooth for $`T0.08`$, whereas there is a jump in $`v`$ at the coexistence pressure for $`T`$ lower than this value. Also apparent from these curves is the existence (as in 2D) of density anomalies for the low density liquid.

The structure factor $`S(k)`$ of the 3D liquid is shown for the two coexisting liquids at $`T=0.075`$ in Fig. 7. The two main peaks corresponding to interparticle distance $`r_0`$ and $`r_1`$ are indicated. There is an abrupt change in the relative weight of these two peaks at the liquid-liquid transition, corresponding to an abrupt change in the mean distance between first neighbor particles.

For the simulations shown in Fig. 6 we did not find evidence of the appearance of any crystalline structure, suggesting that the liquid-liquid coexistence occurs at equilibrium. To be sure, we should investigate the stability of the different possible crystalline structures, as we did in the 2D case. For the 3D system however, this is a much harder work, as we actually do not know all the crystalline structures for this case (see ). Then we take another route, seeking for other indicators of stability of the liquid phase, to rule out the possibility that we are observing a liquid-liquid transition in a metastable region. So we will look for the fulfillment of the Hansen-Verlet and Lowen criteria for crystallization. Although they stand only as approximate, and mainly heuristic, they have been verified in a variety of cases. For simple (i.e. Lennard Jones) systems, the Hansen-Verlet criterion states that freezing occurs when the first peak of $`S(k)`$ reaches a (quasi-)universal value of $`2.85`$. A straightforward generalization of this criterion to ‘anomalous’ fluids tells that this relation should be checked for all peaks in $`S(k)`$. The results for $`S(k)`$ of Fig. 7 show that all peaks are lower than this value, suggesting liquid stability. The Lowen criteria for freezing indicates that this occurs when the relation between long and short times single particle diffusion coefficient drops below a value of $`0.1`$. This dynamical criterion is independent of the nature of the crystalline structure the system freezes to. Although it has been proposed and verified in systems with Brownian dynamics, we have verified in simple cases (i.e. hard spheres) that it is also fulfilled with the Monte Carlo dynamics, and so we apply it tentatively here. The relation between long and short time values of the diffusivity $`D/D_0`$ (calculated with the Monte Carlo dynamics) is shown in Fig. 8 along the $`T=0.075`$ and $`T=0.09`$ isotherms. The relation $`D/D_0`$ is in the whole range well above the value of $`0.1`$ expected at freezing. Then all indicators coincide in telling that the liquid-liquid coexistence we are observing occurs in truly thermodynamic equilibrium. From Fig. 8 we also see that the evolution of $`D`$ on density shows a maximum which is well known to occur in tetrahedrally coordinates materials  and it is also typical of core-softened potentials.

## IV Conclusions

In conclusion, we have shown that a system of identical, spherical, and rigid particles is able to display liquid-liquid coexistence at equilibrium (both in two and three dimensions). We showed explicitly this for an interaction potential that favors two different equilibrium distances between particles. The study of the present (or some related) model displaying liquid-liquid coexistence at equilibrium is likely to contribute to the understanding of the same phenomenon in cases where it occurs in the supercooled region, where it is much harder to study due to metastability against crystallization and long relaxation times. Our findings also open the possibility of observing liquid-liquid coexistence in monodisperse systems of rigid and spherical particles, i.e., colloidal systems, where interaction potentials with more than one equilibrium distance are easily obtained.

## V Acknowledgments

This work was financially supported by Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Argentina. Partial support by Fundación Antorchas is also acknowledged.


format: text

text:
# 1 Introduction

## 1 Introduction

The SLAC B Factory asymmetric $`B`$-factories was approved by President Clinton in October 1993. B-factories The inaugural meeting of the detector collaboration was held at SLAC in December 1993, at which time the 9 nation BaBar collaboration began. First funding for the machine construction was released in January 1994, launching the SLAC/LBNL/LLNL PEP-II project collaboration . The detector technical design report was completed in March 1995 .

PEP-II was completed in July 1998 and first collisions were observed on July 22, 1998. Commissioning without BaBar concluded in February 1999, at which time the BaBar detector commissioning ceased and BaBar was moved from the off-beamline to the on-beamline position. Running commenced in early May 1999, and first collisions with BaBar were observed on May 26, 1999.

In this report, the status of the PEP-II machine and the BaBar detector are presented. Performance will be presented as of the Lepton Photon Conference (August 1999) with subsequent performance levels footnoted as appropriate. This report is made on behalf of the SLAC/LBNL/LLNL team that built and commissioned PEP-II, the SLAC Accelerator Department that commissioned and operate the machine, and the 71 institution, 9 nation detector collaboration that built, commissioned, and operate BaBar.

## 2 PEP-II Machine

The PEP-II PEP-II design parameters are shown in Table 1. The machine consists of two storage rings, vertically stacked—the high-energy (9 GeV) ring stores electrons and the low energy (3 GeV) ring stores positrons (see Fig. 2). Collisions occur in a single interaction region (IR) (see Fig. 2), at which point the upper low energy ring (LER) has been brought down into the plane of the high energy ring (HER).

Collisions are head-on; beam separation is done magnetically. Electrons and positrons are supplied by the SLAC linac. Extraction and injection is done at the collision energy, so that no ramping of the stored beam is required. The RF system components are identical for the two rings—1.2 MW CW klystrons and warm copper cavities, which are actively damped. The on-cavity damping is not sufficient to ensure stable beams at high currents. Accordingly, bunch-by-bunch damping is supplied in both the transverse and longitudinal planes. The damping is primarily needed to ameliorate the effects of multi-bunch instabilities. A much more detailed description of PEP-II can be found in .

The machine group took on aggressive goals with respect to early peak and integrated luminosity performance. To achieve these goals, machine commissioning was phased in time. The high energy and low energy injection systems were completed early. The HER was completed in July 1997, allowing commissioning of the HER in 1997 and also testing of many systems that were common to the LER. The full machine was completed in July 1998 with the LER and IR being the last systems to finish. However, injection into a partial LER was achieved in early 1998. The phased approach to commissioning was an effective mechanism for early success with the full machine. Two weeks after completion of the full machine, the first collisions were observed.

## 3 BaBar Detector

The BaBar BaBar detector design performance specifications are shown in Table 2. An isometric view of the detector is shown in Fig. 4, indicating the silicon vertex detector silicon vertex detector, drift chamber, Cherenkov radiation particle identification system, cesium iodide calorimeter, superconducting 1.5 Tesla magnet, and instrumented flux return . The accelerator magnets in the IR come within 20 cm of the collision point. In fact, the vertex detector is actually mounted on the final accelerator magnets.

The vertex detector, one-half of which is shown in Fig. 4, is comprised of 5 layers of silicon strip detectors read out with radiation-hard electronics. The drift chamber is a conventional, small cell design with 8 layers of axial wires interleaved with 4 layers of $`\pm `$ 30 stereo wires. A low mass helium-based gas reduces multiple scattering. Surrounding the drift chamber is the DIRC Cherenkov detector system, 48 quartz bars running azimuthally. Cherenkov light generated in the quartz bars is propagated along the bars, exiting at the rear end of the detector where the exit angle is mapped by an array of 11,000 phototubes (see Fig. 6). The pattern of phototube hits can be used to measure the Cherenkov angle and thereby identify the radiating particle. Photons and electrons are identified in the cesium diodide calorimeter cesium iodide calorimeter. Because of the asymmetry of the machine, there is no need for electromagnetic calorimeter coverage in the backward region. The steel superstructure used to return the magnetic flux is made in layers interspersed with resistive plate chamers resistive plate chambers. This instrumented flux return is used for muon and $`K_L^0`$ identification. BaBar has chosen a pioneering approach to data acquisition and data analysis software. All software (both the coding and the databases) uses the object oriented approach The accelerator uses photons from the process $`e^+e^{}e^+e^{}\gamma `$ to derive a fast luminosity signal. However, the Bhabha rate is high enough at moderate luminosities to provide an online luminosity measurement from the detector. Luminosities quoted here come from the Bhabhas in BaBar.

## 4 Machine and detector performance

The month of June 1999 was used to integrate the machine and detector programs into a system which could run steadily for physics. The luminosity was improved and machine related backgrounds were understood. The Upsilon (4s) resonance $`\mathrm{{\rm Y}}(4s)`$ resonance was mapped out in a 3-day run as shown in Fig 6. Using about 1/3 of the data, resonance parameters of $`M(\mathrm{{\rm Y}}_{4s})=10.5841\pm 0.0007`$ GeV and $`\mathrm{\Gamma }(\mathrm{{\rm Y}}_{4s})=11.1\pm 3.4`$ MeV were obtained. In July 1999, persistent physics running began, interspersed with prolonged periods for machine development. Performance results for the machine as of the time of LP99 are given in Table 3. In the table, typical running parameters with BaBar are contrasted with the best that was achieved during commissioning without BaBar. As of the end of July 1999, BaBar had logged about 200 pb-1, as can be seen from Fig. 8. As of the end of November 1999, BaBar had logged 2 fb-1, and the record peak luminosity was $`1.43\times 10^{+33}`$ cm-2 sec-1.

Machine related backgrounds were reasonable; at luminosities above $`5\times 10^{32}`$ cm-2 sec-1, all device occupancies were less than 20% of conservative allowances. Radiation dose levels were low; in the case of the highest radiation region of the Silicon Vertex Tracker (SVT), dose rates were 20% of the pre-ordained budget.

As can be seen from Fig. 8, the detector data logging efficiency is high, particularly at this early phase of operation. Typically, BaBar logs $``$ 90% of the delivered luminosity. It is much too early to expect physics results. However, some plots are shown which are representative of the process of tuning up and monitoring the performance of BaBar. This is an evolving process and much work remains to be done.

Figure 9 shows early performance measures for the SVT. If one excludes 8 non-working modules (out of 208), detector efficiencies are seen to be excellent. The hit resolution in $`\varphi `$ and $`z`$ are not far from that expected.

Figure 10 shows the drift chamber single-hit resolution as a function of distance from the wire. This resolution is very close to the design performance of 140 $`\mu `$ averaged over the cell. The current status of the $`dE/dx`$ performance is also shown in Fig. 10. The two dark spots at high momentum are the Bhabha electrons and positrons.

Figure 12 shows the $`\gamma \gamma `$ mass spectrum for di-gamma energies above 500 MeV. A clear $`\pi ^0`$ peak is visible—the resolution on the $`\pi ^0`$ is about 30% larger than design.

Figure 12 shows the effect of Cherenkov detector particle identification using the DIRC system. The $`K^\pm \pi ^{}`$ mass spectrum in the region of the $`D^0`$ is shown without using particle identification and using the DIRC to identify the kaon.

Figure 14 shows the $`D^{}`$-$`D^0`$ mass difference, using the two $`D^0`$ decay modes, $`K^\pm \pi ^{}`$ and $`K^\pm \pi ^{}\pi ^0`$ (no particle ID is used). These plots are indicative of good early performance of the tracking systems.

Figure 14 shows the di-lepton invariant mass spectra in the region of the $`J/\psi `$ particle.

While it is much too early to expect high quality (design level) performance from the detector, these plots demonstrate rather impressive performance this early in the life of the detector.

## 5 Future Plans

As of LP99, the accelerator was complete, but BaBar was missing 2/3 of the quartz bars. These bars were installed in October 1999, thereby completing the detector. The facility will run until Summer 2000, with a short break in December 1999.


text:
# On the Mass Profile of Galaxy Cluster CL 0024+1654 Inferred from Strong Lensing

## 1 Introduction

There has been a lot of recent controversy concerning the density profiles of the dark matter halos of virialized cosmological structures, from dwarf galaxy to galaxy cluster scale. N-body simulations of the formation of “virialized” dark matter halos associated with galaxies and clusters in a Cold Dark Matter model were found to be well-fit by a simple, universal form for the variation of mass density $`\rho `$ with radial distance $`r`$ from the center of mass, given by

$$\rho _{\mathrm{NFW}}(r)=\frac{\rho _S}{(r/r_S)(1+r/r_S)^2},$$
(1)
where $`r_S`$ is some characteristic radius which separates the two asymptotic power-law slopes, $`\rho r^1`$ at $`r<<r_S`$ and $`\rho r^3`$ at $`r>>r_S`$, and $`\rho _S`$ is a characteristic density which reflects the mean density of the universe at the epoch of halo formation (Navarro, Frenk, & White 1997; “NFW”). More recent N-body simulations of higher resolution obtain halo profiles which agree with the NFW profile at large radii but have an even steeper inner slope, $`\rho r^{1.5}`$ (Moore et al. 1999).

This prediction by the standard CDM model of singular density profiles for cosmological halos is apparently in conflict with the observed mass distributions inside dark-matter-dominated halos on two extremes of the halo mass function, dwarf galaxies and galaxy clusters. As a result, these observations and their interpretation have recently come under intense scrutiny. In work going back to Flores & Primack (1994), Moore (1994), Burkert (1995) and more recently in Kravtsov et al. (1998), Moore et al. (1999), and McGaugh & de Blok (1999), these predictions of singular halos were found to be in conflict with the halo density profiles derived from the observed rotation curves of dark-matter-dominated dwarf galaxies. The latter were found, instead, to be better fit by density profiles with flat cores. Recently, the universality of this observational requirement that dark-matter-dominated galactic halos possess a uniform-density core has been challenged on the grounds that the rotation curves are not generally resolved well enough in the centers to distinguish the slowly rising rotation curve which results from a mass profile with a flat core from the more rapid rise which results from a cuspy profile van den Bosch et al. (2000). Initially, this challenge appeared to apply primarily to LSB galaxies, leaving the case for uniform-density cores in dwarf galaxies still strong. For example, the well-resolved rotation curves of the nearby dwarf galaxies DDO 154 and NGC 3109 still demanded a flat density core van den Bosch et al. (2000). Very recently, however, this challenge has been extended to the interpretation of dwarf galaxy rotation curves, as well (van den Bosch & Swaters 2000). While the latter are still generally better fit by halo profiles with a flat core, it can no longer be stated with much confidence that the observations demand such profiles. Until better resolution data become available, profiles with $`\rho r^\alpha `$ near the center, with $`0\alpha 1`$, all seem to yield reasonable fits to the data. Hence, while it appears that a significant conflict remains between the CDM halos predicted by the highest resolution N-body simulations to date (which imply $`\rho r^{1.5}`$ at the center) and the dwarf galaxy rotation curves, the current observations are unable to discriminate effectively between halos with inner profiles as cuspy as the NFW profile and those with a uniform-density core, instead.

On the galaxy cluster scale, the case for a uniform-density core has been made most convincingly using observations of strong gravitational lensing, where the images of background galaxies are distorted by the cluster mass to form arcs and multiple images. To date, the most spectacular example of such observation is a Hubble Space Telescope image of multiple arcs produced by the cluster CL 0024+1654 at $`z=0.39`$. The relaxed structure of CL 0024+1654, the absence of a single, central, dominant cluster galaxy, and the presence of an easily-identified multiply-imaged background galaxy make this cluster a particularly good candidate for a determination of the halo mass profile from lensing analysis. According to Tyson, Kochanski, & Dell’Antonio (1998), the observations of this cluster require a halo mass profile with a uniform-density core, in strong conflict with the predicted cuspy profile of NFW.

These conflicts between the CDM N-body results and the observations of dwarf galaxy and cluster halo profiles have stimulated a vigorous re-examination of the theoretical underpinnings of the CDM model, including a number of suggestions for a variation of the microscopic properties of CDM which might serve to produce halos with uniform-density cores while retaining the more successful aspects of the original CDM model. These include suggestions that the dark matter is nongravitationally self-interacting Spergel & Steinhardt (2000), warm (e.g. Sommer-Larsen & Dolgov 2000; Colín, Avila-Reese, & Valenzuela 2000; Hannestad & Scherrer 2000), fluid Peebles (2000), decaying Cen (2000), repulsive Goodman (2000), fuzzy Hu, Barkana, & Gruzinov (2000), and annihilating Kaplinghat, Knox & Turner (2000). In view of the importance of the dwarf galaxy and cluster halo profile observations in constraining the theoretical models, it is somewhat disappointing that the conclusions based upon dwarf galaxy rotation curves are currently so ambiguous with regard to the question of the uniform-density core. This makes the conclusion of Tyson et al. (1998) regarding the core in CL 0024+1654 all the more critical.

Recently, a new study of this cluster by Broadhurst et al. (2000) has reached a conclusion opposite to that of Tyson et al. (1998) regarding the consistency of the observed mass profile with the density profiles predicted by the N-body simulations of cluster formation in the standard CDM model. They find that the NFW mass profile is consistent with the lensing data. If this is correct, then the case against the standard CDM model is significantly weakened. The purpose of this letter is to point out that the fit by Broadhurst et al. (2000) of the cluster lensing data with a mass distribution which follows the NFW profile implies a cluster velocity dispersion which is much larger then the value measured for this cluster by Dressler et al. (1999) of $`\sigma _V=1150\mathrm{km}\mathrm{s}^1`$.

## 2 Observational Results for the Cluster Mass Profile

The projected dark matter density profile which Tyson et al. (1998) found by modeling the lensing data for CL 0024+1654 within the arcs at radius $`r_{\mathrm{arcs}}100h^1\mathrm{kpc}`$ is well fit by

$$\mathrm{\Sigma }(y)=\frac{K(1+\eta y^2)}{(1+y^2)^{2\eta }},$$
(2)
where $`y=r/r_{\mathrm{core}}`$, $`K=7900\pm 100h\mathrm{M}_{}\mathrm{pc}^2`$, $`r_{\mathrm{core}}=35\pm 3h^1\mathrm{kpc}`$, $`\eta =0.57\pm 0.02`$, and $`h`$ is the Hubble constant in units of $`100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. Additionally, Tyson et al. (1998) rule out at a great confidence level the possibility of a good fit of the observed mass distribution by the NFW profile.

A recent paper by Broadhurst et al. (2000) obtained, instead, a good fit to the lensing data with a total mass distribution given by the NFW profile in equation (1), with $`r_S400h^1\mathrm{kpc}`$ and $`\delta _c8000`$, where $`\delta _c\rho _S/\rho _{\mathrm{crit}}(z)`$, and $`\rho _{\mathrm{crit}}(z)3H^2/(8\pi G)`$, the critical density of the universe at the cluster redshift. This $`\delta _c`$ is directly related to the NFW concentration parameter $`cr_{200}/r_S`$, where $`r_{200}`$ is the radius within which the average density is 200 times this critical density, according to

$$\delta _c=\frac{200}{3}\frac{c^3}{[\mathrm{ln}(1+c)c/(1+c)]},$$
(3)
which yields $`c5`$, and $`r_{200}2h^1`$Mpc. From this they concluded that there is no conflict between the observations of this cluster and the predictions of the standard CDM model.

## 3 Consequences for the Cluster Velocity Dispersion

### 3.1 NFW profile

These mass profiles inferred for the cluster CL 0024+1654 based upon lensing observations have implications for the velocity dispersion of its dark matter and galaxies if the cluster is assumed to be in virial equilibrium. We begin by considering the NFW profile. The same N-body simulations of cluster formation in the CDM universe which indicate that clusters are, indeed, in an approximate virial equilibrium which can be described by the universal mass profile suggested by NFW also yield information about the 1-D velocity dispersion $`\sigma _V`$ and its radial dependence. Over a wide range of radii, the halos obtained in N-body simulations are roughly isothermal (Tormen, Bouchet, & White 1997; Eke, Navarro, & Frenk 1998). For comparison with the velocity dispersion of a cluster like CL 0024+1654 observed within some radius, it is necessary to consider the average velocity dispersion of the NFW halo within a sphere of the same radius. The average velocity dispersion of CL 0024+1654 was measured by Dressler et al. (1999) to be $`\sigma _V=1150\mathrm{km}\mathrm{s}^1`$ within a radius $`r600\mathrm{h}^1\mathrm{kpc}6r_{\mathrm{arcs}}`$, based upon 107 galaxy redshifts, to an uncertainty of roughly less than $`\pm 100\mathrm{km}\mathrm{s}^1`$. For the NFW profile proposed by Broadhurst et al. (2000) for this cluster, this radius corresponds to $`rr_{200}/3`$. It is a relatively simple matter to estimate the predicted velocity dispersion $`\sigma _{\mathrm{V},\mathrm{NFW}}`$ for a given NFW density profile in terms of the circular velocity profile $`v_{\mathrm{c},\mathrm{NFW}}`$ of the same halo, as follows.

The circular velocity profile of the NFW halo is given by

$$v_{\mathrm{c},\mathrm{NFW}}^2(r)\frac{GM(r)}{r}=4\pi G\rho _Sr_S^2\frac{[\mathrm{ln}(x+1)x/(x+1)]}{x},$$
(4)
where $`M(r)`$ is the mass enclosed by radius $`r`$, and $`xr/r_S`$. The maximum value of $`v_{\mathrm{c},\mathrm{NFW}}`$ is $`v_{\mathrm{max},\mathrm{NFW}}0.465(4\pi G\rho _Sr_S^2)^{1/2}`$, which occurs at $`x=x_{\mathrm{max}}2.163`$. The NFW velocity profile for the parameters reported by Broadhurst et al. (2000) is shown in Figure 1. We obtain $`v_{\mathrm{max},\mathrm{NFW}}3340\mathrm{km}\mathrm{s}^1`$ for this cluster. According to the detailed analysis of numerical N-body results by Tormen, Bouchet, & White (1997), the average 1-D velocity dispersion within $`rr_{200}/3`$ for simulated clusters, which are well-fit by the NFW profile, is somewhat lower than $`v_{\mathrm{max},\mathrm{NFW}}`$, but never by more than a factor of $`1.5`$. This factor of 1.5 agrees very well with the aperture-averaged, line-of-sight $`\sigma _V`$ which results from solving the Jeans equation for the variation of the radial $`\sigma _{V,r}`$ with $`r`$ inside the NFW profile, including the possible effects of anisotropic velocities (Lokas & Mamon 2000; Lokas, private communication). Therefore, the NFW profile proposed for CL 0024+1654 by Broadhurst et al. (2000) implies an average velocity dispersion for the cluster within the radius $`r600\mathrm{h}^1\mathrm{kpc}6r_{\mathrm{arcs}}`$ of $`\sigma _{\mathrm{V},\mathrm{NFW}}>v_{\mathrm{max},\mathrm{NFW}}/1.52230\mathrm{km}\mathrm{s}^1`$, which is too large by a factor of $`2`$ to be compatible with the measured velocity dispersion quoted above 111This assumes $`\rho _S0.006h^2\mathrm{M}_{}/\mathrm{pc}^3`$, the value if $`\rho _{\mathrm{crit}}(z)`$ is for an Einstein-de Sitter universe. Our conclusion that $`\sigma _{V,NFW}`$ substantially exceeds the observed $`\sigma _\mathrm{V}`$ remains unchanged if, instead, we consider other background cosmologies. The values of $`v_{\mathrm{max},\mathrm{NFW}}`$ and $`\sigma _{\mathrm{V},\mathrm{NFW}}`$ for a low-density universe ($`\mathrm{\Omega }_0=0.3`$) are only slightly lower than the values quoted above, by 10% for the open case ($`\mathrm{\Lambda }=0`$) and by 25% for the flat ($`\mathrm{\Lambda }0`$) case..

### 3.2 Nonsingular Isothermal Sphere

The Tyson et al. (1998) mass model in equation (2) does not uniquely predict the velocity dispersion $`\sigma _V`$ implied for the cluster even if it is assumed to be in virial equilibrium. For one, it is a fit only  to the projected density within the arcs at $`r_{\mathrm{arcs}}100\mathrm{h}^1\mathrm{kpc}`$ and, therefore, does not constrain the external mass profile very well. For another, even if we restrict our attention to radii within the arcs and assume virial equilibrium, this is still not enough to specify $`\sigma _V`$ and its radial dependence there uniquely. As a result, it is not possible to compare the predicted $`\sigma _V`$ for the Tyson et al. (1998) mass model with the observed $`\sigma _V`$ as directly as we did for the NFW profile above. However, we can at least demonstrate that the mass model is plausibly consistent with the observed $`\sigma _V`$, as follows.

The projected mass density in equation (2) is well fit by that obtained from the truncated, nonsingular, isothermal sphere (“TIS”) profile of Shapiro, Iliev, & Raga (1999). We discuss the details of the fitting procedure and the results elsewhere (Iliev & Shapiro 2000a,b). This TIS represents a particular solution of the Lane-Emden equation which corresponds to the outcome of the collapse and virialization of a top-hat density perturbation. The size and virial temperature which result are unique functions of the mass and redshift of formation of the object for a given background universe. According to this solution, the central density $`\rho _0`$ is roughly more than 500 times the density at the surface, and the core radius $`r_0`$ is about 1/30 of the total size. \[Note: Our definition of the core radius is $`r_{0,\mathrm{TIS}}r_{\mathrm{King}}/3`$, where $`r_{\mathrm{King}}`$ is the “King radius” defined in Binney & Tremaine (1987), p. 228.\] As described elsewhere, this solution is a convenient analytical approximation for the halos which form from more realistic initial conditions in the CDM model, which reproduces many of the average structural properties of the halos found in CDM simulations, except in the very inner profile where the TIS model has a uniform-density core instead of a central cusp. As such, a fit of the TIS profile to the Tyson et al. (1998) mass model provides a plausible, physically-motivated connection between this mass model with a flat density core and the implied cluster velocity dispersion. The projected density profile of a TIS sphere with a central density $`\rho _00.064h^2\mathrm{M}_{}/\mathrm{pc}^3`$ and a core radius $`r_020h^1`$kpc provides a very close match to the Tyson et al. (1998) result for $`\mathrm{\Sigma }(r)`$ discussed above. Based on this best-fit TIS, the velocity dispersion is

$$\sigma _{\mathrm{V},\mathrm{TIS}}=(4\pi G\rho _0r_0^2)^{1/2}1200\mathrm{km}\mathrm{s}^1,$$
(5)
in close agreement with the measured value.

## 4 Conclusion

Two different attempts to invert the observations of strong lensing by cluster CL 0024+1654 to solve for the mass profile of the cluster have reached opposite conclusions regarding the presence of a uniform-density core versus the acceptability of a central cusp like that of the NFW profile. This suggests that, either there is some error in one or both of these analyses or else there is some degeneracy in the inversion process which prevents a clean discrimination between these different solutions. While we do not claim to address the accuracy of either of the two analyses of this cluster, we point out, instead, that there may be additional constraints on the allowed mass models which can aid in distinguishing them. In particular, the measured cluster velocity dispersion, which also reflects the mass distribution of the cluster, should also be required to be consistent with the mass model derived from the inversion of the lensing data. We have shown that this is not the case for the cuspy NFW profile which Broadhurst et al. (2000) report is a good fit to the lensing data. That profile, if it actually corresponds to a halo formed in the standard CDM model, predicts much too high an average velocity dispersion to be consistent with the observed value. By contrast, the mass model of Tyson et al. (1998) does not have this problem. This suggests that the lensing data for CL 0024+1654 still favor a flattening of the density profile at small radii, in conflict with the prediction of cuspy halos by N-body simulations of the standard CDM model.

In view of the importance of this conclusion for the ongoing debate regarding the validity of the standard CDM model, it would be valuable if uncertainties in the gravitational lense models, such as that due to possible departures from spherical symmetry or substructure, were properly quantified. This uncertainty would have to be rather extreme, however, to reconcile the Broadhurst et al. (2000) profile fit with the observed velocity dispersion of CL 0024+1654. The lensing analysis leads to NFW parameters roughly by measuring the mass interior to the arcs, $`M(r_{\mathrm{arcs}})=4\pi \rho _Sr_S^3f(x_{\mathrm{arcs}})`$, where $`f(x)=\mathrm{ln}(1+x)x(1+x)^1`$ and $`xr/r_S`$, and determining $`r_S`$ by matching the logarithmic slope $`\gamma `$ of the density profile at $`r_{\mathrm{arcs}}`$ as derived from the projected mass distribution, according to $`x_{\mathrm{arcs}}=(1+\gamma )/(3+\gamma )`$. In order to adjust $`v_{\mathrm{max},\mathrm{NFW}}`$ downward relative to the value discussed above by a factor large enough to reconcile the NFW profile with the observed velocity dispersion while leaving the measured mass $`M(r_{\mathrm{arcs}})`$ unchanged, therefore, $`r_S`$ must be reduced by factor of order 4, to $`r_S100h^1\mathrm{kpc}`$. This requires an observational uncertainty so large as to change the value $`\gamma 1.3`$ reported by Broadhurst et al. (2000) into $`\gamma =2`$, which seems very unlikely.

We are grateful to Alan Dressler for valuable discussion and to the referee for comments which led us to improve the paper. This work benefited from the partial support of grants NASA ATP NAG5-7821, NSF ASC-9504046, and Texas Advanced Research Program 3658-0624-1999.


text:
# Untitled Document

Characteristic features of anharmonic effects in the lattice dynamics of fcc metals

M. I. Katsnelson1, A. V. Trefilov2, and K. Yu. Khromov2.

1Institute of Metal Physics, Urals Branch of the Russian Academy of Sciences, 620219 Ekaterinburg, Russia

2 Russian Research Center ”Kurchatov Institute”, 123182 Moscow, Russia

Abstract
The dispersion in the entire Brillouin zone and the temperature dependence (right up to the melting temperature) of the anharmonic frequency shift and phonon damping in a number of fcc metals is investigated on the basis of microscopic calculations. It is found that the anharmonic effects depend sharply on the wave vector in the directions $`\mathrm{\Gamma }`$-X, X-W, and $`\mathrm{\Gamma }`$-L and, in contrast to bcc metals, the magnitude of the effects is not due to the softness of the initial phonon spectrum. It is shown that the relative frequency shifts and the phonon damping near melting do not exceed 10-20%. The relative role of various anharmonic processes is examined, and the relation between the results obtained and existing experimental data is discussed.

The investigation of anharmonic effects (AEs) in lattice dynamics is a classic problem of solid-state physics. It is important, specifically, because of the role that such effects they can play in phenomena associated with structural phase transitions and melting in crystals (see, for example, Refs. 1-3). At the same time, obtaining any information about the magnitude and scale of AEs from experiment and theory is a difficult problem. The experimental study of such ”basic” AEs as the frequency shift and damping of phonons is very difficult and leads to a large uncertainty in the results (see, for example, the data presented in Refs. 4 and 5 for bcc and fcc metals, respectively). Up to now first-principles microscopic calculations of AEs have been performed for one point of the Brillouin zone (N) in the bcc phase of Zr and four points (N, P, $`\omega `$, G) in Mo. (Ref. 6) Detailed information about AEs in the entire Brillouin zone and their temperature dependence has been obtained in Refs. 4 and 7 on the basis of pseudopotential theory for the bcc phases of alkali and alkaline-earth metals. For these metals the most striking manifestations of AEs are due to the ”soft-mode behavior” (the anomalous temperature dependence of the phonon frequencies) of the $`\mathrm{\Sigma }_4`$ branch. It is of interest to calculate AEs for the ”general position”, i.e., for crystals which do not possess soft vibrational modes. Such crystals include most metals with close-packed structures, for example, fcc. There is virtually no information about AEs in such metals, not counting Ref. 5, where the damping of phonons in precious metals and Al was calculated. However, the calculations were performed only for individual points of the Brillouin zone and at room temperature (which is much less than the melting temperature $`T_m`$). Moreover, the accuracy of the model used in Ref. 5 for interatomic interactions is unclear, which is probably the reason for the qualitative discrepancy between the theoretical and experimental dependences of the damping on the wave vector. In the present letter the AEs in the lattice dynamics of fcc metals are investigated in detail on the basis of a systematic anharmonic perturbation theory.

The calculations were performed for Ir, the fcc phases of Ca and Sr, and the hypothetical fcc phase of K. This choice of objects was determined by the fact that reliable and at the same simple microscopic models which make it possible to describe a wide range of lattice properties of these metals exist Refs. 4,7-9. Iridium is an example of an fcc crystal with a ”stiff interatomic interaction potential, similar to the Lennard-Jones potential, Ref. 8 while K, Ca, and Sr are characterized by potentials with a softer ”core”. To demonstrate the characteristic features of AEs due to the specific nature of the potential with the same lattice geometry, we present here the results for the ”limiting” cases — Ir and K. The results for Ca and Sr are qualitatively similar to the results for K and will be presented elsewhere.

The parameters determining the interatomic interaction potential are presented in Ref. 8 for Ir and in Ref. 4 for K. The calculations were performed on the basis of the anharmonic perturbation theory taking account of thermal expansion (quasiharmonic contributions, qh) and three- and four-phonon interaction processes to accuracy $`\varkappa ^2`$, where $`\varkappa =(m/M)^{1/4}`$ is the adiabatic smallness parameter, and m and M are the electron and ion masses. The exact formulas actually used in the calculations are presented in Ref. 4. For a qualitative discussion of the results, we present here their high-temperature asymptotic representation for $`T>\mathrm{\Theta }_D`$ ($`\mathrm{\Theta }_D`$ is the Debye temperature):

$$\mathrm{\Delta }_{\lambda 𝐤}=\mathrm{\Delta }_{\lambda 𝐤}^{qh}+\mathrm{\Delta }_{\lambda 𝐤}^3+\mathrm{\Delta }_{\lambda 𝐤}^4$$
(1)
$$\mathrm{\Delta }_{\lambda 𝐤}^{qh}=\gamma _{\lambda 𝐤}\omega _{\lambda 𝐤}\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Omega }},$$
(2)
$$\mathrm{\Delta }_{\lambda 𝐤}^3=\frac{T}{4M^3\omega _{\lambda 𝐤}}\underset{\mu \nu 𝐪}{}\left|V_{\lambda \mu \nu }^{𝐤,𝐪,𝐤+𝐪}\right|^2\frac{1}{\omega _1^2\omega _2^2}\frac{(\omega _1^2\omega _2^2)^2\omega _{\lambda 𝐤}^2(\omega _1^2+\omega _2^2)}{(\omega _1^2+\omega _2^2\omega _{\lambda 𝐤}^2)^24\omega _1^2\omega _2^2},$$
(3)
$$\mathrm{\Delta }_{\lambda 𝐤}^4=\frac{T}{4M^3\omega _{\lambda 𝐤}}\underset{\mu 𝐪}{}W_{\lambda \lambda \mu \mu }^{𝐤,𝐤,𝐪,𝐪}\frac{1}{\omega _{\mu 𝐪}^2},$$
(4)
$$\mathrm{\Gamma }_{\lambda 𝐤}=T\frac{\pi }{8M^3}\underset{\mu \nu 𝐪}{}\left|V_{\lambda \mu \nu }^{𝐤,𝐪,𝐤+𝐪}\right|^2\frac{1}{\omega _1^2\omega _2^2}[\delta (\omega _{\lambda 𝐤}\omega _1\omega _2)+2\delta (\omega _{\lambda 𝐤}+\omega _1\omega _2)]$$
(5)
Here $`\lambda `$, $`\mu `$, and $`\nu `$ are the indices of the phonon branches, k and q are quasimomenta, $`\omega _{\lambda 𝐤}`$, $`\mathrm{\Delta }_{\lambda 𝐤}`$, and $`\mathrm{\Gamma }_{\lambda 𝐤}`$ are, respectively, the initial phonon frequency and its shift and the phonon damping, $`\widehat{V}`$ and $`\widehat{W}`$ are the amplitudes of three- and four-phonon processes (see Ref. 4), $`\gamma _{\lambda 𝐤}=\omega _{\lambda 𝐤}/\mathrm{\Omega }`$ is Gruneisen parameter. $`\mathrm{\Delta }\mathrm{\Omega }`$ is the change in the volume $`\mathrm{\Omega }`$ per atom due to thermal expansion, and $`1(\mu ,𝐪)`$ and $`2(\nu ,𝐤+𝐪)`$.

The basic computational results are presented in Figs. 1-5. It is evident from Figs. 1 and 2 that the high-temperature asymptotic representations (3)-(5) obtain very early (for $`T\mathrm{\Theta }_D/3`$) Comparing the analogous results for bcc metals4 it can be concluded that this property is probably quite general. It follows from Figs. 1-4 that the strongest AEs in Ir occur for longitudinal phonons near the point X. Specifically, the decrease of the frequency and the increase of the damping with temperature are greatest for these phonons The situation differs sharply in this respect from bcc alkali and alkaline-earth metals, Refs. 4,7 where the maximum damping obtains for ”soft” phonons, for which the opposite behavior of the frequency with increasing $`T`$ is characteristic — hardening instead of softening. It is evident from Fig. 5 that this difference is due to the fact that in Ir three-phonon processes dominate over four-phonon processes.

To understand the role of the characteristic features of the interatomic interactions for the same lattice geometry, the computational results tor a hypothetical fcc phase of K are presented in Figs. 3 and 4 The explicit domination of three-phonon over four phonon processes, which causes the frequency shift to be negative over the entire Brillouin zone is also seen in this case and is apparently characteristic for all metals with fcc structure At the same time, the fact that the neighborhood of the point X is distinguished probably occurs only for crystals with a ”hard” interatomic interaction potential of the type in Ir.

It is well known Ref. 8 that the ”hardness” of the potential in metals increases with the effective valence Z (Z= 1 for K and Z=4.5 for Ir). In any case Ir is a striking example of a crystal where AEs are maximum for high-frequency vibrations, in contrast to crystals with soft modes where the AEs are maximum for low-frequency phonons. The latter can be easily explained by the high powers of the frequency in the denominators of the Eqs. (3)-(5), but the example of Ir shows that in the general case the question cannot be solved on the basis of such simple considerations.

According to Fig. 4, an important feature of phonon damping is its nonmonotonic dependence on the wave vector in the $`\mathrm{\Gamma }L`$ and $`\mathrm{\Gamma }X`$ directions. According to the experimental data presented in Ref. 5, such nonmonotonic behavior is observed in Al, Cu, Ag, and Au. It can therefore be supposed that this behavior is typical for all metals with the fcc structure.

In closing we note that, as follows from Figs. 3-5, the AEs have a strong dependence over the Brillouin zone. For this reason, approximations which take them into account ”on average” (approximations of the type used in the recently published Ref. 10), are dangerous. Finally, from the fact, found in the present work, that three-phonon processes predominate over four-phonon processes it follows that the well-known approaches such as the self-consistent phonon approach are inapplicable to metals with a close-packed structure Ref. 11.

references

1. V. G. Vaks, Introduction to the Microscopic Theory of Ferroelectrics (Moscow, Nauka, 1973).

2. V. G. Vaks, S. P. Kravchuk, and A V Trefilov, J. Phys. 10(F) (1980), 2325.

3. M. I. Katsnel’son and A. V. Trefilov, Fiz. Met. Metalloved. 64 (1987), 629.

4. V. G. Vaks, S. P. Kravchuk, and A. V. Trefilov, J. Phys., 10(F) (1980), 2105.

5. M. Zoli, G. Santoro, V. Bortolani et al., Phys. Rev., 41(B) (1990), 7507.

6. Y. Y. Ye, Y. Chen, K. M. Ho el al, Phys. Rev. Lett. 58 (1987), 1769; J. Phys.: Condens. Matter 3 (1991), 9629.

7. V. G. Vaks, G. D. Samolyuk, and A. V. Trefilov, ”On anomalies of anharmonic effects for soft phonons in alkali and bcc alkaline-earth metals”, Phys. Lett., 127(A) (1988), 37.

8. A. S. Ivanov, M. I. Katsnelson, A. G. Mikhin et al., ”Phonon spectra, interatomic interaction potentials and simulation of lattice defects in iridium nad rhodium”, Philos. Mag., 69(B) (1994), 1183.

9. M. I. Katsnelson, I. I. Naumov, A. V Trefilov et al., ”Peculiarities of phonon spectra and lattice heat capacity in Ir and Rh”, Philos. Mag., 75(B) (1997), 389.

10 R. C. Shukla and E. R. Cowley, Phys. Rev., 58(B) (1998), 2596.

11. H. Bottger, Principles of the Theory of Lattice Dynamics (Weinheim, Physik-Verlag, 1983; Moscow, Mir, 1986).

Figure captions

1. Temperature dependence of the frequency shift in Ir at the symmetric points of the Brillouin zone; $`T_m`$ — melting temperature, $`\mathrm{\Theta }_D`$ — Debye temperature. The indices 1 and 2 denote longitudinal and transverse phonons; $`\omega _{pl}`$ — ion plasma frequency ($`\omega _{pl}^{Ir}=871`$ K, see Ref. 8).
2. Temperature dependence of phonon damping in Ir. All notations are the same as in Fig 1.
3. Dispersion of the relative frequency shift $`\delta =\mathrm{\Delta }/\omega `$ at $`T=T_m`$. Solid line - Ir, dotted line - fcc K. The numbers 1, 2 indicate longitudinal and transverse branches.
4. Dispersion of the relative phonon damping $`\eta =\mathrm{\Gamma }/\omega `$ at $`T=T_m`$. Solid line - Ir, dotted line - fcc K. The numbers 1, 2 indicate longitudinal and transverse branches.
5. Variation of the anharmonic contributions $`\mathrm{\Delta }^{qh}`$ (a), $`\mathrm{\Delta }^3`$ (b), and $`\mathrm{\Delta }^4`$ (c) in the $`\mathrm{\Gamma }XW`$ directions in Ir at $`T=T_m`$.


text:
# Optical realization of universal quantum cloning

## I Introduction

A deterministic and perfect quantum cloner for arbitrary quantum states is forbidden by the quantum no-cloning theorem. As an extension, Yuen and D’Ariano showed that a violation of unitarity makes the cloning of two nonorthogonal states also impossible. Then the case of mixed states is considered and it is shown that two noncommuting mixed states cannot be broadcast, even when the state need only be reproduced marginally on the separate systems. For further consideration, Koashi and Imoto generalized the standard no-cloning theorem to the entangled states. Their research showed to what extent the correlation to other systems can be read out from a subsystem without altering its marginal density operators.

After having considered about so many restrictions posed on cloning quantum states by quantum mechanics, people began to think about how to do our best in quantum states cloning under all these restrictions, that is, the inaccurate cloning of quantum states. These approached may be divided into two main categories: universal and state-dependent.

State-dependent quantum cloning machines are designed to clone states belonging to a finite set. It can be divided into three subcategories: deterministic, probabilistic exact and hybrid.

On the other hand, universal quantum cloning machines deterministically generates approximate copies of an unknown quantum state. The word “universal” means that this kind of machines works equally well for any input states, i.e., state-independent. This category can again be divided into two subcategories: symmetric and asymmetric.

Asymmetric universal quantum cloning machines which introduced by Cerf and Bužek et al produce two non-identical imperfect copies of a single state of $`N`$-dimensional quantum system. Symmetric universal quantum cloning machines, first designed by Bužek and Hillery, works in a different way: act on any unknown quantum state and produce identical approximate copies equally well. This kind of Bužek-Hillery cloning machine has been optimized and generalized.

Quantum cloning may be used to improve the performance of some quantum computation tasks and may be heleful in understanding quantum information theory and processing quantum imformation, and might be helpful in eavesdropping versus quantum cryptography. And we can improve the performance of the measurements on observables by performing measurements on clones of the original quantum system.

In the present paper, we realized the one-to-two qubits Bužek-Hillery cloning machine with linear optical devices. This method relies on the representation of several qubits by a single photon. We showed that, the fidelities between the two output qubits and the original qubit are both $`\frac{5}{6}`$ (which proved to be the optimal fidelity of one-to-two qubits universal cloner) for arbitrary input pure states.

## II Universal cloning

The most important quality of symmetric universal cloning machines is how well we can do, i.e., the optimal fidelity. We use the most usual definition of fidelity $`F=\mathrm{\Psi }\left|\rho \right|\mathrm{\Psi }`$. For convenience, we only discuss the case when the input states are pure. Here,$`|\mathrm{\Psi }`$ is the to-be-cloned unknown pure state and $`\rho `$ is the density matrix of one of the output copies.

A $`MN`$ symmetric universal quantum cloning machine is described as the following (for qubit systems):

1. The input system contains $`M`$ original qubits, each in the same unknown pure state $`|\mathrm{\Psi }`$, $`NM`$ “blank” qubits $`|0`$ and $`K`$ auxiliary qubits $`|0`$.

2. The cloning machine acts on the input system and performs a prescribed unitary transformation $`U`$ on it.

3. The output system contains $`N`$ identical copies (usually entangled), each copy described by a reduced density matrix $`\rho _i^{out}`$ $`(i=1,\mathrm{},N)`$, and $`K`$ auxiliary qubits.

The cloning process can be written as

$$U|\mathrm{\Psi }^M|0^{(NM)}|aux=|\mathrm{\Phi }^{out}\text{.}$$
(1)
For $`|\mathrm{\Phi }^{out}`$, there exists $`N`$ qubits such that each of them has the same reduced density matrix $`\rho ^{out}`$.

For the $`MN`$ symmetric universal quantum cloning machine described above, the optimal fidelity has been demonstrated as

$$F^{opt}(M,N)=\frac{MN+N+M}{N(M+2)}\text{.}$$
(2)
From the equation we can see that the optimal $`12`$ fidelity is $`\frac{5}{6}`$.

The special unitary transformation $`U`$ of the one-to-two Bužek-Hillery cloning machine can be written as

$`|0_1|0_2|0_3\stackrel{U}{}\sqrt{{\displaystyle \frac{2}{3}}}|00_{12}|0_3+\sqrt{{\displaystyle \frac{1}{3}}}|+_{12}|1_3\text{,}`$ ()
$`|1_1|0_2|0_3\stackrel{U}{}\sqrt{{\displaystyle \frac{2}{3}}}|11_{12}|1_3+\sqrt{{\displaystyle \frac{1}{3}}}|+_{12}|0_3\text{,}`$ (3)
where

$$|+=\sqrt{\frac{1}{2}}(|10_{12}+|01_{12})$$
(4)
In the cloning process, the quantum information in the original qubit (the first qubit) is copied to qubits $`1`$ and $`2`$ equally well such that the reduced density matrixes at the output are identical, i.e., $`\rho _1^{out}=\rho _2^{out}`$. The third qubit acts as an auxiliary qubit here.

From the specific expression of the unitary transformation $`U`$ given above we can easily calculate the fidelities

$`F_1`$ $`=`$ $`\mathrm{\Psi }|\rho _1^{out}|\mathrm{\Psi }={\displaystyle \frac{5}{6}}\text{,}`$ ()
$`F_2`$ $`=`$ $`\mathrm{\Psi }|\rho _2^{out}|\mathrm{\Psi }={\displaystyle \frac{5}{6}}\text{,}`$ (4)
where $`|\mathrm{\Psi }`$ is an arbitrary unknown pure state

$$|\mathrm{\Psi }=\alpha |0+\beta e^{i\phi }|1,\text{ }\alpha ,\beta real,\text{ }\phi [0,2\pi ).$$
(5)
So, it is shown to be an optimal symmetric universal cloning machine.

In order to realize the transformation $`U`$ in experiment, we first have to construct the quantum logic circuit. To do that, $`U`$ is decomposed into a sequence of basic operations such as the rotation of a single qubit and the controlled-NOT (CNOT) operation of two qubits. The final network of the cloning machine has been shown in figure 1.

We note that the network is divided into two parts (except the state swapping part): preparation and cloning. In the preparation stage, qubit 1 is left unaffected and qubit 2 and 3 are entangled by CNOT gates. That is, there is no “flow” of quantum information from qubit 1 to qubit 2, 3. Then, in the cloning stage, the quantum information in qubit 1 is “redistributed” among the three qubits by a sequence of four CNOT gates.

The unitary transformations in the two stages are described as

$$|\mathrm{\Psi }_1|0_2|0_3\stackrel{U_1}{}|\mathrm{\Psi }|\mathrm{\Psi }_{23}^{prep}\stackrel{U_2}{}|\mathrm{\Psi }_{123}^{out}\text{,}$$
(5)
where $`|\mathrm{\Psi }_{23}^{prep}`$is prepared by three ratations and two CNOT operations

$$|\mathrm{\Psi }_{23}^{prep}=\widehat{R}_2(\theta _3)\widehat{P}_{32}\widehat{R}_3(\theta _2)\widehat{P}_{23}\widehat{R}_2(\theta _1)|0_2|0_3.$$
(6)
The rotation angles $`\theta _j`$ $`(j=1,2,3)`$ are determined by the specific $`|\mathrm{\Psi }_{23}^{prep}`$. In this case the state $`|\mathrm{\Psi }_{23}^{prep}`$ reads

$$|\mathrm{\Psi }_{23}^{prep}=\sqrt{\frac{1}{6}}(2|00_{23}+|01_{23}+|11_{23}).$$
(7)
Although there is no flow of information in the preparation stage, we note that we can controll the distribution of information among the three qubits by the choice of the ratation angles. We can copy qubit 1 to qubits 1 and 2 or qubits 2 and 3, even copy qubit 1 to qubits 1, 2 and 3 when the input state $`|\mathrm{\Psi }`$ is

$$|\mathrm{\Psi }=\alpha |0+\beta |1,\text{ }\alpha ,\beta real.$$
(8)
This is called quantum triplicator.

## III Experiment description

Now we concentrate on constructing the network given in figure 1 with linear optical devices. The basic idea is the representation of several qubits by a single photon. In such a scheme, qubits are represented by different freedom degrees of a single photon, i.e., a photon can carry one polarization qubit and several location qubits. Two orthogonal linear polarization states of the photon, i.e., horizontal and vertical, serve as the basis states of the polarization qubit. Hereafter, they are denoted by $`|H`$ and $`|V`$ corresponding to $`|0`$ and $`|1`$, respectively. Any of the location qubits is characterized by the information about “which path” is taken by the photon when passing a beam splitter. In this way, entanglement of different qubits can be accomplished conveniently, and decoherence is very low because the photons have relatively less interaction with the environment.

In our scheme of universal quantum cloning, three qubits are involved: one polarization (qubit $`1`$) and two location qubits (qubits $`2`$ and $`3`$). Qubit $`1`$ is the initial to-be-cloned qubit in an arbitrary state $`|\mathrm{\Psi }_1`$, which is cloned to qubits $`1`$ and $`2`$ at the output end with qubit $`3`$ serving as an ancilla. Both location qubits are initially in blank states: $`|0_2`$ and $`|0_3`$.

To realize the network in figure 1, we must construct CNOT gates between any two qubits and unitary rotation of any single qubit. With the representation of qubits mentioned above, all these gates can be made with linear devices. However, in previous schemes, a unitary rotation to a single location qubit is realized by a shiftable beam splitter (BS) or a Mach-Zehnder (M-Z) interferometer. Unfortunately, a shiftable BS is infeasible and an M-Z interferometer will inevitably multiply the complexity.

To overcome this difficulty, a trick of state-swapping between qubits $`1`$ and $`2`$ is employed to transform the input state $`|\mathrm{\Psi }_1|0_2|0_3`$ to $`|H_1|\mathrm{\Psi }_2|0_3`$. Then with a definite polarization state, unitary rotations to a single location qubit can be realized by an operation of Controlled-Rotation (polarization controls location), which can be implemented by only using one Polarization Beam Splitter (PBS) and several properly oriented Half-Wavelength Plates (HWP).

After the state-swapping, we can begin to construct the two parts of our network: preparation and cloning. After having finished this work, we may have a complicated construction (composed of many universal logic gates) of such a cloning machine with optical deveces. But by techniques of compiling, the network can be significantly simplified and a feasible setup is shown in figure 1.

Now let’s explain the experiment setup (figure 1) in detail.

The performance of this cloning machine is demonstrated by using a laser beam attenuated to the single-photon level as input source. At the output, the photons are detected using geigermode avalanche photodiodes — Single Photon Counting Modules ($`D1`$-$`D8`$, EG&G #SPCM-AQ, efficiency $`70\%`$, dark count $`50s^1`$). In our experiment, a well polarized light beam is created by a He-Ne laser with a Brewster window (Melles Griot, $`633`$ $`nm`$ ) and its polarization degree is improved by a polarizer. Then a PBS follows, and its vertical output is measured by a power meter $`D_0`$ to confirm the stability of the laser’s power during the experiment. The horizontal output of this PBS is greatly attenuated (by an attenuating plate, ATTN) so that the maximum detection rates are always less than $`20000s^1`$; for the setup passage time of $`5`$ $`ns`$, this means that on average fewer than $`10^4`$ photons are in the setup at any time. Now the input source can be viewed as single-photon source.

After that the to-be-cloned state is prepared by a properly oriented HWP and an adjustable waveplate (AJWP, namely, a Pockels Cell).

To confirm the correct result, in a general way, quantum state tomography at the output end is required, and a simple version by optical interference is considered in the present system. Another auxiliary qubit is utilized to determine the reduced density operators of the two replicas respectively. This auxiliary qubit (qubit aux) is introduced as $`\frac{1}{\sqrt{2}}(|0_{aux}+|1_{aux})`$ by four $`5050`$ BSes. Then a controlled state swapping is performed with qubit aux as the controller: if qubit aux is $`|1_{aux}`$, the states of qubits $`1`$ and $`2`$ interchange, or else they remain unchanged. The four paths with qubit aux as $`|0_{aux}`$ determine the density operator of qubit $`1`$ (polarization): if the relative photon counts and the polarization state of path $`i`$ is $`C_i`$ and $`\rho _i`$ respectively (by measuring four linearly independent states of the polarization qubit: $`|H,`$ $`|V,`$ $`|D,`$ $`|R`$, we can determine each component of the matrix $`\rho _i`$, here $`|D=\frac{1}{\sqrt{2}}(|0+|1),`$ $`|R=\frac{1}{\sqrt{2}}(|0+e^{i\frac{\pi }{2}}|1)`$,and this method is called quantum tomography), the density operator of qubit $`1`$ can be written as $`C_i\rho _i`$. Other four paths determine the density operator of qubit $`2`$ (location) which has been swapped to the polarization qubit with aux as $`|1_{auz}`$, and the method of measurement is the same as above.

Then the reduced density operators of the replicas provide the evaluation of the fidelities,which are shown in Figure 2. We can see that the experiment datas coincide with the theoretical value $`\frac{5}{6}`$ quite well.

Before qubit $`2`$ is under verification, it has been swapped with the polarization one, so errors of both qubits’ fidelities can be analyzed in similar ways. Now we will analyze the error of our experiment. The error of the fidelity $`F`$ originates mainly from two ways:

1. The oscillation of the photon counts of pass $`i`$ (denoted by $`\mathrm{\Delta }C_i`$), which is primarily caused by the small oscillation of the reflexivity and transmissivity of the first HWP (in the stage of state preparation) for different polarization states.

2. The precision limits of the waveplates and polarizers’ orientation angles (denoted by $`\mathrm{\Delta }\theta `$).

It can be deduced that an upper bound of the maximum error of $`F`$ can be expressed as $`\mathrm{\Delta }F=\underset{i=1}{\overset{4}{}}\mathrm{\Delta }C_i+\frac{3}{2}\mathrm{\Delta }\theta `$. In the present system, $`\mathrm{\Delta }\theta 0.0018`$ $`rad`$ ($`0.1^{}`$) and $`\underset{i=1}{\overset{4}{}}\mathrm{\Delta }C_i0.002`$. Thus we get the maximum error of $`F`$ as $`\mathrm{\Delta }F0.005`$, which has been embodied in the result figure as error bars.

A more complicated network can be designed to carry out the $`NM`$ cloning, whose principle has been generalized from the $`12`$ cloning. But in the present scheme, qubits are represented by different freedom degrees of a single photon, and about $`2^n`$ optical elements are needed to represent $`n`$ qubits. So the complexity of the network will increase exponentially with the increase of $`N`$ and $`M`$, and this is a shortcoming of the present method.

In the present experiment, qubits are represented by different freedom degrees of a single photon. With a quantum non-demolition (QND) measurement, the information of the location qubit may be extracted without destroying the polarization qubit, in principle.

## IV Setup and Results

Figure 1 is the setup of the universal quantum cloning machine. The inset is the logic circuit ($``$ and $``$ denote the controller and target of a CNOT operation respectively, $`R`$ is a rotation to a single qubit). We represent the three input qubits by three freedom degrees of a single photon. Qubit $`1`$ (to-be-cloned polarization qubit) is introduced by setting the photons’ initial polarization in an arbitrary state $`|\mathrm{\Psi }(\theta ,\delta )_1=\mathrm{cos}\theta |H_1+e^{i\delta }\mathrm{sin}\theta |V_1`$ $`\left(\frac{\pi }{2}<\theta \frac{\pi }{2},\text{ }0\delta <\pi \right)`$. Qubits $`2`$ and $`3`$ are both location qubits, which are introduced by different choices of “beam paths” with the initial state $`|00_{2,3}`$. By the state-swapping, qubit $`2`$ is transformed to $`|\mathrm{\Psi }(\theta ,\delta )_2=\mathrm{cos}\theta |0_2+e^{i\delta }\mathrm{sin}\theta |1_2`$, while qubit $`1`$ is transformed to $`|H_1`$. In the stage of entanglement preparation, qubits $`1`$ and $`3`$ are prepared into $`|\mathrm{\Psi }_{1,3}^{prep}=\frac{1}{\sqrt{6}}\left(2|H0_{1,3}+|H1_{1,3}+|V1_{1,3}\right)`$. Then several CNOT operations follow to accomplish the cloning, which can be compiled to two independent polarization M-Z interferometers. Quantum state tomography is adopted to confirm the correct result. Another auxiliary qubit enables us to measure the density operators of the replicas separately, which is introduced as $`\frac{1}{\sqrt{2}}\left(|0_{aux}+|1_{aux}\right)`$ by arranging four $`50`$-$`50`$ BSes in each potential paths of the output photon. Then a controlled state swapping is carried out and the four qubits (qubit $`1,`$ $`2,`$ $`3`$ and $`aux`$) evolve to the state $`|\mathrm{\Psi }^{meas}=`$ $`\frac{1}{\sqrt{2}}\left(|0_{aux}|\mathrm{\Psi }_{1,2,3}^{out}+|1_{aux}|\mathrm{\Psi }_{2,1,3}^{out}\right)`$. When the qubit $`aux`$ is $`|0_{aux}`$, the density operator of replica qubit $`1`$ (polarization) can be obtained by detectors $`1`$ to $`4`$ ($`D_1`$-$`D_4`$) and corresponding QUWPs and polarizers, by measuring the mixed state of polarization and relative photon counts of each path. To the other four outputs of the BSes where the qubit $`aux`$ is $`|1_{aux}`$, a state-swapping between qubits $`1`$ and $`2`$ is performed so that the density operator of replica qubit $`2`$ is exactly realized by knowing that of the interchanged polarization qubit, which is measured by $`D_5`$-$`D_8`$ and corresponding QUWPs and polarizers.

Figure 2 is the replicas’ fidelities in the optimal universal cloning machine. The fidelity is defined as $`\mathrm{\Psi }(\theta ,\delta )\left|\rho \right|\mathrm{\Psi }(\theta ,\delta )`$ with $`\rho `$ denoting the density operator of a replica. This machine can deal with any state $`|\mathrm{\Psi }(\theta ,\delta )=\mathrm{cos}\theta |V+e^{i\delta }\mathrm{sin}\theta |H`$ equally well. From down to up, the four subfigures represent the cases $`\delta =0,\frac{\pi }{4},\frac{\pi }{2},\frac{3\pi }{4}`$, respectively. In each subfigure, the state angle $`\theta `$ of the to-be-cloned qubit varies from $`\frac{\pi }{2}`$ to $`\frac{\pi }{2}`$, the dotted lines represent the theoretical optimal fidelities $`\frac{5}{6}`$ $`\left(0.833\right)`$, while the crosses $`(\times )`$ denote the fidelities of the output qubit $`1`$ and the triangles $`()`$ denote the fidelities of qubit $`2`$. As is shown, all the fidelities are close to the theoretical value $`\frac{5}{6}`$, and errors have been analyzed in section III.

## V Conclusion

We have realized a one-to-two qubits Bužek-Hillery cloning machine with linear optical devices. We showed that, the fidelities between the two output qubits and the original qubit are both $`\frac{5}{6}`$ for arbitrary input pure states.

Acknowledgment This work was supported by the National Natural Science Foundation of China.