text:
# Continuous melting of compact polymers

## Semiflexible loop model

In the Flory model a configuration of a polymer chain is described by a random walk on a square lattice (cubic, in three dimensions) which visits every site of the lattice exactly once. The energy associated with a configuration is $`Bϵ>0`$, where $`B`$ is the number of bends in the chain. At large temperatures ($`k_BTϵ`$) bends proliferate for entropic reasons, thus reducing the persistence length of the polymer. In the low-temperature limit ($`k_BTϵ`$) bends are very unlikely and the persistence length is large, of the order of the lattice size, as seen in fig. 1. The central question we address here is the nature of the transition between the disordered high-temperature phase and the ordered low-temperature phase.

The interplay between the entropy and the energy is encoded in the partition function of the Flory model

$$Z_F=\underset{\mathrm{\Gamma }_F}{}w^{B(\mathrm{\Gamma }_F)},$$
(1)
where $`w=\mathrm{exp}(ϵ/k_BT)`$ is the Boltzmann weight associated with a bend, and $`B(\mathrm{\Gamma }_F)`$ is the number of 90 bends in the compact chain configuration $`\mathrm{\Gamma }_F`$. For $`T0`$ ($`w0`$) the equilibrium state is one that minimizes $`B`$, while for $`T\mathrm{}`$ ($`w1`$) entropy dominates and there is a large, macroscopic number of bends in a typical chain configuration; see fig. 1. In order to investigate the transition between the high and low temperature phase we construct a mapping of the polymer problem to a lattice model of fluctuating loops, following Nienhuis nienrev .

We define the semiflexible loop model (SFL), whose configurations $`\mathrm{\Gamma }_L`$ are all possible ways of covering the square lattice with non-intersecting loops. The partition function of the SFL model reads

$$Z_L=\underset{\mathrm{\Gamma }_L}{}n^{N(\mathrm{\Gamma }_L)}w^{B(\mathrm{\Gamma }_L)},$$
(2)
where each of the $`N(\mathrm{\Gamma }_L)`$ loops in the configuration $`\mathrm{\Gamma }_L`$ carries a weight $`n>0`$, and bends are weighted by $`w`$ as in the Flory model. Note that the $`n0`$ limit of $`Z_L/n`$ yields the Flory model partition function for a single closed compact polymer. Instead of explicitly computing $`Z_L`$ the approach we take is to construct a continuum limit of the SFL model in the form of a field theory. Computations of correlation functions in this field theory lead to the phase diagram and exact scaling exponents for any $`n`$, including the polymer ($`n0`$) limit.

## Height model

Since loops in $`\mathrm{\Gamma }_L`$ do not intersect they can be interpreted as contour lines of the height of an interface. In the SFL model, the height is a three component scalar field. To see this, we first note that bonds not covered by loops form so-called ghost loops, if periodic boundary conditions are imposed. Then we assign to each real and ghost loop an orientation (clockwise or anticlockwise) independently and randomly. By doing so, each bond of the square lattice is found in one of four states, labeled $`𝐀,𝐁,𝐂,`$ or $`𝐃`$, depending on whether it is covered by a real or a ghost loop oriented in one or the other direction. We parametrize $`n=2\mathrm{cos}(\pi e)`$, with $`0e<1/2`$, and assign a weight $`\mathrm{exp}(\pm i\pi e)`$ (resp. $`\mathrm{exp}(\pm i\pi /3)`$) to each orientation of a real (resp. ghost) loop. Tracing over the two loop orientations, we recover the desired weights: $`n`$ for the real loops and $`1`$ for the ghost loops.

The heights for the SFL model are defined at the centers of the elementary square plaquettes of the lattice. Height differences between neighboring plaquettes are given by vectors $`𝐀,𝐁,𝐂,`$ or $`𝐃`$. Since each vertex is shared by four bonds, all necessarily in different states, a single relation, $`𝐀+𝐁+𝐂+𝐃=0`$, results from requiring that the total height increment when traversing a closed lattice path be zero. Thus, we conclude that out of the four vectors three are linearly independent and that the heights are three dimensional. For the remainder of the paper we adopt the following normalization: $`𝐀=(1,+1,+1)`$, $`𝐁=(+1,+1,1)`$, $`𝐂=(1,1,1)`$, and $`𝐃=(+1,1,+1)`$. The partition function of the SFL model can now be rewritten as a sum over height configurations, since a height configuration uniquely specifies an oriented loop configuration.

The continuum limit of the SFL model is obtained by coarse-graining the height field over domains much smaller than the system size and much bigger than a single plaquette. This yields a continuous height function $`𝐡(𝐱)=(h^1(𝐱),h^2(𝐱),h^3(𝐱))`$ defined over the plane $`𝐱`$. For $`n2`$ the height field is assumed rough JK\_npb  and the partition function can be written as a sum (path integral) over all height functions with weights given by $`\mathrm{exp}(S)`$, where

$$S=\frac{1}{2}\mathrm{d}^2𝐱K_{ij}(\mathbf{}h^i\mathbf{}h^j)+\frac{\mathrm{i}}{4\pi }\mathrm{d}^2𝐱(𝐞_0𝐡)\rho (𝐱)+\gamma \mathrm{d}^2𝐱\mathrm{cos}\left(2\pi h^2\right).$$
(3)
This particular form of the coarse-grained free energy $`S`$ is dictated by the symmetries of the loop model which involve translations in the space of heights, and permutations of the height-difference vectors JK\_npb . The value of the background charge

$$𝐞_0=\frac{\pi }{2}(1/3+e,0,1/3e)$$
(4)
follows from imposing the correct weighting of loops that close at infinity ($`\rho (𝐱)`$ is the curvature of the $`𝐱`$-surface), while the non-zero elastic constants

$$K_{11}=\frac{\pi }{8}(5/3e),K_{13}=\frac{\pi }{8}(e1/3),K_{22}=K_{22}(e,w)$$
(5)
are partially fixed by the loop ansatz which states that loop weights do not flow under renormalization jk\_prl .

In the previously studied infinite-temperature limit ($`w1`$), the elastic constant $`K_{22}=\pi (1e)/(53e)`$ is also fixed by the loop ansatz, and the cosine term is absent from $`S`$ JK\_npb . For $`w<1`$ the cosine potential accounts for the energy penalty associated with bends. Its scaling dimension in the $`w=1`$ theory, $`x_B=\pi /K_{22}=(53e)/(1e)`$, is greater than two for $`0e<1/2`$ ($`0<n2`$), and is therefore an irrelevant perturbation. It has the effect of renormalizing the elastic constant $`K_{22}`$ which ultimately drives the transition to the crystalline phase at $`w=w_c`$.

## Phase diagram and critical exponents

The field theory of the SFL model implies the following phase diagram for any $`n2`$, including the polymer limit $`n0`$. For $`w=1`$ the interface described by the height field is rough and the loop model is critical, i.e., it is characterized by a power-law distribution of loop sizes. As $`w`$ is reduced, the bending stiffness of the loops increases, and the height remains rough until a critical value $`w=w_\mathrm{c}`$. At $`w_\mathrm{c}`$ there is a roughening transition weeks  to a smooth interface; see fig. 1. In the loop model this corresponds to a phase with a vanishing density of bends. The roughening transition is continuous, in the Kosterlitz-Thouless universality class. At the transition the cosine term in $`S`$ becomes relevant in the renormalization group sense weeks . In other words, the scaling dimension of the operator $`\mathrm{cos}(2\pi h^2)`$ at $`w=w_\mathrm{c}`$ is 2. Computing this dimension using the field theory defined by $`S`$ then leads to the prediction $`K_{22}(e,w_\mathrm{c})=\pi /2`$. One important test of this result comes from the $`n=1`$ ($`e=1/3`$) SFL model, which maps to the exactly solvable six-vertex model. In this case $`w_\mathrm{c}=1/2`$ and $`K_{22}=\mathrm{arcsin}(1/2w)`$ have been computed previously baxter .

In the polymer limit the transition described by the field theory takes us from a critical geometry of a disordered compact polymer, at low bending stiffness, to an ordered state, at high bending stiffness; see fig. 1. The critical phase is characterized by scaling exponents $`\nu `$ and $`\theta `$. Namely, the probability distribution $`p(r,l)`$, for the end-to-end distance $`r`$ of a polymer of length $`l`$, has a scaling form $`p(r,l)=r^\theta f(r/l^\nu )`$ degennes . Using the field theory defined by eq. (3) these exponents can be computed exactly JK\_npb . We find

$$\nu =\frac{1}{2},\mathrm{and}\theta =\frac{3\pi 16K_{22}}{8\pi };$$
(6)
note that $`\theta `$ varies continuously with $`w`$ via the elastic constant $`K_{22}`$. The calculation of $`\nu `$ gives another important check on the field theory since $`\nu =1/2`$ follows directly from the fact that compact polymers completely fill the plane. Finally, exact values for $`K_{22}`$ in the $`w=1`$ and $`w=w_\mathrm{c}`$ case, lead to predictions $`\theta (w=1)=5/56`$ and $`\theta (w=w_\mathrm{c})=5/8`$. These values of $`\theta `$ have a simple physical interpretation. Namely, in the absence of bending stiffness ($`w=1`$) the ends of the chain repel (positive $`\theta `$), as is the case for open chain conformations of polymers in a good solvent. Somewhat surprisingly, this entropic repulsion decreases and turns to attraction (negative $`\theta `$) as the chain stiffness is increased.

## Numerical transfer matrix results

The main result of this paper, the prediction of continuous melting of compact polymers in two-dimensions, rests on the validity of the proposed field theory of the Flory model, eq. (3). To further test the field theory we have made use of a numerical transfer matrix approach JK\_npb  to compute the elastic constant $`K_{22}`$ and the exponent $`\theta `$. Results are shown in table 1. We see that eq. (6) is satisfied to a very good approximation and that the exact values for the two quantities at $`w=1`$ and $`w=w_\mathrm{c}`$ are confirmed by the numerical computations.

## Experiments

DNA absorbed on a lipid bilayer provides a laboratory for testing theories of two-dimensional polymers. Previous experiments DNA  on this system have measured the theoretically predicted scaling exponents derived from the two-dimensional self-avoiding walk model. These measurements were done at low DNA concentrations. As the amount of DNA on the surface is increased it was observed that single chain conformations become more compact. To study the transition discussed here compaction to areas of linear dimensions approaching the DNA persistence length would be needed. It would then be interesting to monitor the two ends of the DNA using fluorescent labels, as the amount of compaction is increased. This would provide direct information about $`p(r,l)`$ and the exponent $`\theta `$ computed above.

## Acknowledgements

We are grateful to K.A. Dill for introducing us to the Flory model. JK further thanks the KITP in Santa Barbara, where this work was initiated, for hospitality. JK is supported by the NSF under grant number DMR-9984471 and is a Cottrell Scholar of Research Corporation.


format: text

text:
# Squeezed-light source for the superresolving microscopy

quant-ph/0401025

## Abstract

We propose a source of multimode squeezed light that can be used for the superresolving microscopy beyond the standard quantum limit. This source is an optical parametric amplifier with a properly chosen diaphragm on its output and a Fourier lens. We demonstrate that such an arrangement produces squeezed prolate spheroidal waves which are the eigen modes of the optical imaging scheme used in microscopy. The degree of squeezing and the number of spatial modes in illuminating light, necessary for the effective object field reconstruction, are evaluated.

Spatial behavior of nonclassical light is presently attracting an increasing interest both in theory and experiment Kolobov99 ; EPJD . Quantum effects in optical imaging and other transverse spatial phenomena are studied within a framework of a European project ”Quantum Imaging” QUANTIM . One of the problems recently addressed in this context is about the ultimate quantum limit of the optical resolution Kolobov00 . A classical criterion of resolution was formulated at the end of the last century by Abbe and Rayleigh. It states that optical resolution is limited by diffraction on the system pupil. However, though it might be difficult to accept, the diffraction limit is not a fundamental limit like, for example the Heisenberg uncertainty relation. Recent publications have theoretically demonstrated that in principle this limit can be beaten for both writing and reading of optical information.

Sub-diffraction-limited optical recording of information, coined as ”quantum lithography”, uses a specially designed nonclassical light, $`N`$-photon entangled state, together with a $`N`$-photon absorber which allows to reach theoretically a resolution of $`\lambda /N`$ Boto00 ; Kok02 . To go beyond the diffraction limit in read-out of optical information like in imaging, is in principle possible having some a priori information about the object using the so-called ”superresolution” techniques. As shown recently in Ref. Kolobov00  the ultimate limit in such superresolution is set not by diffraction but by the quantum fluctuations of light illuminating the object of finite size and by the vacuum fluctuations outside the object. These quantum fluctuations of the illuminating light together with the vacuum fluctuations outside the object set up the standard quantum limit of superresolution which can be much smaller than the diffraction limit. Moreover, one can go beyond this standard quantum limit using a specially designed nonclassical light for illumination of the object. In recent experiment Treps02  it was demonstrated that using multimode squeezed light allows to measure a displacement of 2.9 Å of a laser beam with wavelength $`\lambda =1064`$ nm.

In Ref. Kolobov00  Kolobov and Fabre have suggested a scheme that allows to improve the superresolution beyond the standard quantum limit in reconstruction of an optical object using a multimode squeezed light. This scheme was formulated in terms of the so-called prolate spheroidal waves which are the eigenfunctions of the optical imaging scheme. To achieve the superresolution with multimode squeezed light one has to prepare these prolate waves in squeezed state. The question remained on how to produce such squeezed prolate waves.

In this letter we provide the answer to this question. Precisely, we demonstrate that an optical parametric amplifier (OPA) with a properly chosen diaphragm on its output and a Fourier lens, produces squeezed prolate spheroidal waves used in the superresolving microscopy. We investigate the quantum statistics of the squeezed prolate spheroidal waves in our scheme in dependence of the physical parameters of the OPA and of the optical configuration. We formulate simple estimates on the number of the ”object elements” to be reconstructed in connection with the number of degrees of freedom in the nonclassical illuminating light.

The scheme of optical imaging with multimode squeezed light is shown in Fig. 1. For simplicity we consider a one-dimensional case.

The part of the scheme to the right of the object plane performs diffraction-limited imaging of an object of finite size $`X`$ located in the object plane Kolobov00 . The first lens $`L_1`$ performs the Fourier transform of the object field into the pupil plane with a pupil of size $`d`$. After the second Fourier transform by the lens $`L_2`$ an image is created in the image plane. This image is a diffraction-limited copy of the object due to the finite size of the pupil.

The part to the left of the object plane is an illumination scheme. It consists of a traveling-wave OPA placed in the source plane and a Fourier lens $`L`$. It is well-known from the literature that a traveling-wave OPA with plane-wave pump and nonlinear crystal with large transverse area creates multimode squeezed vacuum on its output Kolobov99 . A new feature of our scheme is a diaphragm of size $`d_s`$ on the output of the OPA which serves for selection of the transverse modes in squeezed state. As we demonstrate below, when the size of this diaphragm matches the size of the pupil, $`d_sd`$, this setup squeezes exactly the prolate spheroidal waves which are the eigen modes of the imaging scheme. This result can be easily understood qualitatively. Indeed, when all three lenses in the scheme have the same focal distance $`f`$, as in Fig. 1, the lenses $`L`$ and $`L_1`$ create a geometrical image of the diaphragm $`d_s`$ in the pupil plane. Therefore, it is intuitively clear that one has to match the diaphragm size $`d_s`$ and the pupil size $`d`$ to select the modes of the source which will pass through the imaging scheme.

Let us introduce the dimensionless coordinates in the object and the image plane as $`s=2x/X`$, and the dimensionless coordinates in the source and the pupil plane as $`\xi =2y/d`$ (see Fig. 1). The dimensionless photon annihilation operators in the source, object, pupil, and image planes are denoted respectively as $`\widehat{c}(\xi ),\widehat{a}(s),\widehat{f}(\xi )`$, and $`\widehat{e}(s)`$. These operators obey the standard commutation relations,

$$[\widehat{c}(\xi ),\widehat{c}^{}(\xi ^{})]=\delta (\xi \xi ^{}),[\widehat{a}(s),\widehat{a}^{}(s^{})]=\delta (ss^{}),$$
(1)
and similar for $`\widehat{f}(\xi )`$, and $`\widehat{e}(s)`$. Naturally, the annihilation and creation operators at different planes commute with each other. These operators are normalized so that $`\widehat{c}^{}(\xi )\widehat{c}(\xi )`$, for example, gives the mean photon number per unit dimensionless length in the source plane. Similar normalization is valid for the other planes. The Fourier transform $`\left(T\widehat{c}\right)(s)`$ physically performed by the lens $`L`$ in Fig. 1 reads as follows,

$$\widehat{a}(s)=\left(T\widehat{c}\right)(s)=\sqrt{\frac{c}{2\pi }}_{\mathrm{}}^{\mathrm{}}𝑑\xi e^{ics\xi }\widehat{c}(\xi ),$$
(2)
where $`c={\displaystyle \frac{\pi }{2}}{\displaystyle \frac{dX}{\lambda f}}`$ is the space-bandwidth product of the imaging system.

The quantum theory of the imaging scheme in Fig. 1 was formulated in Ref. Kolobov00  in terms of the prolate spheroidal functions $`\psi _k(s)`$ Slepian61 ; Frieden71 (for some examples see Ref. Kolobov01 ). These are the eigen functions of the imaging operator of the scheme, orthonormal on the interval $`\mathrm{}<s<\mathrm{}`$. The photon annihilation operator $`\widehat{a}(s)`$ in the object plane can be written as decomposition over $`\psi _k(s)`$,

$$\widehat{a}(s)=\underset{k=0}{\overset{\mathrm{}}{}}\widehat{A}_k\psi _k(s)+\widehat{N}(s),$$
(3)
where the operator-valued coefficients $`\widehat{A}_k`$ are evaluated as

$$\widehat{A}_k=_{\mathrm{}}^{\mathrm{}}𝑑s\widehat{a}(s)\psi _k(s).$$
(4)
The operators $`\widehat{A}_k`$ and $`\widehat{A}_k^{}`$ satisfy the standard commutation relations of the photon annihilation and creation operators for discrete modes,

$$[\widehat{A}_k,\widehat{A}_k^{}^{}]=\delta _{k,k^{}},[\widehat{A}_k,\widehat{A}_k^{}]=0.$$
(5)
The Fourier transform $`\stackrel{~}{\psi }_k(\xi )`$ of the prolate spheroidal functions $`\psi _k(s)`$, performed by the lens $`L_1`$, is zero outside the interval $`|\xi |1`$. Consequently, the set of functions $`\{\psi _k(s)\}`$ is not complete in the Hilbert space $`L^2(\mathrm{},\mathrm{})`$, and to satisfy the commutation relations (1) one has to add an additional term $`\widehat{N}(s)`$. This term has zero Fourier spectrum in the interval $`|\xi |1`$ and does not contribute to the coefficients $`\widehat{A}_k`$,

$$_{\mathrm{}}^{\mathrm{}}𝑑s\psi _k(s)\widehat{N}(s)=0.$$
(6)
The field $`\widehat{N}(s)`$ can be decomposed over the complementary set of prolate functions $`\{\theta _k(s)\}`$, orthogonal to $`\{\psi _k(s)\}`$.

Physically speaking, the first term in Eq. (3) is the object field component that propagates in our scheme through the pupil to the image plane. The second object field component $`\widehat{N}(s)`$ is absorbed outside the pupil and is not observed. Therefore, in what follows we shall omit $`\widehat{N}(s)`$ in the object field (3).

To obtain the canonical transformation of the photon annihilation and creation operators in the imaging scheme on Fig. 1, one has to split the coordinates $`s`$ and $`\xi `$ into two regions, the ”core”, $`|s|1`$ and $`|\xi |1`$, corresponding to the central regions of the object (image) and the source (pupil), and the ”wings”, $`|s|>1`$ and $`|\xi |>1`$, outside these areas. The orthonormal bases in these regions of the object (image) plane are given by

$$\phi _k(s)=\{\begin{array}{cc}\frac{1}{\sqrt{\lambda _k}}\psi _k(s)& |s|1,\hfill \\ 0& |s|>1,\hfill \end{array}\chi _k(s)=\{\begin{array}{cc}0& |s|1,\hfill \\ \frac{1}{\sqrt{1\lambda _k}}\psi _k(s)& |s|>1,\hfill \end{array}$$
(7)
where $`\lambda _k`$ are the eigenvalues of the corresponding prolate spheroidal functions $`\psi _k(s)`$, depending on the space-bandwidth product $`c`$. We should note that the functions $`\phi _k(s)`$ are complete in the Hilbert space $`L^2(1,1)`$.

It follows from (7) that the functions $`\psi _k(s)`$ and $`\theta _k(s)`$ can be written in the form

$$\psi _k(s)=\sqrt{\lambda _k}\phi _k(s)+\sqrt{1\lambda _k}\chi _k(s),\theta _k(s)=\sqrt{1\lambda _k}\phi _k(s)\sqrt{\lambda _k}\chi _k(s).$$
(8)
Similar relations take place in the source (pupil) plane.

In terms of two sets $`\{\phi _k(s)\}`$ and $`\{\chi _k(s)\}`$ we can write the annihilation operators in the source and the object planes as

$$\widehat{c}(\xi )=\underset{k=0}{\overset{\mathrm{}}{}}\widehat{c}_k\phi _k(\xi )+\underset{k=0}{\overset{\mathrm{}}{}}\widehat{d}_k\chi _k(\xi ),$$
(9)
and

$$\widehat{a}(s)=\underset{k=0}{\overset{\mathrm{}}{}}\widehat{a}_k\phi _k(s)+\underset{k=0}{\overset{\mathrm{}}{}}\widehat{b}_k\chi _k(s).$$
(10)
Here $`\widehat{c}_k`$ and $`\widehat{a}_k`$ are the annihilation operators of the prolate modes $`\phi _k`$ in the core region of the source and the object planes, while $`\widehat{d}_k`$ and $`\widehat{b}_k`$ are the annihilation operators of the prolate modes $`\chi _k`$ in the wings regions. The operators $`\widehat{c}_k`$ and $`\widehat{d}_k`$ are expressed through the field operator $`\widehat{c}(\xi )`$ by

$$\widehat{c}_k=_{\mathrm{}}^{\mathrm{}}𝑑\xi \widehat{c}(\xi )\phi _k(\xi ),\widehat{d}_k=_{\mathrm{}}^{\mathrm{}}𝑑\xi \widehat{c}(\xi )\chi _k(\xi ).$$
(11)
Similar relations hold for $`\widehat{a}_k`$, $`\widehat{b}_k`$ and $`\widehat{a}(s)`$.

In our analysis we shall use the following property of prolate spheroidal functions Frieden71 ,

$$_1^1𝑑\xi \psi _k(\xi )e^{iq\xi }=(i)^k\sqrt{\frac{2\pi \lambda _k}{c}}\psi _k(q/c).$$
(12)
Using (7), (8), the field transform (2) between the source and the object plane and this equation, one can find the following propagation relations for the core and wings of the light wave emitted by the source:

$$\left(T\phi _k\right)(s)=(i)^k\left[\sqrt{\lambda _k}\phi _k(s)+\sqrt{1\lambda _k}\chi _k(s)\right]=(i)^k\psi _k(s),$$

$$\left(T\chi _k\right)(s)=(i)^k\left[\sqrt{1\lambda _k}\phi _k(s)\sqrt{\lambda _k}\chi _k(s)\right]=(i)^k\theta _k(s).$$
(13)
Taking into account (9) and (10) we obtain,

$$\widehat{a}(s)=\left(T\widehat{c}\right)(s)=\underset{k=0}{\overset{\mathrm{}}{}}(i)^k\left[\widehat{c}_k\psi _k(s)+\widehat{d}_k\theta _k(s)\right],$$
(14)
and

$$\widehat{A}_k=(i)^k\widehat{c}_k.$$
(15)
As expected, there is no contribution into $`\widehat{A}_k`$ from the second sum in Eq. (9) containing operators $`\widehat{d}_k`$ and describing the illumination coming from the wings area of the source. If the source diaphragm is larger or matches the size of the pupil, $`d_sd`$, it has no effect on the operator amplitudes $`\widehat{c}_k`$ and $`\widehat{A}_k`$, see (11), and we shall neglect the diaphragm in the calculation of $`\widehat{A}_k`$.

To obtain explicitly the coefficients $`\widehat{A}_k`$ in case of illumination of the scheme by a traveling-wave OPA we shall use simplified description of an OPA with a plane-wave undepleted pump (see Ref. Kolobov99 ). In this approximation one can find analytically the spatial Fourier amplitudes $`\stackrel{~}{c}(q)`$,

$$\stackrel{~}{c}(q)=_{\mathrm{}}^{\mathrm{}}𝑑\xi e^{iq\xi }\widehat{c}(\xi ),$$
(16)
of the field at the output of the crystal as a linear transformation of the corresponding input Fourier amplitudes $`\stackrel{~}{c}_{in}(q)`$ and $`\stackrel{~}{c}_{in}^{}(q)`$,

$$\stackrel{~}{c}(q)=U(q)\stackrel{~}{c}_{in}(q)+V(q)\stackrel{~}{c}_{in}^{}(q).$$
(17)
Here the operators $`\stackrel{~}{c}_{in}(q)`$ at the input of the OPA are in the vacuum state. The complex coefficients $`U(q)`$ and $`V(q)`$ depend on the nonlinear susceptibility of the crystal, its length, and the matching conditions in the OPA. These coefficients have the property,

$$|U(q)|^2|V(q)|^2=1,$$
(18)
that guarantees the preservation of the commutation relations (1). Using the standard parameters of multimode squeezing (see Ref. Kolobov99 ) we obtain:

$$U(q)\pm V^{}(q)=e^{i\phi (q)}\left[e^{\pm r(q)}\mathrm{cos}\theta (q)+ie^{r(q)}\mathrm{sin}\theta (q)\right].$$
(19)
The phase of the amplified (stretched) quadrature amplitude of the OPA output field $`\stackrel{~}{c}(q)`$ is $`\theta (q)`$. The phase factor $`\mathrm{exp}(i\phi (q))`$ specifies the phase of the quadrature amplitudes of the input field $`\stackrel{~}{c}_{in}(q)`$ which are squeezed or stretched. For the vacuum input field of the OPA the last phase is irrelevant.

The operator amplitudes $`\widehat{A}_k`$ are found explicitly with the use of Eqs. (15) and (11). The quantities $`\widehat{c}_{in}(\xi )`$ and $`\phi _k(\xi )`$ are expressed through their Fourier transforms (17) and (13). After some calculation we obtain,

$$\widehat{A}_k=\frac{1}{\sqrt{2\pi c}}_{\mathrm{}}^{\mathrm{}}𝑑q\psi _k(q/c)\left[U(q)\stackrel{~}{c}_{in}(q)+V(q)\stackrel{~}{c}_{in}^{}(q)\right].$$
(20)
We shall introduce the real quadrature components of the field amplitudes $`\widehat{A}_k`$ in the object plane as

$$\widehat{A}_k=\widehat{A}_{1k}+i\widehat{A}_{2k}.$$
(21)
For the variances of these quadrature components we obtain

$$\left(\mathrm{\Delta }\widehat{A}_{1k}\right)^2=\frac{1}{4c}_{\mathrm{}}^{\mathrm{}}𝑑q\psi _k^2(q/c)\left[e^{\pm 2r(q)}\mathrm{cos}^2\theta (q)+e^{2r(q)}\mathrm{sin}^2\theta (q)\right],$$
(22)
$$\left(\mathrm{\Delta }\widehat{A}_{2k}\right)^2=\frac{1}{4c}_{\mathrm{}}^{\mathrm{}}𝑑q\psi _k^2(q/c)\left[e^{2r(q)}\mathrm{cos}^2\theta (q)+e^{\pm 2r(q)}\mathrm{sin}^2\theta (q)\right].$$
(23)
Here the upper and the lower sign correspond respectively to the even and the odd prolate spheroidal functions $`\psi _k`$.

It follows from this result that the prolate spheroidal waves in the object illumination can be prepared in squeezed state. By proper choice of the squeezing phase $`\theta (q)`$ at low spatial frequencies $`q`$ one can minimize quantum fluctuations in one of the quadrature amplitudes $`A_{\sigma k}`$, $`\sigma =1,2`$ (namely, in the one detected in the image plane of our scheme). Taking the degree and the phase of squeezing as constant, $`r(q)=r`$, $`\theta (q)=\theta `$, we can estimate the variance of the squeezed quadrature amplitude:

$$\left(\mathrm{\Delta }A_{\sigma k}\right)^2e^{2r}/4.$$
(24)
In reality the spatial-frequency band of multimode squeezing is limited by the phase-matching condition in the OPA. One has to take into account the spatial-frequency dispersion of squeezing, that is, the frequency dependence of the squeezing phase $`\theta (q)`$ due to diffraction in free space and inside the OPA. Both phenomena deteriorate squeezing of prolate spheroidal waves.

As shown in Kolobov89 ; Kolobov99 ; Sokolov99 , the effect of diffraction on squeezing can be almost perfectly compensated by means of adjustment of the lens array. If this is done, it follows from (22), (23), that the degree of squeezing of prolate spheroidal waves depends on the overlap in the object plane of two areas: (i) the area illuminated by the squeezed plane waves $`\stackrel{~}{c}(q)`$, which are focused by lens $`L`$ to the points $`s=q/c`$, and (ii) the area of the energy localization $`\psi _k^2(s)`$ of the prolate spheroidal waves.

In analogy to other phenomena with spatially-multimode squeezed light, in the even and odd components of the object field the different quadrature amplitudes are squeezed; for, say, $`\theta =0`$, the squeezed quadratures are $`A_{1k}`$ and $`A_{2k}`$ for the odd and the even prolate waves respectively. In measurement this implies, e.g., the spatially resolved homodyne detection with numerical evaluation of the odd or even squeezed amplitudes (depending on the local oscillator phase).

Finally, we can formulate the conditions on multimode squeezing in the object illumination in terms of the number $`N`$ of independent degrees of freedom in the light field, propagating through the diaphragm of the OPA (we assume here the optimum size $`d_s=d`$ of the diaphragm). The properties of the OPA emission can be characterized by the coherence length $`l_c`$ of the output field $`\widehat{c}(\xi )`$. The spatial-frequency range $`q_c`$ of effective squeezing is related to the coherence length by the estimate

$$|q_c|\frac{\pi d}{2l_c}.$$
(25)
The minimum requirement on the OPA is that the effectively squeezed waves $`\stackrel{~}{c}(q_c)`$ illuminate, after passing the lens $`L`$, the object region $`|s|1`$. That is, the waves

$$|q_c|/c1,$$
(26)
should be squeezed. This gives the following estimate for the coherence length,

$$l_c\frac{\pi d}{2c},$$
(27)
and for the number $`N`$ of independent degrees of freedom in the illuminating light, emitted from the region $`d`$:

$$Nd/l_c\frac{dX}{\lambda f}=S.$$
(28)
Here $`S`$ is the Shannon number of our optical scheme.

As seen from (8), the wave profiles $`\psi _k(s)`$ can be expanded in terms of two bases, $`\{\phi _k(s)\}`$ and $`\{\chi _k(s)\}`$, representing the core and the wings of the illuminating field in the object plane, where $`\lambda _k1`$ are the eigenvalues of the imaging transformation. As known from the theory of prolate spheroidal functions Slepian61 ; Frieden71 , these eigenvalues are close to 1 only for $`kS`$, where $`S`$ is the Shannon number. For higher values of the index $`k`$ the field energy in prolate spheroidal waves $`\psi _k(s)`$ is concentrated in the wings, that is, outside the object area $`|s|1`$. Hence, the condition $`NS`$ provides an effective squeezing only of the prolate spheroidal waves with $`kS`$. In order to minimize quantum noise of the higher prolate spheroidal waves (with $`k>S`$), it is necessary to use OPA with large number of effectively squeezed spatial modes of radiation,

$$NS,$$
(29)
and to illuminate by non-classical light a spot in the object plane with the size much larger than the object itself.

When the condition $`NS`$ is met, one can illuminate the object (within the area $`|s|1`$) by a bright classical wave with a properly chosen phase. Since for $`\lambda _k1`$ our imaging scheme is analogous to a highly non-symmetrical beamsplitter Kolobov00 , the squeezed illumination concentrated mainly in the wings area of the object plane allows for the low-noise measurement of the relevant amplitudes of the image field.

This work was supported by the Network QUANTIM (IST-2000-26019) of the European Union and by the INTAS under Project No. 2001-2097.


text:
# An explanation for metallicity effects on X-ray Binary properties

## 1 Introduction

Globular clusters (GCs) are important laboratories for studying stellar populations, both main sequence and exotic. X-ray binaries are about 100 times overabundant in GCs compared to in the field compared to the field (see e.g. the compilation of X-ray binary properties in Liu, van Paradijs & van den Heuvel 2001) due to dynamical effects such as tidal capture and/or three and four body interactions (e.g. Clark 1975; Fabian, Pringle & Rees 1975; Hills 1976). Furthermore, since detailed population studies of X-ray binaries require more sources than the $`100`$ in the Local Group, the extragalactic case must be studied. Given that the ages and metallicities of individual extragalactic field stars are nearly impossible to determine, but that these properties can be inferred from integrated light of GCs, the GC X-ray binaries in other galaxies represent an ideal place to study how age and metallicity affect X-ray binary populations.

Clearly, the dynamical properties of a GC are most important for predicting whether it will have an X-ray binary (e.g. Pooley et al. 2003; Heinke et al. 2003). In addition, there also exists a strong residual correlation between the number densities of X-ray binaries and the metallicities of the GCs. This was first suggested to be the case in the Local Group (Grindlay 1993; Bellazzini et al. 1995), and was shown more conclusively and dramatically in NGC 4472 (Kundu, Maccarone & Zepf 2002 - KMZ02). The possibility that this represents an age effect rather than a metallicity effect has been tested and found not to be the case using data on NGC 3115 and NGC 4365 (Kundu et al. 2003 - K03). In NGC 3115, where the age spread is small, but the metallicity spread is large, the metallicity effect was found to be at least as strong as that in other early type galaxies. Conversely, in NGC 4365, where the age spread is large, the fraction of globular clusters with LMXBs did not vary much with age and was similar to that of older GCs with similar metallicities in other galaxies. It thus seems most likely that the LMXB population enhancement effects first associated with the metallicity are, in fact, related directly to the metallicity and not to some other correlated property such as age, half-light radius, or distance from the center of the galaxy. It has also been noted that the metal poor Large Magellanic Cloud has a lower ratio of LMXBs to HMXBs than the more metal rich Milky Way (Cowley 1994; Iben, Tutukov & Fedorova 1997 - ITF97), which may be a combination of the metallicity effects and the difference in their star formation histories.

An additional metallicity effect can be found in the differences of the soft X-ray spectra of these sources. It has been found that the spectra of the blue (i.e. metal poor) GCs are harder than those of the red (i.e. metal rich) GCs in the Milky Way and M31 (Irwin & Bregman 1998 - IB) and in NGC 4472 (Maccarone, Kundu & Zepf 2003 - MKZ03). On the other hand, the quiescent LMXBs seem to show no evidence of a metallicity effect on the spectrum (e.g. Heinke et al. 2003), and the effect becomes much weaker or disappears at higher X-ray energies (Trinchieri et al. 1999; Di Stefano et al. 2002; Sidoli et al. 2001; MKZ03).

In this paper, we will show that the current theoretical explanations for the overabundance effect are unlikely to match the typical factor of 3 difference between the probabilities of metal poor and metal rich GCs’ having X-ray LMXBs. We will put forth a new scenario invoking the effects of irradiation induced winds (IIWs) can explain the population difference. We will also show that IIWs can explain the previously unexplained spectral difference effects as a result of the same physical process, with the same parameter values.

## 2 Past Theoretical Work

### 2.1 IMF variations?

Several attempts have been made to explain why metal rich GCs have more LMXBs, but as yet, none seems satisfactory. Grindlay (1993) suggested that the correlation might be due to a flatter initial mass function in higher metallicity GCs. The general consensus seems to be that the IMF is fairly universal (e.g. Kroupa 2002), and star formation theory suggests that any metallicity dependence is likely to be such that metal rich stars have typically lower masses since the Jeans mass will be lower for the more efficiently cooling metal rich gas (Larson 1998). Although the present day mass functions of metal rich Galactic GCs are flatter than those of metal poor Galactic GCs (McClure et al. 1986), this seems to be largely due to the fact that the metal rich systems are located more centrally within the Milky Way and hence experience more extreme dynamical stripping of their low mass stars, as seen observationally (Piotto & Zoccali 1999) and expected theoretically (Vesperini & Heggie 1997; Baumgardt & Makino 2003).

### 2.2 Effects of stellar radius

Later, Bellazzini et al. (1995) suggested that the larger stellar radii of the stars in metal rich globular clusters might contribute to their increased number densities of X-ray binaries. Larger stars should have a higher cross section for tidal captures and should overflow their Roche lobes at larger separations. The physics of tidal capture are a hotly contested issue and some workers have suggested that tidal captures have great difficulty in forming systems similar to the X-ray binaries we observe (see e.g. Rasio & Shapiro 1991; Rasio, Pfahl & Rappaport 2000 and references within), while others have suggested that the bulk of recycled pulsars may have been formed in tidal capture systems (e.g. Di Stefano & Rappaport 1992). The alternative to formation of X-ray binaries by tidal captures is formation through exchange interactions. These are most commonly three-body interactions where a neutron star encounters a binary system composed of two non-compact stars and replaces one of those two stars and forms a binary system with the other (see e.g. Clark 1975; Hills 1976), or four-body interactions, where the neutron star is already in a binary system when it interacts with another binary system (see e.g. Mikkola 1984; see also Fregeau et al. 2003 for a discussion of recent numerical work on 3-and 4-body interactions).

The easier Roche lobe overflow can be calculated using the standard assumption of a binary separation distribution which is logarithmic, i.e. $`P(a)`$$`1/a`$. Taking Kepler’s law and the period-mass relation \[equation (4.9) of Frank, King & Raine (1992 - FKR92)\] we find that for a neutron star accretor and a typical (i.e. a 0.6 $`M_{}`$) main sequence GC star donor, the Roche lobe overflow orbital separation will be about 3 stellar radii, and varies little with the donor star’s mass. The separation at which Roche lobe overflow will occur increases linearly with the stellar radius for a fixed stellar mass. Given the logarithmic distribution of separations and a metallicity dependence of radius such that $`R_{}`$$`Z^{1/8}`$ (estimated from the stellar model interpolation formulae of Tout et al. 1996), the contribution of the term related to the increase in the number of Roche lobe overflowing systems will be $`1+\frac{1}{8}\mathrm{log}(\frac{Z_r}{Z_p})/\mathrm{log}(\frac{r_{max}}{r_{min}})`$, which is about 1.3 for a metallicity ratio of 10. Numerical calculations suggest that the tidal capture rate goes as $`R_{}^{}{}_{}{}^{0.93},`$ (e.g. Lee & Ostriker 1986), so the number of neutron star binaries formed by tidal capture will go as $`Z^{0.12}`$ which gives another factor of about 1.3 for a ratio of metallicities of 10. An additional small increase in the tidal capture probability might come from the higher turnoff masses for the stars in metal rich GCs (see e.g. Chaboyer et al. 1996), but this factor should be no more than about 5-10%. Stellar radius effects thus seem unlikely to make more than a $``$60% difference (i.e. the product of the two 30% effects) in X-ray binary populations as a function of metallicity. In fact, since tidal captures are not likely to produce binary systems that are similar to the observed parameters, it seems more likely that the bulk of the LMXBs are formed in three and four body exchanges (see e.g. Rasio et al. 2000), and so the overabundance factor in red clusters due to larger stellar radii is likely to be closer to the factor of 1.3 than 1.6.

## 3 The irradiation induced wind model

The mass donors in X-ray binaries can absorb and reprocess luminosities comparable to their intrinsic luminosities. The mass donor may then be “puffed up” to larger radii (e.g. Tutukov & Yungelson 1980; Podsiadlowski 1991; Harpaz & Rappaport 1991) and that the extra kinetic energy added to the envelope may drive off an evaporative wind with a velocity of order the escape velocity from the stellar surface (Arons 1973; Ruderman et al. 1989; Tavani & London 1993; Pfahl, Rappaport & Podsiadlowski 2003). A fraction ($``$ 10%) of the gas lost by the mass donor will be accreted by the compact star in much the same way that more typical stellar winds are accreted in HMXBs. The lifetimes of LMXBs will also be accelerated by the extra mass loss in the IIWs (Tavani 1991a).

In addition to driving the evolution of the system, the irradiation induced wind will also leave a large amount of gas in the environment of the X-ray binary. Most high mass X-ray binaries show clear orbital modulations of the X-ray flux, especially at soft X-ray energies, and this is taken to be evidence of internal absorption by the material in the stellar wind. An IIW from a low mass star will likely have a velocity much closer to the orbital velocities in the binary system, so the orbital modulation will not necessarily be as strong as in the HMXBs. Still, the increased absorption will have an effect on the spectrum.

Let us summarize a few key past results on irradiation induced winds. It is generally found that these winds should be more important in low mass X-ray binaries than in high mass X-ray binaries, since in low mass X-ray binaries the amount of radiation absorbed by the mass donor from the accretion flow may far exceed the nuclear energy generation rate and may hence have a significant effect on the structure of the donor star (see e.g. Tutukov & Yungelson 1980; Podsiadlowski 1991; ITF97); IIWs may affect HMXBs as well (see e.g. Day & Stevens 1993), but have been much less well studied and are far less likely to dominate the overall mass transfer. Strong coronal winds are likely to result from IIWs and may be self-sustaining even if the mass donor does not fill its Roche lobe (see e.g. Basko & Sunyaev 1973; Arons 1973); in fact, ITF97 have found that self-sustaining winds can produce large luminosities even if the mass donor fills only $`80\%`$ of its Roche lobe, in agreement with past results that a mass donor need not fill its Roche lobe if it is sufficiently irradiated (Tavani, Ruderman & Shaham 1989). Some initial Roche lobe overflow is likely to be required in order to start accretion, but if the irradiation induced winds cause the mass loss to be substantially faster than what would be caused by orbital and stellar evolution, then the star may cease to fill its Roche lobe even as mass loss and accretion continue. Metal rich stars can dissipate much of their absorbed energy through line cooling, while metal poor stars dissipate this energy primarily through IIWs. Interpolating from the stellar cooling rate tables of Sutherland & Dopita (1993), the mass loss rate due to the IIWs should scale as $`Z^{0.35}`$; the system lifetime should scale as $`Z^{0.35}`$ if the mass loss is dominated by these winds and the lifetime is determined by the timescale for the mass donor to lose all its mass (ITF97).

As a caveat, we note that the treatment of ITF97, is based in part on the analytical irradiation treatment in Iben, Tutukov & Yungelson (1995); while the results do agree well with the numerical work of Tavani & London (1993) within the parameter space of their models, but ITF97 extrapolate outside this range. A more sophisticated numerical treatment may be in order for future work, but is clearly outside the scope of this paper.

For a wind emitted at the escape velocity from the mass donor, equation (4.35) of FKR92 shows that the fraction of the wind captured by the compact object will go as:

$$\frac{\dot{M}}{\dot{M}_w}=\frac{1}{4}(\frac{M_{CO}}{M_d})^2(\frac{R_d}{a})^2$$
(1)
where $`\dot{M}`$ is the mass accretion rate of the compact object $`\dot{M}_w`$ is the wind mass loss rate, $`M_{CO}`$ is compact object’s mass, $`M_d`$ and $`R_d`$ are the mass and radius of the donor star, and $`a`$ is the orbital separation. Given typical values of $`a=4\times 10^{11}`$ cm, $`R_d=7\times 10^{10}`$ cm, $`M_d=0.6M_{}`$ and $`M_{CO}=1.4M_{}`$, about 5% of the wind is accreted.

To simplify the calculation, we make one additional assumption, that the red GCs accrete from a combination of Roche lobe overflow and IIWs, while the bright LMXBs in blue GCs accrete from a self-sustaining wind. We will later revisit this assumption and show that it is not necessary in order to reproduce roughly the observations. We now define the notation for the following 5 equations - $`\dot{m}`$ indicates the mass accretion rate onto the compact object, while $`\dot{M}`$ is the mass loss rate. The subscript $`w`$ indicates wind mass loss, $`RL`$ indicates Roche lobe overflow effects, $`p`$ indicates systems with metal poor donors, and $`r`$ indicates systems with metal rich donors, $`N`$ indicates the number of systems and $`Z`$ indicates metallicity. Then

$$\dot{m}=ϵ\dot{M}_w+\dot{M}_{RL},$$
(2)
where $`ϵ`$ is the fraction of the mass lost in the wind that is accreted by the compact object and is given by equation (1). Then, starting from the assumption outlined above that the Roche lobe overflow component is neglible for the metal poor systems, we have:

$$\dot{m}_p=ϵ\dot{M}_{w,p},$$
(3)
while for metal rich stars,

$$\dot{m}_r=ϵ(\frac{Z_r}{Z_p})^{0.35}\dot{M}_{w,p}+\dot{M}_{RL},$$
(4)
where $`\dot{M}_{RL}`$ is the mass loss rate due to the Roche lobe overflowing component of the accretion flow, and is assumed to be much smaller than the irradiation driven mass loss rate for the case of a metal poor donor. Combining equations (3) & (4), we find that, for the same luminosity,

$$\frac{\dot{M}_r}{\dot{M}_p}=(\frac{Z_r}{Z_p})^{0.35}+ϵϵ(\frac{Z_r}{Z_p})^{0.35}.$$
(5)
The fraction of the mass loss in the metal rich systems coming from the IIW is $`(\frac{Z_r}{Z_p})^{0.35}`$.

The number of X-ray binaries should scale as the formation rate times the lifetime. The formation rate effects have been studied by Bellazzini et al. (1995), and given their lines of argument, we found in Section 2.2 that the stellar radius effects should produce a difference by a factor of about 1.3 (if exchange interactions dominate) to 1.6 (if tidal captures dominate).

The effects of IIWs are predominantly on the system lifetimes. For systems at a given luminosity, equation 5 shows the difference in mass loss rate, the inverse of which gives the ratio of source lifetimes. Thus we find that

$$\frac{N_r}{N_p}=\frac{1}{(\frac{Z_r}{Z_p})^{0.35}+ϵϵ(\frac{Z_r}{Z_p})^{0.35}}.$$
(6)
To compute the actual ratio of the number of red to blue GC X-ray sources, it is necessary to multiply the factor of 1.3 to 1.6 from the stellar radius effects by the value from equation 6, which should be about 2.1 for the typical parameters $`ϵ=0.05`$ and $`\frac{Z_r}{Z_p}`$=10. This gives a factor between about 2.6 and 3.4, although this factor is likely a slight overestimate, because even the X-ray sources in the most metal poor GCs are likely to have at least some Roche lobe overflow contribution to their X-ray luminosities.

It has been assumed that the stellar wind velocity is equal to the escape velocity from the surface of the star; this need not be the case. The extra wind energy for the metal poor stars may be dissipated as a higher velocity wind rather than as a more dense wind. The fraction of the mass lost that is accreted scales as $`v_{wind}^4`$ (FKR92), which alternatively scales as the inverse of the wind power squared, for a constant mass loss rate. The luminosity of the LMXBs in metal poor GCs would then be surpressed by a factor of about $`(\frac{Z_r}{Z_p})^{0.70}`$, while the lifetimes of the two classes of systems would be about the same. Because the luminosity function has a slope of about $`0.55`$ (KMZ02), the ratio of the number of metal rich and metal poor cluster X-ray systems would then be $`(\frac{Z_r}{Z_p})^{0.39}`$, which, for the canonical factor of 10 difference between the two modes gives a factor of about 2.5 difference in the expected number of observed systems. We do note that the there might be systematic variations in the slope of the luminosity functions as a function of metallicity, but absent measurements or a theoretical model, we assume they will be the same.

We wish now to estimate the contribution to the column density with which these systems will typically be observed due to the IIW. We find an average density of mass in a sphere around the mass donor with radius equal to the diameter of the orbit. For a path length of the orbital radius, the column density $`N_H`$ is then:

$$N_H=\frac{\dot{M}}{8v_wr_{orb}\mu },$$
(7)
where $`\mu `$ is the mean molecular weight of the gas in the wind, $`v_w`$ is the wind velocity, and $`r_{orb}`$ is the orbital separation. This value, $`6\times 10^{21}`$cm-2, should give a good approximation over orbital phase and inclination angle for the column density observed. Edge-on sources might be expected to have higher column densities, but the inclination angles for the GC LMXBs are not well constrained.

## 4 Observational evidence for the scenario

### 4.1 NGC 4472 - number of sources

NGC 4472 is the first galaxy where the metal rich mode was shown to have a higher fraction of GCs with LMXBs than the metal poor mode (KMZ02). It seems unlikely that the GCs in NGC 4472 span a wide range of ages; they are all likely to be within a factor of 1.5 to 2 in age (Beasley et al. 2000; Cohen, Blakeslee & Côté 2003). The metallicity effects on the number of globular cluster X-ray sources thus are not likely due to an age difference between the metal rich and metal poor GCs. The age measurements are rather sensitive to the stellar populations models for the Balmer lines and there could still be a rather substantial age difference between the two samples. More strict constraints are the age measurements of Puzia et al. (2002) which confirm that the correlation between metallicity and LMXB specific frequency in NGC 3115 is due to metallicity and not age (K03), but this system has fewer X-ray sources, so the ratio of the number of LMXBs in metal rich and poor clusters cannot be as well determined. Finally there is the case of NGC 4365, where the ages do span a rather wide range, but do not seem to be strongly correlated with the LMXB number density (K03).

The two modes in color for NGC 4472 peak at $`VI`$ of 0.98 and 1.23 (KMZ02), corresponding to values of \[Fe/H\] of -1.26 and -0.08, respectively, according to the scaling law of Kundu & Whitmore (1998). Defining the metal rich/metal poor mode boundary to be $`VI`$=1.10, we find that 23 of the 450 metal rich GCs and only 7 of the 370 metal poor GCs contain X-ray sources. The metal rich GCs are thus $`2.7\pm 1.2`$ times as likely to contain X-ray sources as the metal poor GCs.

### 4.2 NGC 4472 - source spectra

The spectra also show a difference as a function of metallicity. While the individual spectra cannot be easily measured because of the low count rates, we have found that the summed spectra in NGC 4472 are harder in the blue GCs than in the red ones (MKZ03). If we hold the neutral hydrogen column in both cases to the Galactic value of 1.6$`\times 10^{20}`$ cm-2 and fit a power law model to the data, we find a spectral index of 1.02$`\pm `$.27 for the blue GCs while we find a spectral index of 1.46$`\pm `$.10 for the red GCs (90% error contours). Allowing the column to float freely for the red GCs, we find $`N_H`$ to be 4.8$`\times 10^{20}`$ cm-2 and the power law index to be 1.57. Then, we fix the power law index for the blue GCs and find that the data is best fit with a column density of 1.1$`\times 10^{21}`$ cm-2.

We note that the solution to the spectral difference problem is not unique and is prone to numerous systematic uncertainties. The photoelectric absorption models we have used assume a solar composition for the absorbing medium, and that the medium is cold (i.e. completely un-ionized). The underlying spectrum for accreting neutron stars and black holes in the 0.5-8 keV range is unlikely to be a single power law. Finally, we have averaged over many systems with different values of $`N_H`$ and with different underlying spectra. Still, the rough information given, that the metal poor GCs have an intrinsic absorption of about $`10^{21}`$ cm-2 and that the column density is about 3 times as large for the metal poor GCs as it is for the metal rich GCs seems to be a reasonable inference to draw from the data. The theoretical model predicts a higher value for the blue clusters, but the linear averaging tends to overestimate the effects on the spectrum, so the fact that our crude calculation over-predicts the amount of absorption is to be expected. Furthermore, the fact that the gas is likely to be partly ionized and will make the fitted value of the $`N_H`$ less than the actual value, and also some of the gas mass will condense into a geometrically thin accretion disk, and hence will have an effect of absorbing X-rays only if the inclination angle is very low. That the values are on the same order of magnitude and that the blue clusters have about 3 times as much absorption in the fits is about as good an agreement as can be expected given the crude modelling and the considerable theoretical uncertainties in the models of IIWs.

### 4.3 Local Group Sources

Both the Milky Way and the Magellenic Clouds have been rather well studied in terms of their X-ray stellar populations, and they, show a metallicity difference on the same order as that between the metal rich and metal poor modes for GCs - about a factor of 10. As only a small fraction of the X-ray sources in any of these galaxies is in a GC, the star formation rate has a substantial effect on the relative number densities of X-ray binaries, so merely comparing number counts per unit stellar mass is not likely to prove frutiful. However, one can be fairly confident that the high mass X-ray binary population is not heavily affected by IIWs because the luminosities of high mass stars are much larger than the intercepted and absorbed luminosities. Therefore, the ratio of LMXBs to HMXBs might give a rough estimate of how important the IIWs are. The suggestion of ITF97 that the difference of this ratio might be indicating that irradiation induced winds are playing an important role is therefore additional evidence in favor of this scenario.

## 5 Potential observational tests

This model makes several testable predictions. The first is that there should be a monotonic dependence between the number of LMXBs per unit stellar mass and the metallicity; given enough statistics we should see a difference in the LMXB specific frequency as a function of the metallicity itself, and not just as a function of whether a cluster is in the metal rich or metal poor mode. Much new data has recently entered the Chandra archives, so it is now possible to test this prediction. There is already a correlation over a range of metallicities in the ROSAT spectral indices of GC X-ray sources which does not show a “critical metallicity” (IB), so our model seems to pass this test so far.

This scenario also predicts that neutron star LMXBs will be affected far more than other types of “dynamically interesting” sources. Blue stragglers and cataclysmic variables will not generate high enough X-ray luminosities to excite substantial IIWs. Black hole systems, because of the higher mass compact objects, will accrete the IIW gas much more efficiently.

IIWs have been suggested to explain why systems such as 4U 1820-30 show different period evolution than would be expected from conservative mass transfer driven by gravitational radition (Tavani 1991b). An alternative is that the system is being effected by the interactions with the GC potential (van der Klis et al. 1993). Gravitational wave observations from future missions such as LISA may help break this degeneracy.

Additionally, our scenario predicts that the metallicity effects should be essentially the same for field sources as they are for globular cluster sources. Given two galaxies with similar star formation histories, the more metal rich galaxy should have more field X-ray binaries. The natural way to test this hypothesis would be to look at the field X-ray binary populations of elliptical galaxies, as (1) they have very little recent star formation and (2) metallicity tends to scale with galaxy mass. A potential problem with this approach is that a fraction (and indeed, perhaps a large fraction) of the field X-ray binaries may have been created through stellar interactions in GCs and released into the field through dynamical ejections or through tidal destruction of the GCs (see MKZ03 and references within; see also Grindlay 1988). Given that both the tidal destruction rate and the metallicity are likely to be correlated with the mass, applying this test is not straightforward. On the other hand, the field sources of elliptical galaxies should show a metallicity effect on their energy spectra regardless of concerns over formation processes.

A better test might then be to extend the suggestion of ITF97 that the difference in the ratio of LMXBs to HMXBs in the Milky Way and the LMC is due to the metallicity difference. Since most HMXBs are accretion powered pulsars, a reasonably good separation between the bright ends of the luminosity distributions of HMXBs and LMXBs should be possibly given good Chandra spectra of nearby spiral galaxies. We note that this is a generic prediction of any model in which the metallicity effects are strictly due to metallicity, but will provide a way to distinguish between true metallicity effects and effects of metallicity being correlated with more difficult to measure parameters, such as the dynamics of the system.

Past work on detailed spectral fitting on Milky Way GCs has been suggested to indicate that there is little evidence for intrinsic absorption is these systems (Sidoli et al. 2001). We note that this is not in conflict with our model. Sidoli et al (2001) did not obtain a satisfctory spectral fit to the data for M 15, the most metal poor of the Milky Way’s globular clusters with an X-ray source, probably because BeppoSax was not capable of resolving the two bright X-ray sources in the cluster (White & Angelini 2001). The other two most metal poor clusters, NGC 1851 and NGC 6712 do show excess absorption in the X-rays compared with the optical, and the other globular clusters in the Milky Way all have optical extinctions significantly higher than the excess predicted by our model, so the fits would not be very sensitive to intrinsic absorption. Additionally, we note that the fitting of Sidoli et al. (2001) was done using the standard assumption that the absorber would be cold material of solar composition. Since BeppoSax is sensitive to absorption edges in the $``$ few keV range, this may cause a systematic error in the fitted absorption value, as noted above. The results of Sidoli et al. (2001) certainly place upper limits on the amount of intrinsic absorption in the Milky Way’s LMXBs, but these upper limits are mostly too high to place strong constraints on our model. A more sensitive test would come from XMM-Newton spectra of M31 globular cluster X-ray sources, where there will be little non-intrinsic absorption since the globular clusters will not be viewed through the disk of the Galaxy. In fact, many of the M31 globular cluster X-ray sources show evidence for intrinsic absorption, and the ones that do are predominantly in metal poor clusters (Irwin & Bregman 1999).

Differences in the luminosity functions between red and blue globular cluster X-ray sources should also, in principle, provide a way to discriminate between models for formation and evolution of their X-ray binaries. Unfortunately, it is difficult at this time to make a prediction from our scenario. It is not clear on theoretical grounds whether the extra wind energy in metal poor systems manifests itself as a higher mass loss rate, yielding probably slightly higher luminosities, albeit for much shorter amounts of time, or as higher wind velocities, in which case the efficiency of wind capture is lower, so the luminosity will be lower at a given mass loss rate, or as some combination of the two. Furthermore, there is not yet a sufficiently large sample of X-ray binaries in globular clusters for making a good comparison of the luminosity functions. This does remain a good test to bear in mind for future work.

## 6 Conclusions

We have outlined a scenario whereby the two metallicity effects seen in LMXBs in globular clusters, higher number density in metal rich clusters, and harder low energy X-ray spectral in metal poor clusters, can be explained via the same mechanism – irradiation induced stellar winds. We have presented additional feasible observational tests of this picture. While we have shown that the physics of the irradiation induced winds required to reproduce the observations is consistent with the most recent theoretical work, we also note that this is a rather complicated problem which is deserving of considerable additional attention by experts in binary stellar structure and evolution. We hope this paper will help to stimulate such work in the future.

## 7 Acknowledgments

We acknowledge funding from NASA LTSA grants NAG5-12975 and NAG5-11319. We thank Michael Sipior, Simon Portegies Zwart and Jasinta Dewi for useful comments on this manuscript. We thank the referee, Josh Grindlay, for suggesting that we improve the clarity of the paper and the depth of the background discussion.


text:
# Condensation Transitions in Nonequilibrium systems

## 1 Introduction

Broadly speaking, one can consider two types of nonequilibrium systems: those relaxing towards thermal equilibrium and those held far from thermal equilibrium e.g. by the system being driven by some external field. In the latter case the steady state of the system will not be described by usual Gibbs-Boltzmann statistical weights rather it will be a nonequilibrium steady state. A natural way to construct a nonequilibrium steady state is to drive the system by forcing a current of some conserved quantity, for example energy or mass, through the system. Such systems are known as driven diffusive systems (DDS) .

In recent years the possibility of phase transitions and phase separation in one-dimensional nonequilibrium systems has been explored and some examples are by now well studied. To appreciate the significance one should recall the general dictum that in one-dimensional equilibrium systems phase ordering and phase transitions do not occur (except in the limit of zero-temperature, or with long range interaction) .

Let us briefly review work on one-dimensional phase transitions in driven systems. A very simple one-dimensional driven diffusive system is the asymmetric simple exclusion process (ASEP). Here particles hop in a preferred direction on a one-dimensional lattice with hard-core exclusion (at most one particle can be at any given site). Indicating the presence of a particle by a 1 and an empty site (hole) by 0 the dynamics comprises the following exchanges at nearest neighbour sites

$`\mathrm{1\; 0}`$ $``$ $`\mathrm{0\; 1}\text{with rate}1`$
$`\mathrm{0\; 1}`$ $``$ $`\mathrm{1\; 0}\text{with rate}q`$ (1)
The open system was studied by Krug and boundary induced phase transitions shown to be possible. Specifically one considers a lattice of $`N`$ sites where at the left boundary site (site 1) a particle is introduced with rate $`\alpha `$ if that site is empty, and at the right boundary site (site $`N`$) any particle present is removed with rate $`\beta `$. Thus the dynamical processes at the boundaries are

$`\text{at site }10`$ $``$ $`1\text{with rate}\alpha `$
$`\text{at site }N1`$ $``$ $`0\text{with rate}\beta .`$ (2)
These boundary conditions force a steady state current of particles $`J`$ through the system. Phase transitions occur when $`lim_N\mathrm{}J`$ exhibits non-analyticities. The steady state of this system was solved exactly for the totally asymmetric case  and more recently for the general $`q`$ case . When $`q<1`$ the phase diagram comprises three phases: a high-density phase where the current is controlled by a low exit rate $`\beta `$—one can think of this is queue of cars at a traffic light that doesn’t let many cars through; a low-density phase where the current controlled by a low injection rate $`\alpha `$—think of this as a traffic light that does not let many cars onto an open road; a maximal-current phase where both $`\alpha ,\beta `$ are high ($`\alpha ,\beta >(1q)/2`$) and the current is $`J=(1q)/4`$. Note that since increasing $`\alpha `$ and $`\beta `$ doesn’t increase the current, the current is saturated. In the maximal current phase generic long-range correlations exist, an example being the decay of particle density from the left boundary to the bulk value $`1/2`$ which is a power law $`1/x^{1/2}`$ where $`x`$ is distance from the left boundary.

On the line $`\alpha =\beta <(1q)/2`$ which separates the high and low density phases one finds coexistence between a region of low density in the left part of lattice and a region of high density on the right separated by a ‘shock’ where the density changes sharply over a microscopic distance.

Perodic systems (i.e. a ring of sites) can also exhibit phase separation when inhomogeneities or defects are introduced. A very simple example is to introduce into the asymmetric exclusion process a ‘slow bond’ through which particles hop with a reduced rate. Then in the steady state one can obtain phase separation between a region of high density behind the slow bond and a region of low density in front of the slow bond. Moving defects (i.e. particles with dynamics different from that of the others) have also been considered and exact solutions obtained . One can think of a slow agricultural vehicle on a country road with a large queue of cars behind it and open road in front of it.

A further question is whether systems related to the hopping particle models described so far, but without inhomogeneities, can exhibit phase ordering. A very simple model was introduced in  comprising three species of conserved particles, amongst which all possible exchanges are allowed. However a key feature is that the dynamics has a cyclic symmetry i.e. $`A`$ particles move preferentially to the left of $`B`$ particles which move preferentially to the left of $`C`$ particles which in turn move preferentially to the left of $`A`$ particles. The model exhibits strong phase separation into pure domains of $`A`$ $`B`$ $`C`$. Similar strong phase separation occurs in other related models .

A final class of transitions in one-dimensional hopping particle models is that involving spatial condensation, whereby a finite fraction of the particles condenses onto the same site. Examples include the appearance of a large aggregate in models of aggregation and fragmentation and the emergence of a single flock in dynamical models of flocking . We will analyse a simple example of a condensation transition which occurs in the zero-range process which we now define.

## 2 The zero-range process

The zero-range process was introduced by Spitzer  and recent applications and developments have been reviewed in . We consider a one-dimensional lattice of $`M`$ sites with sites labelled $`\mu =1\mathrm{}M`$ and periodic boundary conditions (more generally one can consider the zero-range process on a lattice of arbitrary dimension). Each site can hold an integer number of indistinguishable particles. The configuration of the system is specified by the occupation numbers $`n_\mu `$ of each site $`\mu `$. The total number of particles is denoted by $`L`$ and is conserved under the dynamics. The dynamics of the system is given by the rates at which a particle leaves a site $`\mu `$ (one can think of it as the topmost particle—see Figure 1a) and moves to the left nearest neighbour site $`\mu 1`$. The hopping rates $`u(n)`$ are a function of $`n`$ the number of particles at the site of departure. Some particular cases are: if $`u(n)=n`$ then the dynamics of each particle is independent of the others; if $`u(n)=\mathrm{const}`$ for $`n>0`$ then the rate at which a particle leaves a site is unaffected by the number of particles at the site (as long as it is greater than zero).

The important attribute of the zero-range process is that it has a steady state described by a product measure. By this it is meant that the steady state probability $`P(\{n_\mu \})`$ of finding the system in configuration $`\{n_1,n_2\mathrm{}n_M\}`$ is given by a product of factors $`f(n_\mu )`$

$$P(\{n_\mu \})=\frac{1}{Z(M,L)}\underset{\mu =1}{\overset{M}{}}f(n_\mu ).$$
(3)
Here the normalisation $`Z(M,L)`$ is introduced so that the sum of the probabilities for all configurations, with the correct number of particles $`L`$, is one.

In the basic model described above, $`f(n)`$ is given by

$`f(n)`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{n}{}}}{\displaystyle \frac{1}{u(m)}}\text{for}n1`$ (4)
$`=`$ $`1\text{for}n=0`$
To prove (3,4) one simply considers the stationarity condition on the probability of a configuration (probability current out of the configuration due to hops is equal to probability current into the configuration due to hops):

$$\underset{\mu }{}\theta (n_\mu )u(n_\mu )P(n_1\mathrm{}n_\mu \mathrm{}n_L)=\underset{\mu }{}\theta (n_\mu )u(n_{\mu +1}+1)P(n_1\mathrm{}n_\mu 1,n_{\mu +1}+1\mathrm{}n_L).$$
(5)
The Heaviside function $`\theta (n_\mu )`$ highlights that it is the sites with $`n>1`$ that allow exit from the configuration (lhs of (5)) but also allow entry to the configuration (rhs of (5)). Equating the terms $`\mu `$ on both sides of (5) and cancelling common factors assuming (3), results in

$$u(n_\mu )f(n_{\mu 1})f(n_\mu )=u(n_{\mu +1}+1)f(n_\mu 1)f(n_{\mu +1}+1)$$
(6)
This equality can be recast as

$$u(n_\mu )\frac{f(n_\mu )}{f(n_\mu 1)}=u(n_{\mu +1}+1)\frac{f(n_{\mu +1}+1)}{f(n_{\mu +1})}=\text{constant}$$
(7)
Setting the constant equal to unity implies

$$f(n_\mu )=\frac{f(n_\mu 1)}{u(n_\mu )}$$
(8)
and iterating (8) leads to (4).

We can easily generalise to consider an inhomogeneous system by which we mean the hopping rates are site dependent: the hopping rate out of site $`\mu `$ when it contains $`n_\mu `$ particles is $`u_\mu (n_\mu )`$. It is easy to check that the steady state is simply modified to

$$P(\{n_\mu \})=\frac{1}{Z(M,L)}\underset{\mu =1}{\overset{L}{}}f_\mu (n_\mu )$$
(9)
where $`f_\mu `$ are given by

$`f_\mu (n)`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{n}{}}}{\displaystyle \frac{1}{u_\mu (m)}}\text{for}n1`$ (10)
$`=`$ $`1\text{for}n=0`$
The proof is identical to that for the homogeneous case, with the replacement of $`u(n_\mu )`$ by $`u_\mu (n_\mu )`$

There exists an exact mapping from a zero-range process to an asymmetric exclusion process. This is illustrated in Figure 1. The idea is to consider the particles of the zero-range process as the holes (empty sites) of the exclusion process. Then the sites of the zero-range process become the moving particles of the exclusion process. Note that in the exclusion process we have $`M`$ particles hopping on a lattice of $`M+L`$ sites. A hopping rate in the zero range process $`u(m)`$ which is dependent on $`m`$ corresponds to a hopping rate in the exclusion process which depends on the gap to the particle in front. So the particles can feel each other’s presence and one can have a long-range interaction.

## 3 Condensation Transitions

We now proceed to analyse the steady states of form (9) and the condensation transition that may occur. The important quantity to consider is the normalisation $`Z(M,L)`$ as it plays the role of the partition sum. The normalisation is defined through the condition

$$Z(M,L)=\underset{n_1,n_2\mathrm{}n_M}{}\delta (\underset{\mu }{}n_\mu L)\underset{\mu =1}{\overset{M}{}}f_\mu (n_\mu )$$
(11)
where the $`\delta `$ function enforces the constraint of $`L`$ particles. The normalisation may be considered as the analogue of a canonical partition function of a thermodynamic system.

We define the ‘speed’ $`v`$ as the average hopping rate out of a site

$`v`$ $`=`$ $`{\displaystyle \frac{1}{Z(M,L)}}{\displaystyle \underset{n_1,n_2\mathrm{}n_M}{}}\delta ({\displaystyle \underset{\mu }{}}n_\mu L)u(n_1){\displaystyle \underset{\mu =1}{\overset{M}{}}}f_\mu (n_\mu )`$ (12)
$`=`$ $`{\displaystyle \frac{Z(M,L1)}{Z(M,L)}}`$
where we have used (9,10). Note that (12) tells us that the speed is independent of site and thus may be considered a conserved quantity in the steady state of the system. In the totally asymmetric system considered in Section 2 the speed is equal to the current of particles flowing between neighbouring sites. The speed is a ratio of partition functions of different system sizes (12) and corresponds to a fugacity.

We now use the integral representation of the delta function to write the partition function as

$$Z(M,L)=\frac{dz}{2\pi i}z^{(L+1)}\underset{\mu =1}{\overset{M}{}}F_\mu (z),$$
(13)
where

$$F_\mu (z)=\underset{m=0}{\overset{\mathrm{}}{}}z^mf_\mu (m).$$
(14)
For large $`M,L`$ (13) is dominated by the saddle point of the integral and the value of $`z`$ at the saddle point is the fugacity. The equation for the saddle point reduces to

$$\frac{L}{M}=\frac{z}{M}\underset{\mu =1}{\overset{M}{}}\frac{}{z}\mathrm{ln}F_\mu (z)$$
(15)
which, defining $`\varphi =L/M`$, can be written as

$$\varphi =\frac{z}{M}\underset{\mu =1}{\overset{M}{}}\frac{F_\mu ^{}(z)}{F_\mu (z)}.$$
(16)
In the thermodynamic limit,

$$M\mathrm{}\text{with}L=\varphi M,$$
(17)
where the density $`\varphi `$ is held fixed, the question is whether a valid saddle point value of z can be found from (16). We expect that for low $`\varphi `$ the saddle point is valid but, as we shall discuss, there exists a maximum value of $`z`$ and if at this maximum value the rhs of (16) is finite, then for large $`\varphi `$ (16) cannot be satisfied. We now consider how condensation may occur in the inhomogeneous and the homogeneous case.

### 3.1 Inhomogeneous case

To give an idea of how a condensation transition may occur we consider the case $`u_\mu (m)=u_\mu `$ for $`m>0`$ i.e. the hopping rate does not depend on the number of particles at a site. $`f_\mu `$ is given by

$$f_\mu (n)=\left(\frac{1}{u_\mu }\right)^{n_\mu }$$
(18)
and the probability of occupancies $`\{n_1,n_2,\mathrm{},n_M\}`$ is

$$P(\{n_1,n_2,\mathrm{},n_M\})=\frac{1}{Z(M,L)}\underset{\mu =1}{\overset{M}{}}\left(\frac{1}{u_\mu }\right)^{n_\mu }.$$
(19)
The mapping to an ideal Bose gas is evident: the $`L`$ particles of the zero-range process are viewed as Bosons which may reside in $`M`$ states with energies $`E_\mu `$ determined by the site hopping rates: $`\mathrm{exp}(\beta E_\mu )=1/u_\mu `$. Thus the ground state corresponds to the site with the lowest hopping rate. The normalisation $`Z(M,L)`$ is equivalent to the canonical partition function of the Bose gas. We can sum the geometric series (14) to obtain $`F_\mu `$ and $`F_\mu ^{}`$ then taking the large $`M`$ limit allows the sum over $`\mu `$ to be written as an integral

$$\varphi =_{u_{\mathrm{min}}}^{\mathrm{}}𝑑u𝒫(u)\frac{z}{uz}$$
(20)
where $`𝒫(u)`$ is the probability distribution of site hopping rates with $`u_{\mathrm{min}}`$ the lowest possible site hopping rate. Interpreting $`𝒫(u)`$ as a density of states, equation (20) corresponds to the condition that in the grand canonical ensemble of an ideal Bose gas the number of Bosons per state is $`\varphi `$. The theory of Bose condensation tells us that when certain conditions on the density of low energy states pertain we can have a condensation transition. Then (16) can no longer be satisfied and we have a condensation of particles into the ground state, which is here the site with the slowest hopping rate.

A very simple example is to have just one ‘slow site’ i.e. $`u_1=p`$ while the other $`M1`$ sites have hopping rates $`u_\mu =1`$ when $`\mu >1`$. Using the mapping to an exclusion process, this corresponds one slow particle i.e. agricultural vehicle example described earlier. One can show  that for a high density of particles in the zero range process (low density of particles in the corresponding asymmetric exclusion process) we have a condensate since site 1 contains a finite fraction of the particles. In the low density phase the particles are evenly spread between all sites.

### 3.2 Homogeneous case

We now consider the homogeneous zero-range process where the hopping rates $`u(n)`$ are site independent. Then (14) is independent of $`\mu `$ and reads

$$F(z)=\underset{n=0}{\overset{\mathrm{}}{}}\underset{m=1}{\overset{n}{}}\left[\frac{z}{u(m)}\right]$$
(21)
The fugacity $`z`$ must be chosen so that $`F`$ converges or else we could not have performed (14). Therefore $`z`$ is restricted to $`z\beta `$ where we define $`\beta `$ to be the radius of convergence of $`F(z)`$. From (21) we see that $`\beta `$ is the limiting value of the $`u(m)`$ i.e. the limiting value of the hopping rate out of a site for a large number of particles at a site. We interpret (16) as giving a relation between the density of holes (number of holes per site) and the fugacity $`z`$. The saddle point condition (16) becomes

$$\varphi =\frac{zF^{}(z)}{F(z)}$$
(22)
Given that the rhs of (22) is a monotonically increasing function of $`z`$ we deduce that density of particle increases with fugacity. However if at $`z=\beta `$, the maximum allowed value of $`z`$, the rhs of (22) is still finite then one can no longer solve for the density and one must have a condensation transition. Physically, the condensation would correspond to a spontaneous symmetry breaking where one of the sites is spontaneously selected to hold a finite fraction of the particles.

Thus, for condensation to occur (i.e. when $`\varphi `$ is large enough for (22) not to have a solution for the allowed values of $`z`$) we require

$$\underset{z\beta }{lim}\frac{F^{}(z)}{F(z)}<\mathrm{}.$$
(23)
We now assume that $`u(n)`$ decreases uniformly to $`\beta `$ in the large $`n`$ limit as

$$u(n)=\beta (1+\zeta (n))$$
(24)
where $`\zeta (n)`$ is a monotonically decreasing function. Analysis of the series

$`F(\beta )`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{exp}\left\{{\displaystyle \underset{m=1}{\overset{n}{}}}\mathrm{ln}\left[1+\zeta (m)\right]\right\}`$
$`F^{}(\beta )`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}n\mathrm{exp}\left\{{\displaystyle \underset{m=1}{\overset{n}{}}}\mathrm{ln}\left[1+\zeta (m)\right]\right\}`$ (25)
reveals that the condition for condensation is simply that $`F^{}(\beta )`$ is finite and this occurs if $`u(n)`$ decays to $`\beta `$ more slowly than $`\beta (1+2/n)`$. (This is easiest to see by expanding $`\mathrm{ln}\left[1+\zeta \right]`$ and approximating the sum over $`m`$ by an integral in (25).)

It is interesting to translate this result into the language of the exclusion process. In this context we can have condensation if the hop rate of a particle into a gap of size $`n`$ decays as $`\beta (1+2/n)`$ therefore there is an effective long range interaction.

### 3.3 Bus route model

As an example of this let us consider the ‘bus route model’ . The model is defined on a $`1d`$ lattice. Each site (bus-stop) is either empty, contains a bus (a conserved particle) or contains a passenger (non-conserved quantity). The dynamical processes are that passengers arrives at an empty site with rate $`\lambda `$; a bus moves forward to the next stop with rate 1 if that stop is empty; if the next stop contains passengers the bus moves forward with rate $`\beta `$ and removes the passengers.

The model thus defined has not been solved but simulations reveal two regimes. At high bus density the gaps between buses are evenly distributed. However at low bus density there is a condensed regime where the lead bus has a large gap to the next bus in front of it with bus-stops full of passengers in between. The other buses have small gaps between them. Thus the buses form a jam of buses and after a long delay all arrive at a bus-stop at once.

The bus route model can be related to the zero-range process by a mean-field approximation in which we integrate out the non-conserved quantity (passengers). The idea is that a bus-stop, next to bus 1 say, will last have been visited by a bus (bus 2) a mean time ago of $`n/v`$ where $`n`$ is the distance between bus 2 and bus 1 and $`v`$ is the steady state speed. Therefore the mean-field probability that the site next to bus 1 is not occupied by a passenger is $`\mathrm{exp}(\lambda n/v)`$. From this probability an effective hopping rate for a bus into a gap of size $`n`$ is obtained by averaging the two possible hop rates $`1,\beta `$:

$$u(n)=\beta +(1\beta )\mathrm{exp}(\lambda n/v).$$
(26)
We can now see that this mean-field approximation to the bus-route model is equivalent to a homogeneous zero-range process discussed earlier.

Since $`u(n)`$ decays exponentially the condition for a strict phase transition in the thermodynamic limit is not met. However on any finite system for $`\lambda `$ sufficiently small, an apparent condensation will be seen. In the bus route problem this corresponds to the universally irritating situation of all the buses on the route arriving at once.

## 4 Conclusion

We have shown how the zero range process exhibits two kinds of condensation transition. One is due to having an inhomogeneous system i.e. we get condensation of particles onto the site with the slowest hopping rate. Although the condensation is spatial the mechanism is equivalent to Bose condensation in an Ideal Bose Gas.

The other type of condensation occurs on a homogeneous systems and involves the spontaneous selection of a site onto which a finite fraction of the particles condense. Recently this condensation mechanism has been used to understand the existence or non-existence of phase separation in a general class of one dimensional driven systems .


text:
# Dislocation nucleation in shocked fcc solids: effects of temperature and preexisting voids

## Abstract

Quantitative behaviors of shock-induced dislocation nucleation are investigated by means of molecular dynamics simulations on fcc Lennard-Jones solids: a model Argon. In perfect crystals, it is found that the Hugoniot elastic limit (HEL) is a linearly decreasing function of temperature: from near-zero to melting temperatures. In a defective crystal with a void, dislocations are found to nucleate on the void surface. Also HEL drastically decreases to $`15`$ percent of the perfect crystal when the void radius is $`3.4`$ nanometer. The decrease of HEL becomes larger as the void radius increases, but HEL becomes insensitive to temperature.

Mechanical properties of shock-loaded solids are important to materials science since they can reveal reliability of materials under extreme conditions. Shocks in solids generate high pressure almost instantaneously, and create extremely high deformation (strain) rates. Yielding phenomena realized at high strain rates are usually quite different from those of quasistatic ones. Notably, Rohde has found that the HEL is almost independent of temperature rohde ; Kanel and his coworkers have found that the HELs of Al and Cu single crystals are increasing functions of temperature kanel . These experiments are in apparent contrast with quasistatic deformations, where yield strength considerably decreases with increasing temperature. However, contrary to Rohde and Kanel, stainless steel shows that the HEL is a decreasing function of temperature gu . These puzzling results provide us with intriguing problems. The microscopic foundations of plastic deformation are dominated by the dynamics of dislocations. Unfortunately, the nature of dislocation dynamics is very complicated and yet to be fully understood from theoretical point of view. Therefore, it is reasonable to decompose the problem into three essential ingredients of dislocation dynamics: nucleation, mobility, and multiplication. As the first stage of this line of thought, the properties of dislocation nucleation are investigated in this paper.

In a perfect crystal, dislocations nucleate with the help of thermal fluctuations. Holian and Lomdahl have studied the emergence of stacking faults initiated by partial dislocation emission in a shocked perfect crystal holian . Recently, Tanguy et al. tanguy  have performed extensive simulations to confirm that there exists a critical size of dislocation loop below which nucleated loops are energetically unstable and eventually annihilate; this is reminescent of droplet nucleation in a metastable gas or liquid. Although these two studies are enlightening, the role of temperature, which is known to be important to dislocation nucleation, has not been considered. In some studies, temperature plays no role in activating dislocation nucleation xu . However, dislocation nucleation seems to be anyhow a fluctuation-assisted process in the sense that no nucleation is observed in a molecular dynamics simulation of a perfect crystal at zero temperature holian . It is important to understand temperature effects on dislocation nucleation to clarify the controversial experimental results rohde ; kanel ; gu .

This paper discusses temperature effects on dislocation nucleation through molecular dynamics simulations. It turns out that the HEL is a linearly decreasing function of temperature from zero to melting temperatures. Then a defect is introduced into the system as a preexisting void, at which dislocation nucleation is found to be enhanced.

The present simulations use a simple Lennard-Jones (LJ) potential $`U(r)=4ϵ\left[(\sigma /r)^{12}(\sigma /r)^6\right]`$ with a cutoff length of $`4.0\sigma `$. To avoid a discontinutity at the cutoff, the potential is suitably shifted. Note that simulations with shorter cutoffs such as $`2.6\sigma `$ resulted in a visible decrease of HEL, as much as $`90`$ percent. Hereafter, we let $`\sigma =3.41`$ Åand $`ϵ=1.65\times 10^{21}`$ \[J\] in order to make a quantitative comparison with experiments on solid Argon. The Lennard-Jones potential has been well-tested in terms of mechanical and thermal properties of Argon such as the elastic constants and the melting temperature. Hence it is reasonable to adopt the LJ potential as the first computational attempt to understand the anomalous HEL behavior of fcc metals kanel .

In our system the axes of $`x`$, $`y`$, and $`z`$ are taken to coincide with $`100`$, $`010`$, and $`001`$ directions, respectively. Shocks propagate along the $`[001]`$ direction. As usual, periodic boundary conditions are applied to $`x`$ and $`y`$ directions. The whole system consists of $`100\times 100\times 200`$ unit cells. Shock waves in the present simulations are generated by pushing a piston of infinite mass into the still target. The velocity of the piston is denoted by $`u_p`$. We call the piston velocity $`u_p`$ throughout this paper.

In this work, the HEL is defined as the longitudinal stress ($`\sigma _{zz}`$) above which dislocation nucleation is observed. Note that this definition of HEL differs from the one that is based entirely on shockwave behavior: the shock front decomposes into an elastic precursor and a plastic wave. However, these two definitions are identical because the decomposition is also observed in molecular dynamics simulations once dislocations are emitted in shocked solids. Longitudinal stress is calculated from the Hugoniot relation; $`\sigma _{zz}=\rho _0u_su_p`$, where $`\rho _0`$ and $`u_s`$ denote initial density and shock velocity, respectively. Note that the initial pressure is zero. The longitudinal stress $`\sigma _{zz}`$ is approximately proportional to the shear stress holian  which is directly responsible for dislocation emission.

Figure. 1 shows the nearly linear dependence of HEL on temperature. The HEL at $`77`$ K, close to the melting temperature ($`83`$ K) agrawal , is approximately half that $`6`$ K. This is clear evidence that dislocation nucleation is a thermally assisted phenomenon. Also, this temperature dependence is reminiscent of elastic constants. The dashed line in FIG. 1 represents the shear modulus of solid Argon keeler . We can see that the temperature dependences of the HEL and of the shear modulus are qualitatively the same for perfect crystals. That is, HEL is approximately half the shear modulus regardless of temperature. (Indeed, the ratio of HEL to shear modulus is slightly greater than $`0.5`$ for $`T24`$ K.) The critical strain also decreases linearly from $`0.138`$ (at $`6`$ K) down to $`0.118`$ (at $`77`$ K) as the temperature increases. Numerical data are shown in Table 1. Note that the temperature dependence of the shear modulus of Argon has been calculated by molecular dynamics simulation using LJ potential shimizu .

The temperature dependence of HEL obtained by the present simulation is opposite to the experiment on fcc metals kanel . But this is not a contradiction, for the experiment involves real crystals where various defects and impurities preexist. Plastic deformations of real materials are mainly governed by the motion of preexisting dislocations, while in a perfect crystal only the dislocation nucleation is considered here. Hence, the mobility of preexisting dislocations must be analyzed for further understanding of temperature dependence of shock-induced plasticity. This subject will be discussed elsewhere.

Strikingly different results are found when we turn to defective crystals, concentrating on the effects of a void. A void is introduced in a model crystal as a blanked spherical region. The radius is varied from $`0.34`$ nm to $`3.41`$ nm. With a void, dislocations nucleate exclusively on the void surface as shown in FIGS. 2, 3, and 4.

The Hugoniot elastic limit decreases drastically in the presence of a void. Square symbols in FIG. 1 represent the HEL for a defective crystal with a void of $`1.71`$ nm radius. With this size of void, the HEL is approximately $`20`$ percent of the perfect crystal’s. Note that the HEL of a defective crystal also linearly decreases with increasing temperature up to $`42`$ K.

When $`T>42`$ K, approximately half the melting temperature, no dislocation is nucleated and a void collapses after the passage of shock front. Although this phenomenon is similar to ”hotspot” formation hatano , this is not the case since the piston velocity used here is not strong enough to cause a jet which is essential to hotspot. The collapse of the void observed in the present situation rather postulates that voids are unstable to perturbations above half the melting temperature gittus .

Below $`42`$ K, the HEL and the critical piston velocity for the defective crystal with a void become insensitive to temperature. For the case of $`r=1.02`$ nm, the critical piston velocity is $`103\pm 5`$ m/s at $`6`$ K, dropping only to $`95\pm 5`$ m/s at $`42`$ K, while for the perfect crystal the corresponding values are $`289\pm 1.6`$ m/s and $`213\pm 1.6`$ m/s, respectively. The absolute decreases are $`8`$ m/s for the defective crystal and $`76`$ m/s for the perfect crystal. Then the relative changes are $`8/1030.08`$ and $`76/2840.26`$, respectively. The difference between the perfect and the defective crystals is remarkable in terms of both absolute and relative decreases. Thus defects dominate thermal fluctuations in the enhancement of dislocation nucleations. In experiment, the microstructure of the specimen is important.

The void-size dependence of HEL is shown in Fig. 5. At small $`r`$, the HEL decreases as $`r^1`$ followed by a crossover to $`r^{0.5}`$ at $`r_c=1.3\pm 0.3`$ nm. Note that the transition radius $`r_c`$ is approximately the width of a double kink and is also very close to the critical loop radius $`1.5`$ nm above which nucleated dislocation loops grow stably tanguy . The quantitative explanation of this intriguing transition is not clear at this point.

To compare the present simulation with experiment, the Hugoniot values ($`u_p`$ versus $`u_s`$) of Argon for both this simulation and experiment dick  are presented in TABLE 2. We see that the simulation values exceed the experimental ones by approximately $`15`$ percent. This overestimation is comparable to that of Belonoshko’s simulations using a Buckingham potential belonoshko . His deviation from the experimental value exceeds $`10`$ percent, especially for small piston velocities up to $`1.2`$ km/sec. Dubrovinsky has attributed this overestimation to the polycrystalline nature of the experimental specimen dubrovinsky . Since both Belonoshko’s and the present simulations deal only with crystals of $`100`$ orientation, the deviation from the experimental result is not unreasonable.

Finally, two quantitative uncertainties in the present simulation are remarked. One is underestimation of stacking fault energy by the use of LJ potential. It is known that stacking fault energy vanishes when one adopts short-range, two-body intermolecular potentials dubrovinsky . In this regard, it is possible that the numerical values of the HEL obtained by the present simulations are underestimates. But, at least, qualitative results such as the temperature dependence or the void effect are not so influenced by the nature of two-body potentials. In addition, a large cut-off, $`4.0\sigma `$, was used to reduce this effect.

The other uncertainty invloves periodic boudary conditions. The system size ($`57\times 57\times 114`$ nm) is small enough that a dislocation loop may interact with itself through the periodic boudary conditions. It is impossible to remove this uncertainty from molecular dynamics simulations at this point. However, the strain around a nucleated dislocation loop propagates at the sound velocity (about $`900`$ m/s for transverse mode), taking $`45`$ picoseconds to cross the simulation cell. Therefore, this efect may not seriously influence the present results, which involve exclusively nucleation processes that occur within $`10`$ picoseconds.

To summarize, molecular dynamics simulations on both perfect and defective fcc crystals are performed from near-zero to melting temperatures. It is found that the Hugoniot elastic limit is about half the shear modulus and decreases linealy with increasing temperature. The critical piston velocity and critical strain also decrease in the same manner. When a void is introduced in the crystal, dislocations are found to nucleate on the void surface. The critical piston velocity decreases drastically as the void radius $`r`$ increases. The decrease is nonlinear and shows a crossover from $`r^1`$ to $`r^{0.5}`$ around $`r1.3`$ nm, which is approximately the width of a double kink. The Hugoniot elastic limit becomes insensitive to temperature in the presence of a void. That is, the effect of nucleation at void dominates that of thermal fluctuations in shock-induced dislocation nucleation.

The author gratefully acknowledges helpful discussions with H. Kaburaki, F. Shimizu, Y. Satoh, N. Nita, and H. Matsui.