text:
# Effective forces in colloidal mixtures: from depletion attraction to accumulation repulsion

## I Introduction

A fundamental description of colloidal interactions based on statistical mechanics is needed to understand and predict the stability and the phase behavior of colloidal suspensions. Most inter-colloidal forces are effective in the sense that some microscopic degrees of freedom are averaged out. This concept of mean or effective interactions  is crucial to bridge the different length scales involved in colloidal systems and has been exploited in many different circumstances. Examples of the microscopic degrees of freedom include solvent particles, smaller colloidal particles, added polymer coils  or monomers of (grafted) polymer chains as well as counter-ions and salt ions in the case of charged suspensions. The resulting effective forces turn out to exhibit a wide range of features: They can be attractive, repulsive, or oscillatory, and are an important key to understanding colloidal stability as well as flocculation and coagulation .

In the present paper we consider a binary colloidal mixture of big and small colloidal particles and investigate the distance-resolved effective force and potential between two big colloidal spheres surrounded by many small colloidal spheres, as depicted in Fig. 1. The small particles are averaged out so that one is left with an effective pair potential between the big ones, which is superimposed onto the direct big-big interaction. The form of these effective interactions depends sensitively on the basic big-small and small-small interactions.

Most previous studies have focused on a hard-sphere (HS)-like interaction between the colloidal particles, with a special emphasis on two cases: (1) the Asakura-Oosawa model, which was originally designed to describe mixtures of sterically-stabilized colloids and nonadsorbing polymer coils . In this strongly non-additive model, the big-small interaction is HS like, while the small-small interaction is ideal gas like. The effective pair potential between two colloids can be exactly calculated and has an attractive well proportional to the density of the small particles, and a range equal to their diameter. (2) Much work has also focused on the complementary model of additive HS mixtures, inspired in part by a prediction that they might phase-separate. The effective depletion forces were obtained by computer simulation , approximate theories, density functional theory (DFT) , and experiments. Simplified potentials were used to investigate the phase-behavior of binary HS mixtures, where it was found that fluid-fluid phase-separation, when it occurs, was always metastable w.r.t. fluid-solid phase-separation.

In this work we generalize these studies to arbitrary big-small and small-small interactions. Our study is motivated by the fact that real colloidal suspensions typically have soft interactions which are beyond the HS model. These soft interactions can easily be tuned, for example by changing the solvent composition, the amount of added salt, the surface charge, etc…, which provides a route to tailoring the effective interactions between the big particles. Hence it is essential to understand on a full statistical level how the basic big-small and small-small interactions affect the effective big-big interactions in a binary colloidal mixture.

There has been some previous work on how non-HS like interactions affect the effective big-big interactions, see e.g. refs. , which showed that the differences with the pure HS case could be substantial. However, with a few exceptions, all these studies were performed with approximate methods such as integral equations. Our goal here is to provide benchmark exact computer simulation results for a set of $`9`$ different parameters corresponding to combinations of attractive and repulsive small-small and big-small interactions. This allows us to systematically study the trends in the effective forces as induced by the basic interactions. It also allows us to test several theoretical techniques, namely DFT, the superposition approximation, and a mapping to a non-additive HS model.

Generalizing from an additive to a non-additive HS mixture already leads to a much richer class of effective potentials, which can be exactly calculated from an accurate theory of additive HS mixtures. Similarly, we find here that the effective big-big interactions are profoundly influenced by the basic big-small and small-small interactions. The trends can be summarized as follows: Adding a big-small repulsion results in enhanced attraction through the standard depletion mechanism. Adding a big-small attraction leads to an accumulation of the small particles near each big one, which in turn results in a more repulsive big-big interaction (accumulation repulsion). For a given big-small interaction, adding a small-small repulsion also results in an enhanced small particle density near the big particles, and therefore in a more repulsive big-big interaction. Adding a small-small attraction when the big-small interaction is repulsive results in a more attractive effective interaction. All these effects can be qualitatively understood from a mapping to a non-additive HS Model. But when a small-small and a big-small attraction are combined, a repulsion through attraction effect occurs, which is not captured by DFT or by our mapping scheme.

The paper is organized as follows: In Chapter II, we discuss how to formally map a two-component mixture onto an effective one-component one, and define our target quantities. In Chapter III, we describe our model for the basic interactions. We discuss our simulation results in Chapter IV, and the results of several different theories in Chapter V. Our conclusions are stated in Chapter VI.

## II Mapping a binary mixture onto an effective one-component system

In this section we briefly describe how to map a two-component binary mixture onto an effective one-component system. This statistical mechanical procedure was originally developed by McMillan and Mayer to trace out the solvent in a suspension, but it can be used just as well to trace out a smaller component in a binary colloid mixture.

We begin with the Hamiltonian of our two-component mixture

$$H=K+H_{bb}+H_{ss}+H_{bs}.$$
(1)
Here $`K`$ is the total kinetic energy of the mixture and three potential energy contributions are:

$`H_{bb}`$ $`=`$ $`{\displaystyle \underset{i<j}{\overset{N_b}{}}}\mathrm{\Phi }_{bb}(𝐫_i^b𝐫_j^b),`$ (2)
$`H_{ss}`$ $`=`$ $`{\displaystyle \underset{i<j}{\overset{N_s}{}}}\mathrm{\Phi }_{ss}(𝐫_i^s𝐫_j^s),`$ (3)
$`H_{bs}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N_b}{}}}{\displaystyle \underset{j=1}{\overset{N_s}{}}}\mathrm{\Phi }_{bs}(𝐫_i^b𝐫_j^s),`$ (4)
where the $`\mathrm{\Phi }_{ij}`$ denote the pairwise interaction potentials, and $`𝐫_i^b`$ and $`𝐫_j^s`$ denote the coordinates of the centers of the big and small particles respectively.

The formal procedure to reduce the two-component system to an effective one-component one, described in more detail elsewhere, see e.g., consists of tracing out the small particles for a fixed configuration $`\{𝐫_i^b\}`$ of the big particles. This is most conveniently done for a semi-grand ensemble where the small particles in the volume $`V`$ are kept at a fixed chemical potential $`\mu _s`$. The resulting effective one-component Hamiltonian that governs the behavior of the big particles takes the form:

$$H_{bb}^{eff}=H_{bb}+\mathrm{\Omega },$$
(5)
where $`\mathrm{\Omega }=\mathrm{\Omega }(N_b,z_s,V;\{𝐫_i^b\})`$ is the grand potential of a fluid of small particles at fugacity $`z_s=\lambda _s^3\mathrm{exp}(\beta \mu _s)`$, ($`\beta =1/k_BT`$) subjected to the external potential of the fixed configuration $`\{𝐫_i^b\}`$ of the big particles. This grand-potential can be further expanded as a sum of (n-body) terms:

$$\mathrm{\Omega }=\underset{n=0}{\overset{N_b}{}}\mathrm{\Omega }_n.$$
(6)
The first term, $`\mathrm{\Omega }_0`$, is the so-called volume term, which only depends on properties of the small-particles (i.e. $`\mathrm{\Phi }_{ss}`$) and is therefore independent of the configuration $`\{𝐫_i^b\}`$. The one-body term $`\mathrm{\Omega }_1`$ can be related to $`N_b`$ times the free-energy gained by inserting a single big particle into a sea of small particles. It depends on both $`\mathrm{\Phi }_{ss}`$ and $`\mathrm{\Phi }_{bs}`$, but is independent of $`\mathrm{\Phi }_{bb}`$. Of most interest to us here is the two-body term $`\mathrm{\Omega }_2`$, which can be written as:

$$\mathrm{\Omega }_2(N_b,z_s,\{𝐫_i^b\})=\underset{i<j}{\overset{N_b}{}}V_{bb}^{eff}(|𝐫_i𝐫_j|),$$
(7)
where the two-body effective pair-potential $`V_{bb}^{eff}(r)`$ is related to the free-energy difference between two big spheres a distance $`r`$ apart, immersed in a sea of small spheres at fugacity $`z_s`$, and the same system with the two spheres at $`r=\mathrm{}`$. It is important to notice that $`V_{bb}^{eff}(r)`$ itself only depends on $`\mathrm{\Phi }_{bs}`$ and $`\mathrm{\Phi }_{ss}`$, and not directly on $`\mathrm{\Phi }_{bb}`$. This observation was a key to understanding why the effective potentials for non-additive HS mixtures can be calculated by a theory for additive HS mixtures. Nevertheless, through Eq. (5), the effect of $`\mathrm{\Phi }_{bb}(r)`$ is felt in that it sets the range where $`V_{bb}^{eff}(r)`$ is relevant.

In a similar way, higher order interactions can be derived, but for this scheme to remain tractable, one usually truncates at the pair level. As long as the ranges of $`\mathrm{\Phi }_{ss}`$ and $`\mathrm{\Phi }_{bs}`$ are small compared to the range of $`\mathrm{\Phi }_{bb}`$, higher order terms are not expected to be very important.

Thus far we have developed this formal tracing out in the semi-grand ensemble, which helps emphasize that in a mixture with many large particles, the correct effective potential $`V_{bb}^{eff}(r)`$ to be used in Eq. (7) is the one fixed by the chemical potential or fugacity instead of the overall density of the smaller particles, as first emphasized by Lekkerkerker et al.. It is important to keep this in mind when studying phase-behavior. However, for a practical calculation of the pair term, with only two big particles in an infinite sea of small ones, one can just as easily work in the canonical ensemble, keeping in mind that when using the obtained $`V_{bb}^{eff}(r)`$ for a system with many big particles, the input density $`\rho _s`$ is really that of a reservoir of small particles kept at the same chemical potential as the full mixture of big and small particles.

For the case of two big-spheres a distance $`r`$ apart, in a bath of small particles at density $`\rho _s`$, the average force is related to the effective potential by:

$$𝐅_{bb}^{eff}(r)=\frac{}{𝐫}V_{bb}^{eff}(r).$$
(8)
This average force provides an intuitive way of understanding the effective interactions which parallels the more formal derivation above. First we define a one-body density as follows:

$$\rho _s^{(1)}(𝐫^{})=<<\underset{i}{\overset{N_s}{}}\delta (𝐫^{}𝐫_𝐢^𝐬)>>,$$
(9)
where $`<<>>`$ denotes a canonical thermodynamic average over the small particles. If one big particle is fixed at the origin, and one is fixed at r, then we can define a one-particle density of the small particles $`\rho _s^{(1)}(𝐫^{};\mathrm{𝟎},𝐫)`$ for that fixed configuration. The effective force induced between two big spheres by the smaller particles can then be written as:

$$𝐅_{bb}^{eff}(𝐫)=\rho _s^{(1)}(𝐫^{};\mathrm{𝟎},𝐫)\frac{}{𝐫^{}}\varphi _{bs}(r^{})𝑑𝐫^{},$$
(10)
which has the obvious physical interpretation that the total effective or mean force is simply the average of the sum over all the big-small interactions.

In the next sections, we will concentrate on the effect of changing $`\mathrm{\Phi }_{ss}(r)`$ and $`\mathrm{\Phi }_{bs}(r)`$, at fixed reservoir density $`\rho _s`$, on the effective interactions $`F_{bb}^{eff}(r)`$ and $`V_{bb}^{eff}(r)`$.

## III Models

### A Interactions

We model the basic small-small and big-small interactions defined in Eq. (2) as hard-core Yukawa pair potentials:

$$\mathrm{\Phi }_{ij}(r)=\{\begin{array}{cc}\mathrm{}\hfill & \hfill r<\sigma _{ij}\\ \varphi _{ij}(r)\hfill & \hfill \text{otherwise},\end{array}$$
(11)
where in each case $`r`$ denotes the distance between the centers of the relevant particles. Throughout, we take $`\varphi _{bb}(r)=0`$, while the big-small interaction added to the HS repulsion is given by:

$$\varphi _{bs}(r)=\frac{ϵ_{bs}\sigma _{bs}}{r}\mathrm{exp}\left[\kappa _{bs}(r\sigma _{bs})\right],$$
(12)
and the small-small interaction added to the HS repulsion is given by:

$$\varphi _{ss}(r)=\frac{ϵ_{ss}\sigma _{ss}}{r}\mathrm{exp}\left[\kappa _{ss}(r\sigma _{ss})\right],$$
(13)
where $`\sigma _{bs}=\frac{1}{2}(\sigma _{bb}+\sigma _{ss})`$.

The Yukawa tail can be either repulsive or attractive. This gives us the possibility to change the interaction over a broad range from pure HS (which serve as a reference case) to soft repulsions and attractive tails. We note that a hard-core Yukawa interaction has been frequently used in theoretical studies on liquids, see e.g.  and refs therein. Both an attractive and repulsive tail is realized in various colloidal solutions: a short-ranged attractive Yukawa tail has been shown to satisfactorily model the “stickiness” in the interaction between globular protein solutions . A repulsive Yukawa tail, on the other hand, describes charged suspensions  where the range is controlled by the added salt concentration and the amplitude by the colloidal charge.

### B Parameter combinations

Our aim is to systematically investigate the effect of the basic big-small and small-small interactions on the effective big-big interactions. We vary the big-small interaction by changing $`\beta ϵ_{bs}`$, and vary the small-small interaction by changing $`\beta ϵ_{ss}`$. The other parameters are fixed as follows: $`\sigma _{ss}=0.2\sigma _{bb}`$, $`\kappa _{bs}=6/\sigma _{bb}=1.2/\sigma _{ss}`$, and $`\kappa _{ss}=15/\sigma _{bb}=3/\sigma _{ss}`$. The range of the big-small interaction is of the order of the small particle size, while the range of the small-small interactions is significantly less than the small particle size. For all simulations the packing fraction of the small particles in the bulk is set to $`\eta _s=\pi \rho _s\sigma _{ss}^3/6=0.1`$.

In total, we studied $`9`$ different parameter combinations corresponding to taking $`\beta ϵ_{bs}`$ and $`\beta ϵ_{ss}`$ to be positive, negative, or zero. The detailed values are given in table I.

## IV Computer simulations of effective forces and interactions

### A Simulation method

Our simulation set-up contains two big spherical particles in a large cubic simulation box of length $`L`$ with periodic boundary conditions in all three directions. The big particles with their center-to-center separation $`r`$ are fixed along the body diagonal of the cubic box and $`N_s=10000`$ small mobile particles of diameter $`\sigma _{ss}`$ are added to the box. The box length $`L=7.49\sigma _{ss}`$ is sufficiently large to exclude any spurious periodic image effects. We also studied a single big particle in the simulation box to access the one-body small-particle density $`\rho _s(r)`$ around an isolated big sphere. For two big spheres $`\rho _s=N_s/(L^3\frac{1}{3}\pi \sigma _{bb}^3)`$, while for one big sphere $`\rho _s=N_s/(L^3\frac{1}{6}\pi \sigma _{bb}^3)`$. We checked the asymptotic density profiles to confirm that the bulk densities were always the same.

We use a molecular dynamics (MD) code combining the velocity Verlet algorithm with discrete collisions and reflections induced by the hard core of the interaction potentials in order to calculate the trajectories of the small particles. The system was carefully equilibrated and then statistical averages were computed such as the mean force acting on the big particles (see Eq. (10)) or the inhomogeneous density field $`\rho _s(r)`$ of the small particles around a single big particle. Details of our simulation procedure are described in references.

### B One body density profiles

Before we discuss the depletion potentials it is useful to first examine the density profiles $`\rho _s(r)`$ of the small particles around a single sphere.

#### 1 No extra big-small interaction: runs 1-3

In Fig. 2, these are shown for runs $`13`$, where $`ϵ_{bs}=0`$. Compared to the density profile of the pure HS, there is a significantly increased accumulation near the big spheres when a small-small repulsive interaction is added. To first order this can be understood by mapping the $`\varphi _{ss}(r)`$ onto an effective HS diameter $`\sigma _{ss}^{eff}`$: Adding a small-small repulsive interaction increases the effective HS size, and therefore also the effective HS packing fraction, resulting in a more pronounced accumulation of density near the surface. Adding a small-small attraction has the opposite effect. This is because the bulk is now more favorable for the small particles, and they are attracted to it.

#### 2 Added big-small repulsion: runs 4-6

In Fig. 3, the density profiles are shown for runs $`46`$, where $`\beta ϵ_{bs}=0.82`$. As expected, the repulsive big-small interaction leads to a reduction of the density at contact. Again, adding a small-small repulsion increases the density at contact, and adding a small-small attraction leads to a further depletion of the density at contact.

#### 3 Added big-small attraction: runs 7-9

In Fig. 4, the density profiles are compared for runs $`79`$, where $`\beta ϵ_{bs}=0.82`$. As expected, the attractive big-small interaction results in an enhanced density at contact. Similarly to the previous two cases where $`ϵ_{bs}=0`$ and $`ϵ_{bs}>0`$, respectively, adding a small-small repulsion enhances the contact density w.r.t. the case of no small-small repulsion. But in contrast to the two previous cases, where adding a small-small attraction resulted in a depleted density profile w.r.t. the pure HS case, here adding $`ϵ_{ss}<0`$ results in an enhanced total accumulation of the small particles near the big one. Although the contact value is slightly smaller than the case for pure HS, there is a marked accumulation further out, corresponding to a second layer of particles. The relative adsorption of the small particles around the big one is thus considerably larger than for the case of no small-small attractions. This can be understood by the following simple argument: the big-small attraction leads to an accumulation of the small particles near the surface of the big sphere. When this accumulation is large enough, it becomes favorable for the small particles to leave the bulk, and approach the surface of the big particle, where their local density is larger. Thus the two attractions amplify each other in a non-linear fashion.

In conclusion then, when comparing Figs. 2, 3 and 4, it is clear that for the parameters we have chosen, the big-small interaction has the largest relative effect on the density profiles. The effect of adding a small-small repulsion can be qualitatively understood by the larger effective sphere size and concomitant larger packing fraction. The effect of adding a small-small attraction can be qualitatively understood by the fact that the bulk is usually preferred over the surface of the particle, except when the big-small interaction is strong enough to provoke a non-linear enhancement of the density of small particles near a big one.

Having investigated the effect of the interactions on the one-body density profiles, we now turn to the related two-body depletion potentials.

### C Effective pair forces and potentials

In Fig. 5 we compare the effective depletion force between two big spheres for all nine parameter combinations detailed in table I. In Fig. 6 we compare the related depletion potentials. First we note that changing the big-small interaction is the dominant effect: the depletion forces and potentials split naturally into three groups: no added big-small interaction (solid lines, runs $`13`$), big-small attraction (dotted lines, runs $`46`$), and big-small repulsion (dot-dashed lines, runs $`79`$). We will treat each case in turn.

#### 1 No extra big-small interaction: runs 1-3

As was already seen for the one-body profiles, adding a small-small repulsion (run 2) results in a larger effective small-sphere size and packing fraction, which is reflected in more pronounced oscillations compared to the pure HS case (run 1). These are evident both in the effective force and in the effective pair potential. Adding a small-small attraction (run 3) results in a reduced density near a big particle, as seen in Fig. 2. This would imply that each big sphere excludes slightly more free volume than if there were no small-small attractions, which, in turn, implies a slightly more attractive potential, as is observed.

#### 2 Added big-small repulsion: runs 4-6

The dominant effect of adding a big-small repulsion is to make the effective forces and potentials much more attractive. This can be understood with the classical picture of depletion: Adding a big-small repulsion results in a larger depletion layer around each large particle. When two large spheres approach, the amount of doubly excluded volume is therefore larger, resulting in a more attractive effective (depletion) potential or force. Again, adding a small-small repulsion (run 5) results in enhanced layering as compared to the pure HS small particles. Adding a small-small attraction (run 6) has only a weak effect similar to what was seen for run 3.

#### 3 Added big-small attraction: runs 7-9

The dominant effect of adding a big-small attraction is to make the effective forces and potentials much more repulsive. This can be qualitatively understood from the fact that the big-small attraction results in an enhanced density of small particles near a single large sphere. When two such large spheres approach one another, the layers of small particles around each one begin to interact, leading to an enhanced repulsion between them. As demonstrated in Fig. 4 for a given big-small attraction, adding small-small repulsion or attraction both result in a further enhanced density $`\rho _s(r)`$ of the small particles around a single large one, which in turn explains why the effective pair potentials are more repulsive for both $`ϵ_{ss}>0`$ and $`ϵ_{ss}<0`$. Relatively speaking, adding small-small attractions has the largest effect on the effective potentials, which is due to the non-linear coupling between the big-small and the small-small attractions.

## V Theoretical descriptions of the effective forces and potentials

### A Superposition approximation

To calculate the effective forces via Eq. (10), one needs a prescription for calculating the one-body density of the small particles. In section IV, we essentially did this by computer simulations. In this section we approximate the full one-body density by a superposition of the one-body density $`\rho _s(r)`$ around an isolated single sphere:

$$\rho _s^{(1)}(𝐫^{};\mathrm{𝟎},𝐫)=\rho _s(r^{})\rho _s(|𝐫𝐫^{}|)/\rho _s,$$
(14)
an approach similar in spirit to the Kirkwood superposition approximation. Since all input information comes from the (radially symmetric) problem of a single sphere, this superposition approximation greatly simplifies the calculation of the two-body depletion forces. The input $`\rho _s(r)`$ could come from density functional theory or integral equation theory, as was done previously by other authors, but here we will use the $`\rho _s(r)`$ generated by our simulations and depicted in Figs. 24. The results are shown in Figs. 715, where we compare in detail the radial force $`F_{bb}^{eff}(r)=𝐅_{bb}^{eff}(r)𝐫/r`$ and the effective potential $`V_{bb}^{eff}(r)`$ calculated with the superposition approximation, to the results obtained by direct simulations.

Figure 7 (run 1) shows the HS reference case. Here the superposition approximation works quite well. The packing fraction of the small spheres $`\eta _s=0.1`$ is rather low, so we expect that the small-small correlation effects are not very strong. The total one-body density is therefore well approximated by Eq. (14). As demonstrated by other authors for the pure HS case, this superposition approximation begins to break down as $`\eta _s`$ increases and two-body correlation effects become more important. For example, they found that for large $`\eta _s`$, the superposition approximation leads to an overestimate of the strength of the attractive force at contact. This is exactly what is seen in Fig. 8 (run 2), where the effect of an added repulsive $`\varphi _{ss}(r)`$ can be understood in terms of a larger effective packing fraction $`\eta _s`$. Figure 9 (run 3) shows that for $`ϵ_{ss}<0`$ the superposition approximation shows a similar error to what was seen for run 2, i.e. the forces and potentials are too attractive.

In Figs. 1012 (runs 4-6), where the repulsive $`ϵ_{bs}(r)`$ induces much more attractive effective interactions than for the pure HS case, the superposition approximation is seen to work quite well for $`ϵ_{ss}=0`$ and $`ϵ_{ss}>0`$. This is most likely because the effective (depletion) force or potential is dominated by $`\varphi _{bs}(r)`$, which also induces a lower $`\rho _s(r)`$ (see Fig. 3), so that small-small correlations only play a relatively minor role and Eq. (14) is rather accurate. Even for $`ϵ_{ss}<0`$, the results are quite good, although the superposition approximation tends to predict forces and potentials that are too attractive, just as was found for $`\varphi _{bs}(r)=0`$.

In Figs. 1315 (runs 7-9), where the attractive $`\varphi _{bs}(r)`$ induces an increased local accumulation of the small spheres near each big sphere, the superposition approximation performs less well. This is perhaps not surprising, as the enhanced density of small particles near each big sphere results in a more important role for small-small correlations, which are not well treated by the superposition approximation.

In summary then, for HS or other repulsive small-small interactions, the superposition approximation works best for low (effective) packing fractions $`\eta _s`$, where correlations between the small particles do not strongly alter the full one-body density $`\rho _s^{(1)}(𝐫^{};\mathrm{𝟎},𝐫)`$ from the superposition of the densities around an isolated big sphere. The case of a strongly repulsive $`\varphi _{bs}(r)`$, which lowers the effective packing fraction near the big spheres, is particularly well described by the superposition approximation. In contrast, the case of an attractive $`\varphi _{bs}(r)`$, which results in an increased accumulation of small particles around each big sphere, is not as well described.

### B Density Functional Theory

In a recent development, density functional theory (DFT) has been used to derive effective potentials for additive and non-additive HS mixtures. For the additive case, quantitative agreement with computer simulations was achieved. Since the calculations for the non-additive case were shown to be equivalent to the additive ones, a similar quantitative accuracy could be expected there. In brief, the method is based on the following exact relationship between the effective potential and the one-body direct correlation function:

$$\beta V_{bb}^{eff}(𝐫)=\underset{\mu _b\mathrm{}}{lim}(c_b^{(1)}(\mathrm{})c_b^{(1)}(𝐫))$$
(15)
where $`c_b^{(1)}(𝐫)`$ is defined for the case where a big sphere is fixed at the origin, and exerts a field on the small spheres and on a big test particle inserted at $`𝐫`$. DFT can therefore provide a route to the effective potentials since

$$c_b^{(1)}(𝐫)=\beta \frac{\delta _{ex}[\rho _b,\rho _s]}{\delta \rho _b(𝐫)}$$
(16)
where $`_{ex}[\rho _b,\rho _s]`$ is the excess (over ideal) intrinsic Helmholtz free energy functional of the mixture. Given some mixture functional, one can obtain the effective potentials and forces from a radially symmetric calculation of $`c_b^{(1)}(𝐫)`$ in the limit of vanishing density of the big spheres. This is much simpler than say trying to use DFT to directly calculate the cylindrically symmetric one-body density $`\rho _s^{(1)}(𝐫^{};\mathrm{𝟎},𝐫)`$, for use in Eq. (10). All one now needs is some prescription for the mixture functional. Here we use Rosenfeld’s very successful fundamental measure theory, in its form valid for HS mixtures. As found previously, the DFT and simulation results for effective interactions of the pure HS case depicted in Fig. 7 (run 1) agree very well.

This DFT method can easily be extended to an arbitrary $`\varphi _{bs}(r)`$, since this simply corresponds to an additional external field in Eq. (15). Previous DFT calculations found good agreement with earlier simulations with an attractive $`\varphi _{bs}(r)`$. Here we also find very good agreement with the runs $`4`$, and $`7`$, (Figs. 10 and 13) which correspond to a finite $`\varphi _{bs}(r)`$ but no additional $`\varphi _{ss}(r)`$.

Since there is at present no successful two-component DFT for mixtures where the small-particles are not HS like, some approximations must be made. We chose to map the $`\varphi _{ss}(r)`$ onto effective HS diameters $`\sigma _{ij}^{eff}`$ using the Barker-Henderson approach:

$$\sigma _{ij}^{eff}=\sigma _{ij}+_{\sigma _{ij}}^{\mathrm{}}(1\mathrm{exp}[\beta \varphi _{ij}(r)])𝑑r,$$
(17)
where $`\sigma _{ij}`$ is the bare HS diameter. We then calculated effective pair potentials and effective pair forces with our full DFT approach, including an explicit $`\varphi _{bs}(r)`$ when needed, but with the small-spheres mapped onto the effective diameters $`\sigma _{ss}^{eff}=1.447\sigma _{ss}`$ and $`\sigma _{ss}^{eff}=0.658\sigma _{ss}`$ for the repulsive and attractive $`\varphi _{ss}(r)`$ respectively. Similarly the effective packing fraction becomes $`\eta _{ss}^{eff}=0.303`$ for the repulsive and $`\eta _{ss}^{eff}=0.0290`$ for the attractive interaction.

Firstly, for runs $`2`$, $`5`$, and $`8`$ (Figs. 8, 11, and 14, respectively), which all correspond to a repulsive $`\varphi _{ss}(r)`$, we find good agreement for no additional $`\varphi _{bs}(r)`$ (see Fig. 8), but less good agreement for a repulsive or attractive $`\varphi _{bs}(r)`$ (see Fig. 11 and Fig. 14 respectively).

Because of the strength of the small-small attraction, the mapping results in a very low effective packing fraction of the small particles. Overall, the DFT underestimates the effective forces for $`ϵ_{bs}=0`$ and $`ϵ_{bs}>0`$, as can be seen in Figs. 9 and 12. It performs rather poorly for run 9, where the attractive $`\varphi _{bs}(r)`$ results in a non-linear enhancement of the small-particle density profile, an effect not taken into account with our HS mapping. This suggests that a different two-component DFT, which explicitly takes into account the small-small attraction $`\varphi _{ss}(r)`$ needs to be developed before we can use this route to derive accurate effective pair forces and potentials.

### C Mapping to non-additive HS system

In a previous paper two of us proposed that the effects of big-small and small-small interactions on the $`V_{bb}^{eff}(r)`$ could be understood by mapping onto those of non-additive HS systems. These have the advantage that they can be determined by an exact mapping onto the depletion potentials of additive HS mixtures, which, in turn, are well understood and for which a good parameterization exists. Even a small non-additivity was shown to have a large effect on the depletion potentials.

By mapping the big-small and small-small interactions onto effective HS diameters through Eq. (17) one can define the non-additivity in terms of the parameter $`\mathrm{\Delta }`$:

$$\sigma _{bs}^{eff}=\frac{1}{2}\left(\sigma _{bb}+\sigma _{ss}^{eff}\right)\left(1+\mathrm{\Delta }\right).$$
(18)
Four different ways of adding interactions to introduce non-additivity were studied:

(i) repulsive $`\beta \varphi _{ss}(r)`$: $`\sigma _{bs}^{eff}=\sigma _{bs}`$; $`\sigma _{ss}^{eff}>\sigma _{ss}`$; $`\mathrm{\Delta }<0`$

(ii) attractive $`\beta \varphi _{ss}(r)`$: $`\sigma _{bs}^{eff}=\sigma _{bs}`$; $`\sigma _{ss}^{eff}<\sigma _{ss}`$; $`\mathrm{\Delta }>0`$

(iii) repulsive $`\beta \varphi _{bs}(r)`$: $`\sigma _{bs}^{eff}>\sigma _{bs}`$; $`\sigma _{ss}^{eff}=\sigma _{ss}`$; $`\mathrm{\Delta }>0`$

(iv) attractive $`\beta \varphi _{bs}(r)`$: $`\sigma _{bs}^{eff}<\sigma _{bs}`$; $`\sigma _{ss}^{eff}=\sigma _{ss}`$; $`\mathrm{\Delta }<0`$.

For details of how each of the different cases of non-additivity affect the depletion potentials, we refer to refs..

An example where the mapping to non-additivity works well is given in Fig. 8 (run 2), where we mapped the small-small repulsion onto an effective HS diameter (case (i)), so that the DFT result is really that of a non-additive HS mixture with $`\mathrm{\Delta }=0.069`$. In ref.  we also found semi-quantitative agreement with the mapping for weak and short-ranged $`\beta \varphi _{bs}(r)`$. In the present case, where $`\beta \varphi _{bs}(r)`$ is stronger and longer ranged, the agreement is no longer quantitative. Nevertheless, for runs 1-8 the mapping scheme provides a qualitative explanation of the trends. I.e. for a fixed $`\varphi _{bs}(r)`$, adding a repulsive (case (i)) or attractive (case (ii)) $`\varphi _{ss}(r)`$ causes $`V_{bb}(r)`$ to become more repulsive or attractive, respectively. Similarly, for a $`\varphi _{ss}(r)`$, adding a repulsive (case (iii)) or attractive (case (iv)) $`\varphi _{bs}(r)`$ results in a more attractive or repulsive $`V_{bb}(r)`$, respectively. The only case where this scheme breaks down is run 9, where adding $`\varphi _{ss}(r)<0`$ for an attractive $`\varphi _{bs}(r)`$ should fall under case (ii). But instead of inducing more attraction, the effective big-big potential becomes more repulsive. Of course this is not surprising, since the DFT results already showed that a mapping scheme misses the non-linear enhancement of the small particle density profiles.

In conclusion then, the mapping to non-additivity works best for a repulsive small-small interaction. For repulsive and attractive big-small interactions, the mapping is only qualitative. The much better semi-quantitative agreement found in ref.  can be traced to the much weaker effect of the $`\varphi _{bs}(r)`$ used there. Just as was seen for the direct DFT methods, it is the case of small-small attraction combined with big-small attraction which seems most difficult to capture within our mapping scheme.

## VI conclusions

In conclusion, we have shown how the basic interactions in a colloidal mixture control the resulting effective interactions between the big particles. This knowledge may be exploited to stabilize colloidal particles against coagulation and to tailor the colloidal phase diagrams.

Adding a repulsive $`\varphi _{bs}(r)`$ results in a strongly enhanced attraction through the standard depletion mechanism. We also found at least two ways to obtain significantly more repulsive effective interactions caused by accumulation of the small particles (accumulation repulsion): (1) Adding repulsions between the small particles results in an enhanced accumulation near the surface of the large particles; when two large particles approach each other, this results in an effective repulsion between them. (2) Adding an attraction between the large and the small particles also results in an enhanced accumulation near the surface of the large particles and therefore in repulsive effective interactions. Furthermore, we found that for an attractive big-small interaction, adding small-small attractions resulted in even more effective repulsion. This “repulsion through attraction” effect is caused by a coupling between $`\varphi _{ss}(r)`$ and $`\varphi _{bs}(r)`$.

These predictions could in principle be verified in experiments that directly measure the effective forces of colloidal suspensions, such as optical tweezers or Total Internal Reflection Microscopy (TIRM). They could also be verified indirectly through measurements of phase-behavior and coagulation. Measurements of the second osmotic virial coefficient might also be very sensitive probes of the effective interactions. Systems where these interactions could be tuned include for example ternary suspensions where the small-small attraction is generated by depletion attraction of an even smaller third colloidal component or highly salted charged suspensions where the van-der-Waals attraction dominates.

In some interesting recent experiments, a colloidal suspension of neutral big particles was stabilized by the addition of highly charged small nanoparticles. The proposed mechanism was termed “nanoparticle halos”, and is very similar to the mechanism we observe, for example, in our run $`2`$, where a repulsive $`\varphi _{ss}(r)`$ of the screened Coulomb (Yukawa) form was added to the smaller particles, resulting in an increased accumulation of small particles near each big one, and an effective repulsion between the big particles. The effects in the experiments may also be enhanced by a small attractive $`\varphi _{bs}(r)`$ (like our run $`8`$). We are currently actively pursuing a more detailed comparison with these experiments.

Another possible application of this work is to supercritical solvents, which have important applications in industrial processes. The question of how the effective force on the big particles depends on the interactions with a low density solvent is encountered there as well, sometimes for similar size ratios.

More generally to make useful predictions for practical applications, one needs to more thoroughly explore the rather large parameter space, which includes the $`ϵ_{ij}`$, the $`\kappa _{ij}`$, $`\eta _s`$ and the ratio $`\sigma _{ss}/\sigma _{bb}`$. This would be very tedious with simulations – trustworthy theoretical techniques would be more practical. We attempted a number of theoretical descriptions of the effective interactions. The superposition approximation, which has the advantage of only needing one-body input, works best for a strong big-small repulsion, but becomes much less accurate when the (local) small particle packing is high. When the smaller component is purely HS like, then our DFT approach is very accurate, just as was found for additive and non-additive HS mixtures. When a repulsive interaction is added between the small particles, DFT does not perform quite as well, and when an attraction is added between the small particles, the differences are even more important. By mapping the small-small and big-small interactions onto effective HS diameters, we can map onto an effective non-additive HS model. This explains the qualitative trends for most of our parameter combinations. It is quantitative if only a small-small repulsion is added, but breaks down when both a small-small and a big-small interaction couple together to induce an enhanced small particle density.

All three tested theoretical methods perform less well for added small-small attractions. Constructing a reliable theory to treat this very interesting case is therefore a challenging problem. One possible way to extend the DFT calculations to an attractive small-small interaction would be to add a mean-field attractive term. For example, mean-field like functionals based on thermodynamic perturbation theory around a HS reference system have been successfully applied to the density profiles and phase-behavior in systems with attractive potentials .

New physics is expected when wetting or drying phenomena control the density of the small particles between the big ones. This is relevant if the bulk fluid of the small particles is close to liquid-gas phase coexistence. A wetting transition is expected to have a profound impact on the effective interaction as well, one important effect is a liquid “bridge” of small particles between the big one which has been recently studied in more detail .

We thank Jean-Pierre Hansen, Geoff Maitland, Edo Boek, and Bob Evans for illuminating discussions. AAL acknowledges support from the Isaac Newton Trust, Cambridge, HL acknowledges support from Schlumberger Cambridge Research, through the Schlumberger visiting professor scheme.

i


format: text

text:
# New mechanism of membrane fusion

## I Introduction

Fusion of membranes is involved in basic biological processes but its mechanism remains poorly understood. That different proteins trigger different fusion events obscures the possibility of a stage common to them, one which depends only on the properties of membrane bilayers themselves. The existence of such a universality is made plausible by the realization that, in all fusion processes, membranes must merge, and the properties of the membranes will most likely determine the nature of this stage, the fusion intermediate. Time and length scales of fusion events are of the order of microseconds  and nanometers, so that direct measurements of intermediate structure have not been possible. Limited theoretical treatment of this problem, based on membrane elasticity theory, has focused on the so-called stalk mechanism. In this scenario, the leaves of the two cis membranes, closest to one another, fuse forming a stalk. The stalk expands radially and thins so that the trans leaves make contact and form a single bilayer. This stage is denoted hemifusion. Hole formation in this bilayer completes the fusion pore. Due to the lack of direct experimental confirmation, this mechanism, although plausible, remains hypothetical.

To obtain a direct view of the fusion mechanism we have undertaken the first Monte Carlo simulation of membrane fusion, one in which membranes are treated on the molecular level. In contrast to the stalk hypothesis, we find a very different fusion mechanism. Although stalk formation does take place in the initial stage of fusion, it does not result directly in the fusion pore. Instead, the stalk destabilizes the contacting bilayer membranes and promotes formation of transient holes next to it. When two holes are nucleated in apposing bilayers in the vicinity of the stalk, the latter grows along the hole edges, sealing their periphery like a zipper, and eventually closes thereby forming the fusion pore.

## II The model

Our model describes membranes formed by single-chain amphiphiles, like block copolymers which exhibit the same phases as do biological lipids, and also form vesicles. It has the advantage that it has been well studied, permits detailed analysis of molecular configurations, and is well suited to processes occurring on the small time and length scales characteristic of fusion. The amphiphiles are treated using the bond fluctuation model  in which each molecular segment occupies a cube of a three-dimensional lattice. The eight lattice sites defining the cube cannot be occupied by another segment centered on neighboring sites. Segments along an amphiphile are connected by one of 108 bond vectors of lengths $`2,\sqrt{3},\sqrt{5},3`$ or $`\sqrt{10}`$, measured in units of the lattice spacing $`a_0`$. Mapping this model onto lipids in solution, we find the lattice spacing to correspond to approximately 1Å. The large number of bond vectors and the extended segment shape allow a rather faithful approximation of continuous space, while retaining the computational advantages of lattice models. The amphiphilic molecules consist of $`N=32`$ segments, of which 10 are hydrophilic and 22 are hydrophobic. This particular choice of the ratio of hydrophilic and hydrophobic segments results in the diblock system being close to coexistence of the lamellar and inverted hexagonal phases. The solvent is represented by a homopolymer, chains consisting of 32 hydrophilic segments. Like segments attract each other and unlike segments repel each other via a square well potential which comprises the nearest 54 lattice sites. Each contact changes the energy by an amount $`ϵ=0.177k_BT`$. The particular choice of the interaction parameter $`ϵ`$ guarantees that the interfacial width between hydrophilic and hydrophobic segments is not too small to be comparable to the lattice spacing, and at the same time results in well-defined bilayers.

The simulation cell is $`L\times L`$ in the $`x,y`$ directions and of length $`D`$ in the $`z`$ direction, with $`L=156a_0`$ and $`D=96a_0`$. Periodic boundary conditions are utilized in all three directions. The monomer density of the system is $`\rho =1/(16a_0^3)`$, corresponding to 146,016 segments within the volume, or $`2376`$ amphiphiles and $`2187`$ homopolymers. To encourage the fusion of the bilayers that we study, we prepare them under tension by providing only enough amphiphiles to form bilayers of thickness $`25a_0`$, thinner than their equilibrium thickness of $`25.2a_0`$, which we determined independently. Given that the system is nearly incompressible, this corresponds to a fractional increase of area of less than 0.8%, one easily sustained by biological membranes. Two such bilayers are created parallel to the $`xy`$ plane, stacked one upon the other, thereby mimicking the dehydration which permits close bilayer contact, a circumstance known to promote fusion. Sixty-four independent starting configurations were prepared. Monte Carlo simulations were performed in the canonical ensemble. The conformations are updated by local monomer displacements and slithering–snake like movements. The different moves are applied with a ratio $`1:3`$. The latter moves do not mimic the realistic dynamics of lipid molecules and we cannot identify the number of Monte Carlo steps with time. However, the density of hydrophilic and hydrophobic segments is conserved and the molecules diffuse. Consequently, we expect the time sequence on length scales larger than a single molecule to resemble qualitatively those of a simulation with more realistic dynamics.

To make sure that an isolated membrane is stable, we have simulated a single bilayer under the above conditions. The system was extremely stable with respect to hole formation, apparently due to a very high energy cost to form the hole’s edge, and exhibited the usual capillary wave fluctuations. In contrast, the system with two apposing membranes resulted in a range of structural transformations that eventually led to fusion pore formation.

## III Results

### A Qualitative observations

Initially the apposed bilayers formed a number of local connections, the so-called stalks. Formation of these stalks was promoted by the system being near bulk coexistence between lamellar and inverted hexagonal phases. In the stalk mechanism, these stalks grow radially and eventually form transmembrane contacts. In contrast, we observed anisotropic growth of the stalks, which resulted in the formation of structures whose length was a few times that of their width. These may be precursors of line defects observed during the lamellar to inverted-hexagonal transition.

Formation of the elongated stalks led to local destabilization of the bilayers. Holes formed next to these defects. We stress again that single, isolated, membranes were very stable. Thus hole formation in the apposed bilayers is apparently due to inter-bilayer interactions expressed in the stalk formation. It is easily seen that a hole formed next to a stalk has a smaller edge energy, in comparison to an isolated hole due to the reduction in curvature energy. This reduction is roughly proportional to the length of the edge adjacent to the stalk.

Hole formation was a transient process, i.e. numerous holes were observed to open and close reversibly. Only when two holes, one in each bilayer, nucleated next to, and on the same side of, a given stalk did the final stage of the fusion process commence. The stalk then propagated along the edges of the two holes bringing them together, like a zipper, to complete the fusion pore.

We show what one such partially formed pore looks like in FIG. 1. Panel (b) presents the top view of the two membranes. Only tail segments are shown. The pore through the two bilayers is clearly visible. Panel (c) presents the top view of the layer between the two membranes. Hydrophobic segments are light grey, while hydrophilic segments are dark grey. Connectivity is established along a portion of the pore’s rim by the stalk, while the two hydrophobic sheets are still unconnected elsewhere. This is confirmed by a side view through the pore (d). A schematic of this intermediate is shown in panel (a).

### B Quantitative analysis

After 62,500 attempted local monomer displacements per segment and 187,500 slithering–snake movements per molecule, we find that the initially sharp interface of the bilayers have widened, and profiles across the membranes have adopted locally their stationary form. Profiles at a later stage are broadened due to fluctuations of the local membrane position, and to the formation of holes in the membranes. In order to reduce the effect of membrane fluctuations on our observation of the local profiles, we define at each $`x`$ and $`y`$ a local midpoint, $`z_{\mathrm{middle}}(x,y)`$, of the two-bilayer complex as follows. For each $`x`$, $`y`$ and $`z`$, we determine the total number of amphiphile segments, $`n_a(x,y,z)`$, in a volume centered at these coordinates, of width $`B=28a_0`$ in the lateral directions, and $`50a_0`$ vertically, the thickness of two bilayers. The value of the coordinate $`z`$ which maximizes $`n_a(x,y,z)`$ defines $`z_{\mathrm{middle}}(x,y)`$. We have chosen a rather large lateral length scale, $`B`$, so that a small hole in a single bilayer does not significantly alter our estimate for the midpoint of the two bilayers.

To analyze structural changes in the system, we looked not only at the fusion pores traversing both bilayers, but also at the holes in each bilayer which we found to occur about an order of magnitude more frequently than pores. We did this by calculating the local density of hydrophobic tail segments coarse grained over $`4a_0\times 4a_0`$ square columns centered at $`x`$ and $`y`$. Holes and pores should have a very low, or vanishing, density of tails compared to the areal density of two normal bilayers. If the local value of the coarse-grained density of tail segments within a column extending through the whole $`z`$-range falls below 25% of the normal areal density of the two bilayers, we define this column as belonging to a pore traversing them. To locate a hole in only one of the bilayers, the integration over the $`z`$ coordinate, rather than extending over the whole simulation box width, $`D`$, is limited to either the upper or lower half of it starting at the point between the two bilayers. If the local coarse-grained tail density falls below 25% of the normal areal density of a single bilayer, the column belongs to a hole. Pores and holes are defined as the aggregate of nearest neighbor $`4a_0\times 4a_0`$ plaquettes with low hydrophobic segment density.

The probability of hole sizes in the two apposed bilayers and the isolated bilayer are presented in FIG. 2. Not only is the absolute frequency of hole formation very much less in the isolated membrane than in the two apposed bilayers, but its decay with increasing hole size is also much more rapid. This indicates that the edge free energy of holes in the isolated bilayer is larger and corroborates the observation that stalk formation promotes hole formation when two bilayers are in contact.

To determine unambiguously whether fusion pore formation involves a stage in which the trans leaves make contact, one unavoidable in the standard stalk mechanism, we perform a statistical analysis of the fusion intermediates observed in our simulations. We examine in each simulation run the first pore which exceeds a certain radius, $`R_{min}=5.5a_0`$, one sufficiently large that the pore must be open; that is, there is a channel connecting upper and lower parts of the system accessible to solvent molecules.

Our strategy is to look at the region along the rim, or edge, of the pore in the pore’s mid plane. If fusion proceeds via hemifusion this region has to be essentially completely hydrophobic. In contrast, if fusion proceeds through the formation of two transient holes in apposing bilayers, then we may see pores in an intermediate stage in which the holes are aligned, but the rim of the pore has not fully formed. In such a case, the region along the rim of the pore would be at least partially hydrophilic. We measure the density of segments in a series of columns of height $`5a_0`$ and base $`4a_0\times 4a_0`$ centered along the rim of the pore, which is defined as the points in the $`xy`$-mid plane with a distance $`14a_0\pm 2a_0`$ from the center of the pore. This distance corresponds to a sum of radius of the pore, $`5.5a_0`$, and half of the hydrophobic thickness of the bilayer, $`8.5a_0`$.

To quantify the connectivity of the pore’s rim, we extract the distribution of the density difference $`\mathrm{\Delta }\varphi `$ between hydrophilic and hydrophobic segments in this region. On average there are $`\varphi _0=5`$ monomers in a column whose volume is $`4a_0\times 4a_0\times 5a_0`$. If the connectivity of the hydrophobic regions of the two bilayers were complete, as in the stalk mechanism, then we expect the rim region to be hydrophobic, so that the distribution of $`\mathrm{\Delta }\varphi `$ from all these columns would exhibit a single peak around $`\mathrm{\Delta }\varphi /\varphi _0=1`$. This is not at all what we see. Instead, only 12 pores have this single peaked distribution, whereas the rest are either bimodal (37 configurations) with a second peak centered about $`\mathrm{\Delta }\varphi /\varphi _0=+1`$ or have a peak at $`\mathrm{\Delta }\varphi /\varphi _0=1`$ with a long “tail” (15 configurations) toward positive $`\mathrm{\Delta }\varphi `$. From this we can conclude that in the majority of cases there is a considerable portion of the rim which is still hydrophilic; that is, the rims of the two holes are only partially connected. The average of all distributions which we obtain are shown by the solid circles in FIG. 3. The distribution of the values of $`\mathrm{\Delta }\varphi `$ obtained from the series of columns along the rim of the particular pore in FIG. 1 are shown as a solid line.

## IV Discussion

The fusion mechanism we propose is similar to that considered for electrofusion, with the addition that we specify the form of the fusion intermediate. It is in accord with the recent observation in molecular dynamics simulations that lipid tails take on a much greater range of configurations, including those sampling the hydrophilic region, when membranes are brought into close contact. Our mechanism is similar in spirit to the initial stages of the mechanism of hemagglutinin-mediated viral fusion recently proposed by Bonnafous and Stegmann who suggest that before fusion can occur, a hole must first form in the target membrane. However they go on to suggest that this hole leads to fusion of the cis layers producing a hemifusion diaphragm in which a hole is nucleated. In contrast, hemifusion plays no role in the scenario we propose. We do observe strong mixing of amphiphiles in the cis leaves before pore formation. From such mixing, it is often inferred that hemifusion precedes complete fusion. Instead, we find that this exchange results simply due to formation of the stalks.

The destabilization of the membranes by the stalks that leads to formation of the transient holes in both bilayers, observed by us, is consistent with the work of Cevc and Richardsen who emphasize the fact that membrane fusion is strongly promoted by defects in the bilayer structure. It seems clear that the additional destabilization of bilayers by the defects embedded in the membranes prior to fusion facilitates formation of the transient holes necessary for fusion.

The destabilizing role of the stalk intermediates, observed in our simulations, is missing from the current phenomenological models of membrane fusion. Although we do not observe the standard stalk mechanism, our simulations do not rule it out under different thermodynamic conditions. Extensions of phenomenological approaches that include the new observed intermediates should clarify this point. We expect the fusion intermediate we have seen in our simulation to be similar to that in biological membranes. Its elucidation should facilitate the control and modification of the fusion process itself.

Note: While preparing this manuscript, we learned that Noguchi & Takasu, studying a very different system consisting of rigid amphiphilic molecules of three atoms and containing no solvent, observed fusion behavior over a limited temperature range which is similar to that we observed in our system. This may indicate a universality to the mechanism we have described.

###### Acknowledgements.

Financial support was provided by the National Science Foundation under grant No. DMR 9876864 and the DFG Bi 314/17 in the priority program “wetting and structure formation at interfaces”. Computer time at the NIC Jülich, the HLR Stuttgart and the computing center in Mainz are also gratefully acknowledged.


text:
# Multifractal current distribution in random diode networks

## Abstract

Recently it has been shown analytically that electric currents in a random diode network are distributed in a multifractal manner \[O. Stenull and H. K. Janssen, Europhys. Lett. 55, 691 (2001)\]. In the present work we investigate the multifractal properties of a random diode network at the critical point by numerical simulations. We analyze the currents running on a directed percolation cluster and confirm the field-theoretic predictions for the scaling behavior of moments of the current distribution. It is pointed out that a random diode network is a particularly good candidate for a possible experimental realization of directed percolation.

Directed percolation (DP) is a simple model for directed connectivity in a random medium MarroDickman ; Hinrichsen00a . It differs from ordinary isotropic percolation by the additional constraint that activity is restricted to percolate along a distinguished direction in space. Interpreting this direction as a temporal degree of freedom, directed percolation may also be viewed as a reaction-diffusion process $`A2A`$, $`A\mathrm{}`$. In this process particle creation and removal compete one another, leading to a non-equilibrium phase transition form an active (percolating) into an absorbing (non-percolating) phase. Because of its robustness, DP plays the role of a standard universality class of nonequilibrium phase transitions which may be of a similar importance as the Ising model in equilibrium statistical mechanics.

One of the simplest realizations of DP is directed bond percolation on a (tilted) hypercubic lattice. In this model neighboring sites are connected by bonds which are conducting with probability $`p`$ and impermeable otherwise. In contrast to isotropic percolation, the bonds conduct only in one direction and may be realized by randomly distributed diodes, as shown in Fig. 1. Obviously, the effective conductivity of the system along the distinguished direction depends on the percolation probability $`p`$. If $`p`$ is small, all cluster of connected sites are finite so that the conductivity vanishes in the limit of large distances. On the other hand, if $`p`$ exceeds a certain critical threshold $`p_c`$, there is a finite probability to find a directed path between two points, leading to a finite resistance when averaged over many independent samples. At the critical threshold $`p=p_c`$, the system undergoes a continuous phase transition where clusters of conducting paths display a fractal structure. In this case the currents running on the cluster are distributed in a nontrivial way.

The DP transition can be described by field-theoretic renormalization group methods ReggeonFT ; Janssen81 . In this framework the critical behavior near the transition can be characterized by simple scaling laws which involve three independent critical exponents. The order parameter in the percolating phase is the density of sites belonging to an infinite cluster or – using the reaction-diffusion language – the density of active sites. Approaching the transition, this density scales as $`\rho (pp_c)^\beta `$, where $`\beta `$ is the critical exponent associated with the order parameter. Similarly, the correlation length along the distinguished direction scales as $`\xi _{}|pp_c|^\nu _{}`$, while the correlation length perpendicular to this direction diverges as $`\xi _{}|pp_c|^\nu _{}`$. The exponents $`\nu _{}`$ and $`\nu _{}`$ are generally different since introducing a distinguished direction in space breaks the symmetry between parallel and perpendicular degrees of freedom. In a two-loop approximation, these exponents are given by Janssen81 ; ExponentsFT

$$\begin{array}{cc}\hfill \beta & =1ϵ/60.01128ϵ^2+O(ϵ^3)\hfill \\ \hfill \nu _{}& =1+ϵ/12+0.02238ϵ^2+O(ϵ^3)\hfill \\ \hfill \nu _{}& =1/2+ϵ/16+0.02110ϵ^2+O(ϵ^3)\hfill \end{array}$$
(1)
where $`ϵ=4d`$. The exponents have also been estimated by various numerical methods (see Table 1). The triplet of exponents $`(\beta ,\nu _{},\nu _{})`$ labels the universality class of $`d`$+1 dimensional DP, while other critical exponents are usually related to these standard exponents by simple scaling relations. This type of scaling invariance, where only a finite number of independent exponents is involved, is referred to as simple scaling.

Let us now consider an electric current running on a directed percolation cluster according to Kirchhoff’s laws. By introducing such a current the theory is extended by an additional physical concept with non-trivial properties. In fact, even though the cluster itself is characterized by simple scaling laws, the currents on the cluster turn out to be distributed in multifractal manner, meaning that all moments scale with individual exponents which cannot be related to the standard exponents by simple scaling laws. In the case of isotropic percolation, this phenomenon was first conjectured by Rammal et al. rammal\_etal\_85  as well as by Arcangelis et al. arcangelis\_etal\_85 . In the meantime these results have been substantiated by field-theoretic renormalization group studies park\_harris\_lubensky\_87 ; ResistorNetwork  and numerical simulations rammal\_etal\_85 ; Batrouni ; Barthelemy . Recently, Stenull and Janssen StenullJanssen2  investigated the directed case. By employing field theoretic methods they determined the critical exponents of the moments of the current distribution. Their results clearly indicate multifractality. The purpose of the present work is to verify these predictions quantitatively by extensive numerical simulations. To this end we analyze a random diode network as shown in Fig. 1, which is a special case of the model considered in StenullJanssen2 . Again we would like to stress that multifractality is not observed as long as one is concerned exclusively with the underlying DP process. An additional process, here the transport of electrical currents, is necessary for the emergence of multifractality.

The moments of the current distribution in a $`d`$+1-dimensional DP process are defined as follows. Let us denote by $`x_{}`$ the coordinate along the distinguished direction (=time) and by $`x_{}`$ the vector of $`d`$ perpendicular coordinates (=space). An important quantity in the theory of DP is the pair-connectedness function $`C(x_{},x_{};x_{}^{},x_{}^{})`$, which is the probability to find a conducting path running from the point $`(x_{},x_{})`$ to the point $`(x_{}^{},x_{}^{})`$. The union of all possible paths in a given realization of conducting and non-conducting bonds is called the backbone of the correlation function. Note that the backbone is invariant under reversal of all diodes. In the field-theoretic framework this invariance is known as the so-called rapidity reversal symmetry RapidityReversal .

We now assume the diodes to be ideal in the sense that they are perfectly isolating if the applied voltage is negative while they have a fixed resistance $`R`$ if the voltage is positive. If we apply an external voltage at the two points $`\text{x}=(x_{},x_{})`$ and $`\text{x}^{}=(x_{}^{},x_{})`$ there will be a current running on the backbone which can be computed by solving Kirchhoff’s laws. Obviously the currents are distributed in a non-trivial way, especially when the backbone consists of many nested loops, as shown in Fig. 2. This current distribution has multifractal properties which can be studied by considering the moments

$$M_{\mathrm{}}=\underset{b}{}(I_b/I)^2\mathrm{},$$
(2)
where the sum runs over all bonds $`b`$ of the backbone while $`I_b`$ and $`I`$ denote the current passing bond $`b`$ and the total current, respectively. Each of these moments has a specific physical meaning. For example, $`M_0`$ is just the total number of bonds of the backbone while $`M_1`$ is the total resistance between the two points in units of $`R`$. $`M_2`$ is the second cumulant of the resistance fluctuations and may be considered as a measure of the noise in a given realization. Finally, $`M_{\mathrm{}}`$ is the number of so-called red bonds which carry the full current $`I`$. These bonds play a special role since they are shared by all possible paths.

At criticality the moments obey a power-law StenullJanssen2

$$M_{\mathrm{}}(x_{}^{}x_{})^{\psi _{\mathrm{}}/\nu _{}}.$$
(3)
Multifractal behavior manifests itself in a nonlinear dependence of the exponents $`\psi _{\mathrm{}}`$ on their index $`\mathrm{}`$. The exponents $`\psi _{\mathrm{}}`$, just as the geometrical exponents compiled in Eq. (1), are known in a two-loop approximation. Thus, at least in principle, we can gather the $`ϵ`$ expansion for the exponents $`\psi _{\mathrm{}}/\nu _{}`$ from the literature. This $`ϵ`$ expansion, however, would hardly be useful for comparison to our numerical data because it is valid only asymptotically in the vicinity of $`d_c=4`$. To produce reliable predictions for $`d=3`$ or even $`d=2`$ the $`ϵ`$ expansion has to be improved by incorporating additional information. Here, we carry out a so-called rational approximation. To be specific, we incorporate the rigorous feature $`\psi _{\mathrm{}}/\nu _{}(d=0)=1`$ by supplementing the $`ϵ`$ expansion with an appropriate third order term. The so obtained result reads:

$$\frac{\psi _{\mathrm{}}}{\nu _{}}=1ϵ\left(4ϵ\right)\left\{a_{\mathrm{}}+ϵ\left[b_{\mathrm{}}+c_{\mathrm{}}\mathrm{ln}\left(\frac{4}{3}\right)\right]\right\}+O\left(ϵ^4\right)$$
(4)
where the $`a_{\mathrm{}}`$, $`b_{\mathrm{}}`$ and $`c_{\mathrm{}}`$ are $`\mathrm{}`$-dependent coefficients taking the values listed in Table 2.

In order to estimate the exponents numerically, we determine the current distribution averaging over an ensemble of 5000 independent realizations of a conducting backbone in 3+1, 2+1, and 1+1 dimensions. In principle the current distribution can be computed by solving Kirchhoff’s laws at all nodes of the backbone. In the case of isotropic percolation, where the network consists only of resistors, one obtains a set of linear equations which can be solved exactly. In the present case of a diode-resistor network, however, the resulting equations are nonlinear making it impossible to compute the currents analytically. Thus, in order to compute the currents, we have to apply an iterative approximation procedure. To this end each site of the backbone network is connected to a grounded capacitor $`C`$. For a given set of voltages at the lattice sites, we can compute the actual difference between incoming and outgoing currents. This difference then changes the voltage of the capacitor at site $`i`$ by

$$\frac{dU_i}{dt}=\frac{1}{C}\left(I_{\text{incoming}}I_{\text{outgoing}}\right).$$
(5)
After sufficiently long time this nonlinear dynamical system reaches a stationary state where incoming and outgoing currents balance each other. Therefore, by iterating Eq. (5), the stationary current distribution can be approximated accurately, mainly limited by the available CPU time and the precision of floating-point arithmetics. Moreover, the convergence during the iteration procedure can be controlled easily by monitoring the quantity $`_b|\mathrm{\Delta }I_b|`$, which should tend to zero as the stationary state is approached. Using Newton’s iteration scheme, this quantity actually decreases until it reaches a value of about $`10^{10}`$, which is essentially determined by the limited machine precision. At this point we stop the iteration procedure and compute the moments $`M_{\mathrm{}}`$, assuming that the errors caused by the iteration procedure are much smaller than the statistical and finite-size errors generated by the DP process itself.

In the simulations we observe that the average size of the backbone increases rapidly with $`x_{}^{}x_{}`$, especially in 3+1 dimensions. While the numerical effort for generating conducting backbones between two given points plays a minor role, the algorithm is mainly limited by the CPU time needed for the dynamical approximation and by the memory capacity for storing the configuration of the backbone, the currents, and the voltages. Because of these limitations, we could only go up to distances of $`x_{}^{}x_{}75`$ in 3+1, 2+1 and 1+1 dimensions. Clearly, in this limited range corrections to scaling still play an important role. Nevertheless we find that the moments exhibit very clean power laws (see Fig. 3), which allows us to estimate the critical exponents $`\psi _{\mathrm{}}`$ in the limit $`x_{}^{}x_{}\mathrm{}`$ by standard extrapolation techniques. These estimates are summarized in Table 3.

Figure 4 shows the field-theoretic estimates as given in Eq. (4) compared to the numerical results. The error bars include the statistical error of the simulation combined with the expected error of the extrapolation procedure. In 3+1 dimensions we find excellent agreement. The numerical values and the field theoretic predictions almost coincide. As expected, the agreement is less pronounced in 2+1 dimensions. The analytic estimates lie slightly outside the error bars of the simulations. Though not quantitatively correct, field theory yet predicts the right shape of the dependence of $`\psi _{\mathrm{}}/\nu _{}`$ on $`\mathrm{}`$. In 1+1 dimension, the field theory fails to predict the numerical values of the multifractal exponents accurately. We therefore refrain from showing the respective graph. All in all, our simulations support the validity of the field theoretic results, at least for dimensions to that the $`ϵ`$ expansion can be continued reliably. Our simulations clearly confirm the existence of multifractality in current distributions running on DP clusters.

As in the field-theoretic case, we can check the consistency of our results. For example, as shown in StenullJanssen2 , the zeroth moment in a ($`d`$+1)-dimensional system, which is just the mass of the backbone, can be expressed in terms of the backbone dimension $`D_B=d+12\beta /\nu _{}`$ (cf. Ref. Backbonedim ) by

$$\psi _0=\nu _{}(D_B1+z),$$
(6)
where $`z=\nu _{}/\nu _{}`$. Using the values listed in Table 1, we obtain for $`\psi _0/\nu _{}`$ the values 1.11 in three, 1.23 in two, and 1.307 in one dimension, respectively. These values are in perfect agreement with the estimates listed in Table 3.

To summarize, we have studied the moments of currents running in a random diode network at the critical point by numerical simulations. Our results fully confirm recent field-theoretic results, supporting the claim that currents running on a critical DP cluster indeed show multifractal behavior.

We would like to point out that a random diode network might be a good candidate for an experimental realization of directed percolation. In particular the total resistance $`M_1`$ is easy to measure. In an experimental realization various idealized assumptions of the model shown in Fig. 1 can be relaxed. For example, not all diodes need to be polarized along the distinguished direction, rather it is sufficient for them to be preferentially oriented in a statistical sense. In addition, the resistance of conducting diodes may be stochastically distributed. Nevertheless the critical exponents of such a system at the transition are expected to be the same as in DP. Even inhomogeneities, which cannot be avoided in any experiment, are not expected to play a role since they would result into an annealed noise which is irrelevant under renormalization group transformations (see Ref. Hinrichsen00b ). Because of these favorable properties, random diode networks are particularly suitable for a possible experimental realization of directed percolation. It should be emphasized that so far no experiment is known where the critical exponents of directed percolation could be verified in a reliable way Hinrichsen00b ; Hinrichsen99 . Therefore, the experimental realization of a random diode network would be a particularly challenging task.

H. K. J. acknowledges the support by the Sonderforschungsbereich 237 “Unordnung und große Fluktuationen” of the Deutsche Forschungsgemeinschaft. O. S. is grateful for the support by the Emmy-Noether-Programm of the Deutsche Forschungsgemeinschaft.


text:
# Towards the observation of Hawking radiation in Bose–Einstein condensates.

thanks: Supported by the US NSFthanks: Supported by the US DOE
(6 October 2001; -ed )

## Abstract

Acoustic analogues of black holes (dumb holes) are generated when a supersonic fluid flow entrains sound waves and forms a trapped region from which sound cannot escape. The surface of no return, the acoustic horizon, is qualitatively very similar to the event horizon of a general relativity black hole. In particular Hawking radiation (a thermal bath of *phonons* with temperature proportional to the “surface gravity”) is expected to occur. In this note we consider quasi-one-dimensional supersonic flow of a Bose–Einstein condensate (BEC) in a Laval nozzle (converging-diverging nozzle), with a view to finding which experimental settings could magnify this effect and provide an observable signal. We identify an experimentally plausible configuration with a Hawking temperature of order $`70`$ n K; to be contrasted with a condensation temperature of the order of $`90`$ n K.

Laval nozzle, converging-diverging nozzle, supersonic flow, surface gravity

Acoustic analogues of black holes are formed by supersonic fluid flow unruh ; visser . The flow entrains sound waves and forms a trapped region from which sound cannot escape. The surface of no return, the acoustic horizon, is qualitatively very similar to the event horizon of a general relativity black hole; in particular Hawking radiation (in this case a thermal bath of *phonons* with temperature proportional to the “surface gravity”) is expected to occur unruh ; visser . There are at least three physical situations in which acoustic horizons are known to occur: Bondi–Hoyle accretion bondi-hoyle , the Parker wind parker  (coronal outflow from a star), and supersonic wind tunnels. Recent improvements in the creation and control of Bose–Einstein condensates (see *e.g.*Dalf99 ) have lead to a growing interest in these systems as experimental realizations of acoustic analogs of event horizons. In this note we consider supersonic flow of a BEC through a Laval nozzle (converging-diverging nozzle) in a quasi-one-dimensional approximation. We show that this geometry allows the existence of a fluid flow with acoustic horizons without requiring any special external potential, and we study this flow with a view to finding situations in which the Hawking effect is large. We present simple physical estimates for the “surface gravity” and Hawking temperature footnote .

While writing up this paper, we have found that fluid flow in a Laval nozzle geometry has also been considered in reference sakagami . The occurrence of the same Laval nozzle geometry is the only significant point of overlap of our work with sakagami , as that paper only deals with a classical effect, related to the Hawking effect, but does not consider the quantum physics of Hawking radiation itself. (See padmanabhan  for a previous discussion along these lines.)

*Laval nozzle:* A general problem with the realization of acoustic horizons is that most of the background fluid flows so far studied seem to require very special fine-tuned forms for the external potential in order to be realized (see *e.g.*, the Schwarzschild-like geometry in reference unexpected ). In this respect a possible improvement toward the realizability of acoustic horizons is the construction of a trap which “geometrically constrains” the flow in such a way as to replace the need for a special external potential. An example of such a geometry is the so called Laval nozzle. In particular we shall consider a pair of Laval nozzles; this provides a system which includes a region of supersonic flow bounded between two subsonic regions.

Consider such a nozzle pointing along the $`z`$ axis. Let the cross sectional area be denoted $`A(z)`$. We apply, with appropriate modifications and simplifications, the calculation of unexpected . The crucial approximation is that transverse velocities (in the $`x`$ and $`y`$ directions) are small with respect to velocity along the $`z`$ axis. Then, assuming steady flow, we can write the continuity equation in the form

$$\rho (z)A(z)v(z)=J;J=\text{constant.}$$
(1)
The Euler equation (including for the moment possible external body forces $`\mathrm{d}\mathrm{\Phi }/\mathrm{d}z`$, and internal viscous friction $`f_v`$) reduces to

$$\rho v\frac{\mathrm{d}v}{\mathrm{d}z}=\frac{\mathrm{d}p}{\mathrm{d}z}\rho \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}z}+f_v.$$
(2)
Finally, we assume a barotropic equation of state $`\rho =\rho (p)`$, and define $`X^{}=\mathrm{d}X/\mathrm{d}z`$. Then continuity implies

$$\rho ^{}=\rho \frac{(Av)^{}}{(Av)}=\rho \left[\frac{A^{}}{A}+\frac{v^{}}{v}\right]=\rho \left[\frac{A^{}}{A}+\frac{a}{v^2}\right],$$
(3)
while Euler implies

$$\rho a=\frac{\mathrm{d}p}{\mathrm{d}\rho }\rho ^{}\rho \mathrm{\Phi }^{}+f_v.$$
(4)
Defining the speed of sound by $`c^2=\mathrm{d}p/\mathrm{d}\rho `$, and eliminating $`\rho ^{}`$ between these two equations yields a form of the well-known “nozzle equation”

$$a=\frac{v^2}{c^2v^2}\left(c^2\left[\frac{A^{}}{A}\right]\mathrm{\Phi }^{}+\frac{f_v}{\rho }\right)$$
(5)
The presence of the factor $`c^2v^2`$ in the denominator is crucial and leads to several interesting physical effects. For instance, if the physical acceleration is to be finite at the acoustic horizon, we need

$$c^2\left[\frac{A^{}}{A}\right]\mathrm{\Phi }^{}+\frac{f_v}{\rho }\mathrm{\hspace{0.33em}0}.$$
(6)
This is a fine-tuning condition relating the shape of the nozzle (the amount of convergence/divergence at the horizon) to the external body force and specific friction. Experience with wind tunnels has shown that the flow will attempt to self-adjust (in particular, the location of the acoustic horizon will self-adjust) so as to satisfy this fine tuning. We could now calculate the acceleration of the fluid at the acoustic horizon by adopting L’Hôpital’s rule. However, it is more useful to consider the “surface gravity” defined by the limit of the quantity visser

$$g=\frac{1}{2}\frac{d(c^2v^2)}{\mathrm{d}z}$$
(7)
It is this combination $`g`$, rather than the physical acceleration of the fluid $`a`$, that more closely tracks the notion of “surface gravity”, and it is the limit of this quantity as one approaches the acoustic horizon that enters into the Hawking radiation calculation essential . To calculate the limit $`g_H`$ we need to use

$`(c^2)^{}`$ $``$ $`{\displaystyle \frac{\mathrm{d}^2p}{\mathrm{d}\rho ^2}}\rho ^{}=\rho {\displaystyle \frac{\mathrm{d}^2p}{\mathrm{d}\rho ^2}}{\displaystyle \frac{(Av)^{}}{Av}}`$ (8)
$``$ $`\rho _H{\displaystyle \frac{\mathrm{d}^2p}{\mathrm{d}\rho ^2}}|_H\left[{\displaystyle \frac{A_H^{}}{A_H}}+{\displaystyle \frac{a_H}{c_H^2}}\right].`$
This implies that the fine-tuning (6) used to keep $`a`$ finite at the acoustic horizon will also keep $`g`$ finite there. In particular, we can use L’Hôpital’s rule at the acoustic horizon to evaluate

$$g_H=\frac{\rho }{2}\frac{\mathrm{d}^2p}{\mathrm{d}\rho ^2}\frac{A^{}}{A}|_H+\frac{\left[c^2+\frac{\rho }{2}\frac{\mathrm{d}^2p}{\mathrm{d}\rho ^2}\right]\left[c^2\frac{A^{}}{A}\mathrm{\Phi }^{}+\frac{f_v}{\rho }\right]^{}|_H}{2g_H}$$
(9)
This result, which we quote for its generality, is actually more complicated than we really need.

*Free flow:* The Laval nozzle corresponds to the special case in which there is no external body force $`\mathrm{\Phi }^{}=0`$. Additionally, when considering a superfluid flow, there is no viscous friction $`f_v=0`$. Then the nozzle equation (5) reduces to

$$a=\frac{c^2v^2}{c^2v^2}\left[\frac{A^{}}{A}\right].$$
(10)
Regularity now requires the much simpler fine-tuning condition that $`A^{}=0`$ at the horizon. That is, the acoustic horizon occurs at a point of minimum area; exactly the behaviour which is physically seen in a Laval nozzle. (That a horizon cannot form at a maximum of the cross sectional area is established below.) Now apply the L’Hôpital rule at the horizon (using the fact that at this point $`A^{}=0`$)

$$a_H=\frac{c^4A^{\prime \prime }/A}{(c^2)^{}2a_H}|_H.$$
(11)
But from the specialization of (8) we now see

$$a_H^2=\frac{c^4A^{\prime \prime }/A}{2+\rho (\mathrm{d}^2p/\mathrm{d}\rho ^2)/c^2}|_H.$$
(12)
That is

$$a_H=\pm \frac{c^2}{\sqrt{2A}}\sqrt{\frac{A^{\prime \prime }}{1+(1/2)\rho [\mathrm{d}^2p/\mathrm{d}\rho ^2]/c^2}}|_H.$$
(13)
Very similar formulae hold for $`g_H`$:

$$g_H^2=+\frac{1}{2}\left[c^2+\frac{\rho }{2}\frac{\mathrm{d}^2p}{\mathrm{d}\rho ^2}\right]_Hc^2\left[\frac{A_H^{\prime \prime }}{A_H}\right],$$
(14)
and so

$$g_H=\pm \frac{c_H^2}{\sqrt{2A_H}}\sqrt{1+\frac{\rho }{2c^2}\frac{\mathrm{d}^2p}{\mathrm{d}\rho ^2}}|_H\sqrt{A_H^{\prime \prime }}.$$
(15)
The first factor is of order $`c_H^2/R`$, with $`R`$ the minimum radius of the nozzle, while the second and third factors are square roots of dimensionless numbers. This is in accord with our intuition based on dimensional analysis visser ; unexpected . If $`A^{\prime \prime }<0`$, corresponding to a maximum of the cross section, then $`g_H`$ is imaginary which means no event horizon can form there. The two signs $`\pm `$ correspond to either speeding up and slowing down as you cross the horizon, both of these must occur at a minimum of the cross sectional area $`A^{\prime \prime }>0`$. (If the flow accelerates at the horizon this is a black hole horizon \[future horizon\]; if the flow decelerates there it is a white hole horizon \[past horizon\]. See Figure 1.) If the nozzle has a circular cross section, then the quantity $`A_H^{\prime \prime }`$ is related to the longitudinal radius of curvature $`R_c`$ at the throat of the nozzle, in fact $`A_H^{\prime \prime }=\pi R/R_c`$.

*Bose–Einstein condensate:* The use of BECs as a working fluid for acoustic black holes has been advocated by Garay *et al.* garay  (see also barcelo  for a discussion of plausible extensions to that model). The present note can be interpreted as a somewhat different approach to the same physical problem, side-stepping the technical complications of the Bogoliubov equations in favour of a more fluid dynamical point of view. For a standard BEC

$$c^2=\frac{\lambda \rho }{m}.$$
(16)
Then

$$\rho \left[\frac{\mathrm{d}^2p}{\mathrm{d}\rho ^2}\right]=\rho \frac{\mathrm{d}(c^2)}{\mathrm{d}\rho }=c^2,$$
(17)
while

$$1+\frac{1}{2}\frac{\rho }{c^2}\left[\frac{\mathrm{d}^2p}{\mathrm{d}\rho ^2}\right]=\frac{3}{2}.$$
(18)
So we have, rather simply

$$a_H=\pm \frac{c_H^2}{\sqrt{A_H}}\sqrt{A_H^{\prime \prime }/3}.$$
(19)
Similarly

$$g_H=\pm \frac{c_H^2}{\sqrt{A_H}}\sqrt{3A_H^{\prime \prime }/4}.$$
(20)
This implies, at a black hole horizon \[future horizon\], a Hawking temperature unruh ; visser ; essential

$$k_\mathrm{B}T_H=\frac{\mathrm{}g_H}{2\pi c_H}=\mathrm{}\frac{c_H}{2\pi \sqrt{A_H}}\sqrt{\frac{3A_H^{\prime \prime }}{4}}.$$
(21)
Ignoring for now the issue of gray-body factors (they are a refinement on the Hawking effect, not really an essential part of the physics), the phonon spectrum peaks at

$$\omega _{\mathrm{peak}}=\frac{c_H}{2\pi \sqrt{A_H}}\sqrt{\frac{3A_H^{\prime \prime }}{4}},$$
(22)
that is

$$\lambda _{\mathrm{peak}}=4\pi ^2\sqrt{A_H}\sqrt{\frac{4}{3A_H^{\prime \prime }}}.$$
(23)
This extremely simple result relates the Hawking emission to the physical size of the constriction and a factor depending on the flare-out at the narrowest point. Note that you cannot permit $`A_H^{\prime \prime }`$ to become large, since then you would violate the quasi-one-dimensional approximation for the fluid flow that we have been using in this note. (There is of course nothing physically wrong with violating the quasi-one-dimensional approximation, it just means the analysis becomes more complicated. In particular, if there is no external body force and the viscous forces are zero then by slightly adapting the analysis of unexpected  the acoustic horizon \[more precisely the ergo-surface\] is a minimal surface of zero extrinsic curvature.) The preceding argument suggests strongly that the best we can realistically hope for is that the spectrum peaks at wavelength

$$\lambda _{\mathrm{peak}}\sqrt{A_H}.$$
(24)
You can (in principle) try to adjust the equation of state to make the second factor in (15) larger, but this is unlikely to be technologically feasible.

(Note that this is the analog, in the context of acoustic black holes, of the fact that the Hawking flux from general relativity black holes is expected to peak at wavelengths near the physical diameter of the black hole, its Schwarzschild radius — up to numerical factors depending on charge and angular momentum.)

*Discussion:* It is the fact that the peak wavelength of the Hawking radiation is of order the physical dimensions of the system under consideration that makes the effect so difficult to detect. In particular, in BECs it is common to have a sound speed of order $`6\text{mm/s}`$. If one then chooses a nozzle of diameter about 1 micron, and a flare-out of $`A_H^{\prime \prime }1`$, then $`T_H7\text{n K}`$. Compare this to the condensation temperature required to form the BEC

$$T_{\mathrm{condensate}}90\text{n K}.$$
(25)
We see that in this situation the Hawking effect, although tiny, is at least comparable in magnitude to other relevant temperature scales. Moreover recent experiments indicate that it is likely that these figures can be improved. In particular, the scattering length for the condensate can be tuned by making use of the so called Feshbach resonance feshbach . This effect can be used to increment the scattering length; factors of up to 100 have been experimentally obtained feshbach2 . Therefore the acoustic propagation speed, which scales as the square root of the scattering length, could thereby be enhanced by a factor up to 10. This suggests that it might be experimentally possible to achieve $`c_H6\text{cm/s}`$, and so

$$T_H70\text{n K},$$
(26)
which places us much closer to the condensation temperature. The speed of sound can also be enhanced by increasing the density of the condensate (propagation speed scales as the square root of the density). In all of these situations there is a trade-off: For fixed nozzle geometry the Hawking temperature scales as the speed of sound, so larger sound speed gives a bigger effect but conversely makes it more difficult to set up the supersonic flow.

So if the Hawking effect can be experimentally realized in these situations, it may be sufficiently large to disrupt the condensate configuration, or even the condensate itself, providing in this way a clear signal. Let us elaborate this point: The Garay *et al.* analysis garay  shows that it should be possible to create classically stable BEC configurations with the presence of acoustic horizons. In their particular analysis they found that these classically stable configurations are surrounded by unstable regions; (at this stage, we don’t know how general that result is). Once the system is engineered to be in a classically stable configuration one can look for the purely quantum effect of Hawking emission. The power loss due to Hawking radiation would be:

$$P=\sigma T_H^4A_H=\frac{3\mathrm{}c_H^4}{5120\pi ^2A_H}(A_H^{\prime \prime })^2.$$
(27)
Numerically (including the effect of the Feshbach resonance), the emitted power is extremely small $`P10^{48}\text{ W}`$, but it should be noted that some sizable fraction of the Hawking phonons will be sufficiently energetic to knock atoms out of the condensate phase.

Additionally, we want to stress that the present analysis of the BEC is purely “hydrodynamic” (superfluid approximation of the BEC), and does not seek to deal with the “quantum potential”(see *e.g.*Dalf99 ; barcelo ). The latter is responsible for the fact that the dispersion relation for perturbations in the BEC is modified at high momenta in such a way as to recover “infinite” propagation speed (this is the so-called Bogoliubov dispersion relation barcelo ; cpt01 ). This issue has relevance to the so-called trans-Planckian problem (which in this BEC condensate context becomes a trans-Bohrian problem). Fortunately it is known, thanks to model calculations in field theories with explicit high-momentum cutoffs, that the low energy physics of the emitted radiation is largely insensitive to the nature and specific features of the cutoff. However, because of this high-frequency “superluminal” dispersion, one can have additional hopes of detecting a signal because Hawking radiation could provide a new type of instability, disrupting the classically stable configuration by leading to a run-away production of phonons (see the discussion of “black hole lasers” in CT ).

To summarize what we have learned: The present note complements the analysis by Garay *et al.* garay , in that it provides a rationale for simple physical estimates of the Hawking radiation temperature without having to solve the full Bogoliubov equations. Additionally, the current analysis provides simple numerical estimates of the size of the effect and identifies several specific physical mechanisms by which the Hawking temperature can be manipulated: via the speed of sound, the nozzle radius, the equation of state, and the degree of flare-out at the throat. In this manner, we have identified an experimentally plausible configuration with a Hawking temperature of order $`70`$ n K; to be contrasted with a condensation temperature of the order of $`90`$ n K.

*Acknowledgements:* S. Liberati wishes to thank B.L. Hu and T.A. Jacobson for illuminating discussions and remarks. M. Visser wishes to thank M. Stone for stimulating comments regarding the nozzle equation.