text:
# A square-well model for the structural and thermodynamic properties of simple colloidal systems

## I Introduction

In the simplest model of a colloidal dispersion, the interactions among the (large) solute molecules and the excluded volume effects of the (small) solvent molecules, which lead to intercolloidal solvation forces, are incorporated by treating the solute particles as hard spheres (HS) of diameter $`\sigma `$. This allows one to take advantage of the analytical solution of the Percus-Yevick (PY) integral equation for the HS potential. On the other hand, it is well known that the solvent particles (e.g., macromolecules) can induce an effective short-range attraction between two colloidal particles due to an unbalanced osmotic pressure arising from depletion of the solvent particles in the region between the two colloidal ones.  This explains the widespread use of Baxter’s model of sticky hard spheres (SHS) to represent the properties of colloidal suspensions.  This interaction model represents the attractive part of the potential as an infinitely deep, infinitely narrow well. In addition to the hard-core diameter $`\sigma `$, the model incorporates a “stickyness” parameter $`\tau _{\text{SHS}}^1`$ (essentially the deviation of the second virial coefficient from the HS value) that can be understood as a measure of the temperature: the smaller the temperature the larger the stickyness. The SHS potential has the advantage of lending itself to an analytical solution in the PY approximation and in the Mean Spherical Approximation (MSA).  On the other hand, the SHS potential presents two limitations. First, the system of monodisperse SHS is not thermodynamically stable; however this pathology, which is not captured by Baxter’s solution to the PY equation, can be remedied by including some degree of polydispersity in the system. The most important limitation of the SHS model as a representation of a realistic short-range interaction lies in the fact that it is unable to distinguish two situations with the same “stickyness” (i.e. same second virial coefficient) and density, but different temperature and/or range.

In order to mimic the particle interaction in colloidal systems in a more realistic way, a simple choice is to assume that the particles interact via the square-well (SW) potential

$$\phi (r)=\{\begin{array}{cc}\mathrm{},\hfill & r<\sigma \hfill \\ ϵ,\hfill & \sigma <r<\lambda \sigma \hfill \\ 0,\hfill & r>\lambda \sigma ,\hfill \end{array}$$
(1)
where $`\sigma `$ is the diameter of the hard core, $`ϵ`$ is the well depth, and $`(\lambda 1)\sigma `$ is the well width. The equilibrium properties of a SW fluid depend on the values of three dimensionless parameters: the reduced number density $`\rho ^{}=\rho \sigma ^3`$, the reduced temperature $`T^{}=k_BT/ϵ`$ ($`k_B`$ being the Boltzmann constant), and the width parameter $`\lambda `$. In the limits $`ϵ0`$ (i.e. $`T^{}\mathrm{}`$) and/or $`\lambda 1`$, the SW fluid becomes the HS fluid. In addition, the SHS fluid  is obtained by taking the limits $`ϵ\mathrm{}`$ (i.e. $`T^{}0`$) and $`\lambda 1`$, while keeping the stickyness parameter $`\tau _{\text{SHS}}^1=12(\lambda 1)(e^{1/T^{}}1)`$ constant. If we define in the original SW fluid a generalized stickyness parameter $`\tau _{\text{SW}}^1=4(\lambda ^31)\left(e^{1/T^{}}1\right)`$ as being proportional to the deviation of the second virial coefficient from the HS value, then at a given density $`\rho ^{}`$ the parameter space can be taken as the plane $`(\tau _{\text{SW}}^1,\lambda )`$. The SHS limit explores the line $`(\tau _{\text{SW}}^1,\lambda =1)`$ only, while the HS model corresponds to the point $`(\tau _{\text{SW}}^1,\lambda )=(0,1)`$. This geometrical picture illustrates why the SW interaction model can be useful to uncover a much richer spectrum of values for the relevant parameters of the problem, even if the attraction range is relatively short.

Despite the mathematical simplicity of the SW potential, no analytical solution of the conventional integral equations for fluids (YBG, HNC, PY, …) is known.  The mean spherical model approximation of Sharma and Sharma provides an analytical expression for the structure factor, but it is not consistent with the hard core exclusion constraint. Most of the available theoretical information about the SW fluid for variable width comes from perturbation theory. In general, perturbation theory is based on an expansion of the relevant physical quantities in powers of the inverse temperature. For instance, the radial distribution function $`g(r;\rho ^{},T^{})`$ of the SW fluid is expressed as

$$g(r;\rho ^{},T^{})=g_0(r;\rho ^{})+T_{}^{}{}_{}{}^{1}g_1(r;\rho ^{})+\mathrm{},$$
(2)
where $`g_0(r;\rho ^{})`$ is the radial distribution function of the reference HS fluid and $`g_1(r;\rho ^{})`$ represents a first order correction. Good analytical approximations for $`g_0`$ are known, such as Wertheim-Thiele’s solution of the PY equation,  Verlet-Weiss parameterization,  or the Generalized Mean Spherical Approximation.  Thus, it is in the choice of $`g_1`$ where different versions of perturbation theory essentially differ. A few years ago, Tang and Lu (TL) proposed an analytical expression (in Laplace space) for $`g_1`$, based on the MSA. Comparison with Monte Carlo (MC) simulation data  showed a general good agreement for the cases considered, but the quality of the agreement worsened as smaller values of $`\lambda `$ and/or $`T^{}`$ were taken. The expectation that perturbation theory becomes less accurate as the well width and the temperature decrease has been already reported elsewhere. In fact, perturbation theory tends to overestimate the critical temperature $`T_c^{}`$, especially for small values of $`\lambda 1`$. All these limitations are significantly apparent in the SHS limit ($`\lambda 1`$, $`T^{}0`$), in which case the expansion (2) becomes meaningless. By following a completely different approach, Nezbeda  proposed to approximate $`rg(r)`$ in the interval $`\sigma r\lambda \sigma `$ by a quadratic function of $`r`$, within the context of the PY theory. The coefficients of the polynomial were then determined analytically by imposing the continuity of the cavity function $`y(r)g(r)e^{\phi (r)/k_BT}`$ and its first two derivatives at $`r=\sigma `$. Notwithstanding the merits of this theory, its main limitation is that it is only applicable for very narrow wells, typically $`\lambda 10.01`$. In addition, Nezbeda’s theory fails to predict a thermodynamic critical point, except at $`\lambda =1`$, in which case Baxter’s solution for SHS is recovered.

In order to provide a simple theory that, while including the SHS case as a special limit, could also be applicable to SW fluids of variable width, Yuste and Santos proposed a model by assuming an explicit functional form for the Laplace transform $`G(t)`$ of $`rg(r)`$. That functional form was suggested by the exact virial expansion of the radial distribution function  and by the property $`lim_{r\sigma ^+}g(r)=\text{finite}`$. The parameters were subsequently determined as functions of $`\rho ^{}`$, $`T^{}`$, and $`\lambda `$ by imposing the condition $`S(0)=\text{finite}`$, where $`S(q)`$ is the structure function, as well as the continuity of $`y(r)`$ at $`r=\lambda \sigma `$. The structural properties predicted by the model showed a good agreement with MC simulation results, not only for narrow square wells, but even for relatively wide ones ($`\lambda 1.5`$) up to densities slightly above the critical density. On the other hand, the model was not fully analytical because the continuity condition of $`y(r)`$ at $`r=\lambda \sigma `$ gives rise to a transcendent equation that must be solved numerically. In fact, the exact solution in the case of one-dimensional SW fluids involves a similar transcendent equation.

The aim of this paper is to propose a simpler version of the model introduced by Yuste and Santos in Ref. . While we keep the same functional structure of $`G(t)`$ and enforce the conditions $`g(\sigma ^+)=\text{finite}`$ and $`S(0)=\text{finite}`$, we replace the transcendent equation stemming from the continuity of $`y(r)`$ at $`r=\lambda \sigma `$ by a quadratic equation suggested by the SHS limit. The resulting model has therefore a degree of algebraic simplicity similar to that of the PY solution for SHS (but now the parameters have a $`\lambda `$-dependence beyond the one captured by the stickyness coefficient) and reduces to it in the appropriate limit. The structural properties $`g(r)`$ and $`S(q)`$ exhibit a fairly good agreement with MC simulations,  similar to that found in the original, more complicated version of the model. Nevertheless, it is in the calculation of the thermodynamic properties (which were not addressed in Ref. ) where the present model becomes especially advantageous. In particular, the isothermal compressibility is obtained as an explicit function of density, temperature, and well width. This allows us to get the $`\lambda `$-dependence of the critical temperature $`T_c^{}`$ and density $`\rho _c^{}`$. Comparison with computer simulation estimates  shows that the model predictions for $`T_c^{}(\lambda )`$ are remarkably good, even for values of $`\lambda `$ as large as $`\lambda =1.75`$. On the other hand, the predicted values of $`\rho _c^{}(\lambda )`$ are typically 30–45% smaller than the simulation ones, a fact that can be traced back to the solution of the PY equation for SHS. We also compare the predictions of the model for the compressibility factor at $`\lambda =1.125`$ and $`\lambda =1.4`$ with simulation data  and with the TL perturbation theory. In both cases the model presents a better agreement than the perturbation theory, except for temperatures larger than about twice the critical temperature.

The paper is organized as follows. The model is introduced and worked out in Sec. II, with some technicalities being relegated to the appendices. Section III deals with the comparison with simulation results and perturbation theory. The paper ends with a brief discussion in Sec. IV

## II The model

### A Basic requirements

The radial distribution function $`g(r)`$ of a fluid is directly related to the probability of finding two particles separated by a distance $`r`$. It can be measured from neutron or x-ray diffraction experiments through the static structure factor $`S(q)`$. Both quantities are related by Fourier transforms:

$`S(q)`$ $`=`$ $`1+\rho {\displaystyle 𝑑𝐫e^{i𝐪𝐫}[g(r)1]}`$ (3)
$`=`$ $`12\pi \rho {\displaystyle \frac{G(t)G(t)}{t}}|_{t=iq},`$ (4)
where $`\rho `$ is the number density and

$$G(t)=_0^{\mathrm{}}𝑑re^{rt}rg(r)$$
(5)
is the Laplace transform of $`rg(r)`$. The isothermal compressibility of the fluid, $`\kappa _T=\rho ^1\left(\rho /p\right)_T`$, is directly related to the long-wavelength limit of the structure function:

$$\chi _T\rho k_BT\kappa _T=S(0).$$
(6)
Thus, all the physically relevant information about the equilibrium state of the system is contained in $`g(r)`$ or, equivalently, in $`G(t)`$.

Now we particularize to the SW interaction potential, Eq. (1), so $`g(r)=0`$ for $`r<\sigma `$. Henceforth, we will take the hard-core diameter $`\sigma =1`$ as the length unit and the well depth $`ϵ/k_B=1`$ as the temperature unit, so that the asterisks in $`\rho ^{}`$ and $`T^{}`$ will be dropped. It is convenient to define an auxiliary function $`F(t)`$ through the relation

$`G(t)`$ $`=`$ $`t{\displaystyle \frac{F(t)e^t}{1+12\eta F(t)e^t}}`$ (7)
$`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(12\eta )^{n1}t[F(t)]^ne^{nt},`$ (8)
where $`\eta \frac{\pi }{6}\rho \sigma ^3`$ is the packing fraction. Laplace inversion of Eq. (7) provides an useful representation of the radial distribution function:

$$g(r)=r^1\underset{n=1}{\overset{\mathrm{}}{}}(12\eta )^{n1}f_n(rn)\mathrm{\Theta }(rn),$$
(9)
where $`f_n(r)`$ is the inverse Laplace transform of $`t[F(t)]^n`$ and $`\mathrm{\Theta }(r)`$ is Heaviside’s step function. Thus, the knowledge of $`F(t)`$ is fully equivalent to that of $`g(r)`$ or $`S(q)`$. In particular, the value of $`g(r)`$ at contact, $`g(1^+)`$, is given by the asymptotic behavior of $`F(t)`$ for large $`t`$:

$$g(1^+)=f_1(0)=\underset{t\mathrm{}}{lim}t^2F(t).$$
(10)
Since $`g(1^+)`$ must be finite and different from zero, we get the condition

$$F(t)t^2,t\mathrm{}.$$
(11)
On the other hand, according to Eq. (3), the behavior of $`G(t)`$ for small $`t`$ determines the value of $`S(0)`$:

$$G(t)=t^2+\text{const}+\frac{1S(0)}{24\eta }t+o(t^2).$$
(12)
Insertion of Eq. (12) into the first equality of Eq. (7) yields the first five terms in the expansion of $`F(t)`$ in powers of $`t`$:

$$F(t)=\frac{1}{12\eta }\left(1+t+\frac{1}{2}t^2+\frac{1+2\eta }{12\eta }t^3+\frac{2+\eta }{24\eta }t^4\right)+𝒪(t^5).$$
(13)
The value of $`S(0)`$ is related to the coefficients of $`t^5`$ and $`t^6`$ by

$$S(0)=\frac{24}{5}\eta ^3\left[6\frac{d^5F(t)}{dt^5}|_{t=0}\frac{d^6F(t)}{dt^6}|_{t=0}\right]1+8\eta +2\eta ^2.$$
(14)
So far, all the expressions apply to any density and any hard-core potential. As is well known, in the limit of zero density the cavity function $`y(r)g(r)e^{\phi (r)/k_BT}`$ is equal to 1. In the special case of the SW potential this translates into

$$\underset{\eta 0}{lim}g(r)=(1+x)\mathrm{\Theta }(r1)x\mathrm{\Theta }(r\lambda ),$$
(15)
where $`xe^{1/T}1`$. Equation (15) implies that

$$\underset{\eta 0}{lim}F(t)=(1+x)(t^2+t^3)xe^{(\lambda 1)t}(\lambda t^2+t^3).$$
(16)
### B Construction of the model

Any meaningful approximation of $`F(t)`$ for the SW potential must comply with Eqs. (11), (13), and (16). Let us decompose $`F(t)`$ as

$$F(t)=R(t)\overline{R}(t)e^{(\lambda 1)t}.$$
(17)
The model proposed by Yuste and Santos  consists of assuming the following rational forms for $`R(t)`$ and $`\overline{R}(t)`$:

$$R(t)=\frac{A_0+A_1t}{1+S_1t+S_2t^2+S_3t^3},\overline{R}(t)=\frac{\overline{A}_0+\overline{A}_1t}{1+S_1t+S_2t^2+S_3t^3}.$$
(18)
This forms are compatible with (16). In addition, condition (11) is satisfied by construction. In fact, the contact value is, according to Eq. (10),

$$g(1^+)=\frac{A_1}{S_3}.$$
(19)
In order to ease the proof that the HS and SHS cases are included in Eq. (18), it is convenient to introduce the new parameters $`A12\eta \overline{A}_0`$, $`L_112\eta \left[\overline{A}_0(\lambda 1)+A_1\overline{A}_1\right]`$, and $`L_212\eta \overline{A}_1(\lambda 1)`$. With these changes, Eqs. (17) and (18) can be recast into

$$F(t)=\frac{1}{12\eta }\frac{1+A+[L_1+L_2(\lambda 1)^1A(\lambda 1)]t[A+L_2(\lambda 1)^1t]e^{(\lambda 1)t}}{1+S_1t+S_2t^2+S_3t^3},$$
(20)
where we have already taken into account the property $`F(0)=1/12\eta `$, according to Eq. (13). The model (20) contains six parameters to be determined. The exact expansion (13) imposes four constraints among them. Thus, we can express four of the parameters in terms of, for instance, $`A`$ and $`L_2`$. The result is

$$L_1=\frac{1}{1+2\eta }\left[1+\frac{1}{2}\eta +2\eta (1+\lambda +\lambda ^2)L_2\frac{1}{2}\eta (3+2\lambda +\lambda ^2)A^{}\right],$$
(21)
$$S_1=\frac{\eta }{1+2\eta }\left[\frac{3}{2}+2(1+\lambda +\lambda ^2)L_2\frac{1}{2}(3+2\lambda +\lambda ^2)A^{}\right],$$
(22)
$$S_2=\frac{1}{2(1+2\eta )}\left\{1+\eta +2[12\eta \lambda (1+\lambda )]L_2[1\eta (1+\lambda )^2]A^{}\right\},$$
(23)
$$S_3=\frac{1}{1+2\eta }\left\{\frac{(1\eta )^2}{12\eta }\left[\frac{1}{2}(\lambda +1)\eta \lambda ^2\right]L_2+\frac{1}{12}[4+2\lambda \eta (3\lambda ^2+2\lambda +1)]A^{}\right\},$$
(24)
where $`A^{}A(\lambda 1)^2`$. These four parameters are linear functions of $`A`$ and $`L_2`$. Taking into account Eq. (14), the value of $`S(0)`$ can be expressed as a quadratic function of $`A`$ and $`L_2`$. Its explicit expression is given in Appendix A.

Two additional constraints are still needed to determine $`A`$ and $`L_2`$. First note that in the zero-density limit we have $`L_11`$, $`S_10`$, $`S_2\text{finite}`$, $`S_3(12\eta )^1`$. Thus, Eq. (20) is consistent with Eq. (16) provided that

$$\underset{\eta 0}{lim}A=x,\underset{\eta 0}{lim}L_2=x\lambda (\lambda 1).$$
(25)
In the original formulation of the model, the parameter $`A`$ was assumed to be independent of density, so that $`A=x`$. As for $`L_2`$, it was determined by imposing the (exact) continuity condition of the function $`y(r)`$ at $`r=\lambda `$, which implies

$$g(\lambda ^{})=(1+x)g(\lambda ^+).$$
(26)
The implementation of this condition in the model (20) leads to a transcendent equation that must be solved numerically. In this paper we will be concerned with a simpler version of the model in which the strong condition (26) is replaced by a weaker one. To that end, let us introduce a parameter $`\tau `$ as

$$\tau \frac{1}{12}\left(\frac{L_1}{L_2}\frac{S_2}{S_3}\right).$$
(27)
From (25) and (27) it follows that

$$\underset{\eta 0}{lim}\tau =[12x\lambda (\lambda 1)]^1.$$
(28)
The definition (27) is suggested by the fact that, as proved in Appendix B of Ref. , Eq. (26) is equivalent to $`\tau =\tau _{\text{SHS}}`$ in the SHS limit, namely $`x\mathrm{}`$, $`\lambda 1`$, $`x(\lambda 1)(12\tau _{\text{SHS}})^1=\text{finite}`$. Therefore, we may expect that a simple prescription for $`\tau `$ (such that $`\tau ^112x(\lambda 1)`$ in the SHS limit) can be a good substitute for the transcendent equation arising from (26), at least for wells relatively narrow.

Since $`L_1`$, $`S_2`$, and $`S_3`$ are linear functions of $`L_2`$, Eq. (27) shows that $`\tau `$ is the ratio of two quadratic functions of $`L_2`$. This relation is easily inverted to get

$$L_2=\frac{1+\alpha _1\eta +\alpha _2\eta ^2+\alpha _3\eta ^3+(1+2\eta )\left[1+\beta _1\eta +\beta _2\eta ^2+\beta _3\eta ^3+\beta _4\eta ^4\right]^{1/2}}{12\eta (\gamma _0+\gamma _1\eta +\gamma _2\eta ^2)},$$
(29)
where the expressions of the coefficients $`\alpha _i`$, $`\beta _i`$, and $`\gamma _i`$ as functions of $`\lambda `$, $`A`$, and $`\tau `$ are given in Appendix A.

So far, we are free to fix $`\tau `$ as a function of $`\eta `$, $`x`$, and $`\lambda `$ by following any criterion we wish. In particular, we can enforce Eq. (26), as done in Ref. . Analogously, the parameter $`A(\eta ,x,\lambda )`$ (which was taken $`A=x`$ in Ref. ) can be freely chosen. In the simplified version of the model we consider here, we assume that both $`A`$ and $`\tau `$ are independent of density, so they take values needed to satisfy the requirement (16). Those values are simply

$$A=x,\tau =[12x\lambda (\lambda 1)]^1.$$
(30)
Note that this $`\tau `$ is only slightly different from the parameter $`\tau _{\text{SW}}`$ introduced below Eq. (1), both of them becoming identical to $`\tau _{\text{SHS}}`$ in the limit $`\lambda 1`$. The choice (30) closes the construction of the model. The pair distribution function in Laplace space is given by Eqs. (7) and (20), where the expressions for the coefficients are (21)–(24), (29), and (30). The model provides the quantity $`G(t)`$ as an explicit function of the Laplace variable $`t`$, the packing fraction $`\eta `$, the well width $`\lambda `$, and the temperature parameter $`xe^{1/T}1`$. Since the poles of $`F(t)`$ are the roots of a cubic equation, the inverse Laplace transforms of $`t[F(t)]^n`$ are analytically derived and then the radial distribution function is readily obtained from the representation (9). From Eq. (30) it follows that the relationship between the temperature and the parameter $`\tau `$ is

$$T=1/\mathrm{ln}[1+\tau ^1/12\lambda (\lambda 1)].$$
(31)
The predictions of our model to first order in the packing fraction are compared with the exact results in Appendix B. As an illustrative example, Fig. 1 shows the quantity

$$\mathrm{\Delta }(r)=\underset{\eta 0}{lim}\frac{1}{\eta }\frac{g(r)g_{\text{exact}}(r)}{g_{\text{exact}}(r)}$$
(32)
for $`\lambda =1.1`$ and $`\lambda =1.125`$ and for the temperatures $`T=0.5`$, 0.67, and 1. The function $`\mathrm{\Delta }(r)`$ is different from zero in the interval $`1<r<\lambda `$ only, where our model slightly overestimates the value of $`g(r)`$. Note that the relative deviation of our $`g(r)`$ from the exact radial distribution function is $`\mathrm{\Delta }(r)\eta `$ for small densities. Of course, the linear growth of this deviation with increasing density is valid in this low-density regime only, as comparison with simulation values at finite densities shows \[cf. Figs. 4 and 5\].

Appendix C shows that the model reduces to Wertheim-Thiele’s and Baxter’s analytical solutions of the PY equation in the HS and SHS limits, respectively. If we consider the parameter space $`\lambda `$$`\tau ^1`$, then the HS potential corresponds to the line $`\tau ^1=0`$ (in which case the physical properties are independent of $`\lambda `$), while the SHS potential corresponds to the line $`\lambda =1`$ (the physical properties being $`\tau `$-dependent). What our model does is to extend the above picture to the entire plane $`\tau ^10`$, $`\lambda 1`$, without compromising the mathematical simplicity present in the analytical solutions of the PY integral equation for HS and SHS.

## III Comparison with simulation and other theories

### A Structural properties

The structural properties obtained from the original version of the model \[i.e. with $`L_2`$ determined by solving the transcendent equation stemming from Eq. (26)\] were profusely compared with simulation data  in Ref. . We have checked that the simplified version of the model presented here \[cf. (29), (30)\] gives results very close to those of the original model, especially for narrow wells. Consequently, we will present only a brief comparison with simulations in this Subsection.

The structure of a fluid is usually determined by neutron or x-ray scattering experiments, which measure the structure factor $`S(q)`$. This quantity is directly related to the Laplace transform $`G(t)`$ by Eq. (3) and so can be obtained explicitly in our model. In 1984, Huang et al.  performed MC simulations of SW fluids with $`\lambda =1.1`$ to reproduce the main features of the structure factor of micellar suspensions. Figure 2 shows $`S(q)`$ obtained from simulation  for $`\lambda =1.1`$, $`T^1=0.92`$, and $`\eta =0.07`$, as compared with our model, the PY solution for SHS (with the conventional choice of $`\tau _{\text{SW}}^1`$ as stickyness parameter), the Tang-Lu (TL) perturbation theory, and the PY solution for hard spheres. The deviations of the simulation data from the PY-HS curve are essentially a measure of effects associated with the non-zero values of the square well width and of the inverse temperature. These are qualitatively described by the TL perturbation theory, except for small wavenumbers. This region is well represented by the PY-SHS curve and by our model, but the latter is better near the first maximum. The value of the critical temperature for $`\lambda =1.1`$ can be estimated to be $`T_c0.5`$ \[cf. Table I\]. Consequently, the case considered in Fig. 2 corresponds to a rather hot gas. That is why the simulation data are not too far from the HS values (except in the region $`q3`$). In order to highlight the differences that can be expected at a smaller temperature, the case $`T=0.5`$ is considered in Fig. 3. Now, the structure factor predicted by our model and by the PY-SHS solution is very different from that predicted by the TL perturbation theory, the latter being very close to the PY-HS solution. This is not surprising, since in any perturbation theory the quantities are expanded in powers of $`T^1`$ and obviously the value $`T^1=2`$ is beyond its range of applicability. On the other hand, the curves corresponding to the PY-SHS solution and our model are relatively close (except for a slight phase shift), as expected from the fact that the well width is rather small. Moreover, since the stickyness parameter has been chosen as to reproduce the correct second virial coefficient, the PY-SHS solution does a general good job at this very low density.

Now we consider the radial distribution function itself. As a representative example of a width not extremely small, we take the case $`\lambda =1.125`$, for which MC simulations are available. Figure 4 shows $`g(r)`$ for $`\lambda =1.125`$ at the thermodynamic state $`T=1`$, $`\rho =0.8`$. As seen in the figure, the original model of Ref.  exhibits a remarkable agreement with the simulation data. We also observe that our present model captures reasonably well the behavior of $`g(r)`$ at this high density, while the TL approximation predicts a too small contact value $`g(1^+)`$. The contact values are plotted in Fig. 5 as a function of $`1/T`$ for $`\lambda =1.125`$ and $`\rho =0.4`$, 0.6, and 0.8. At $`1/T=0`$ the system corresponds to HS and then our model and the TL theory reduce to Wertheim-Thiele’s solution of the PY equation, which tends to underestimate $`g(1^+)`$ at high densities. As the inverse temperature increases, the TL theory predicts a linear increase of $`g(1^+)`$ \[cf. Eq. (2)\] that clearly deviates from the MC data. On the other hand, the temperature dependence of $`g(1^+)`$ is well described by both the original and the simplified versions of the model, except at $`1/T=2`$.

### B Thermodynamic properties

The compressibility equation of state is obtained from Eqs. (6) and (A12). This route is preferable to the virial route because the latter is known to yield an unphysical critical behavior in the SHS limit. From Eq. (6), we have

$$Z(\eta ,T)\frac{p}{\rho k_BT}=\eta ^1_0^\eta 𝑑\eta ^{}\chi _T^1(\eta ^{},T),$$
(33)
where the density dependence of $`\chi _T`$ is given in our model by Eqs. (A12) and (29). Although this dependence is known explicitly, it does not allow us to perform the integration in Eq. (33) analytically, so that the compressibility factor $`Z`$ is obtained by numerical integration.

Before delving into the thermodynamic predictions of our model, let us compare it with the original one  in the case of a moderately narrow well, namely $`\lambda =1.125`$. Figure 6 shows the density dependence of the inverse susceptibility for $`T=0.5`$, 0.67, and 0.9. As the temperature increases, our simplified model is seen to overestimate the compressibility of the fluid, especially for large densities. On the other hand, at $`T=0.5`$ (which is practically the critical temperature, cf. Table I), both versions of the model yield undistinguishable results. This is quite encouraging, since it is in the domain of low temperatures (or, equivalently, high stickyness) where our model is expected to correct the deficiencies of perturbation theories and become useful.

One of the advantages of our simplified model is that it allows us to derive explicit expressions for the coordinates of the critical point. The spinodal line is the locus of points where the isothermal compressibility diverges. Equation (A12) clearly shows that $`\chi _T=S(0)\mathrm{}`$ if and only if $`L_2\mathrm{}`$. Thus, according to Eq. (29), the spinodal line is given by a solution to the quadratic equation $`\gamma _0+\gamma _1\eta +\gamma _2\eta ^2=0`$. This equation has two real solutions, $`\eta _\pm (\tau )`$, only if $`\tau `$ is a smaller than a critical value $`\tau _c`$ given by

$$\tau _c^1=12\frac{3+\lambda +2\sqrt{2\lambda }}{92\lambda +\lambda ^2}.$$
(34)
The critical temperature $`T_c`$ is obtained by setting $`\tau =\tau _c`$ in Eq. (31). If $`\tau <\tau _c`$, the smallest root, $`\eta _{}(\tau )`$, does not define the vapor branch of the spinodal line because it is also a root of the numerator of Eq. (29) and so $`L_2`$ remains finite at $`\eta =\eta _{}(\tau )`$. Therefore, only the liquid branch of the spinodal line, $`\eta _+(\tau )`$, exists. It is given by

$$\eta _+(T)=\frac{(1+\lambda +\lambda ^2)\sqrt{1\frac{1}{6}(3+\lambda )\tau ^1+\frac{1}{144}(92\lambda +\lambda ^2)\tau ^2}+\frac{1}{12}(3+\lambda ^3)\tau ^11\lambda +\lambda ^2}{\lambda \left[\frac{1}{3}(2+3\lambda +\lambda ^2+\lambda ^3)\tau ^14\lambda \right]}.$$
(35)
The critical density is $`\eta _c=\eta _+(\tau _c)`$, i.e.

$$\eta _c=\frac{\frac{1}{12}(3+\lambda ^3)\tau ^11\lambda +\lambda ^2}{\lambda \left[\frac{1}{3}(2+3\lambda +\lambda ^2+\lambda ^3)\tau ^14\lambda \right]}.$$
(36)
As $`\lambda `$ increases, $`T_c`$ increases, while $`\eta _c`$ decreases. Equations (34)–(36) generalize to $`\lambda >1`$ the results predicted by the PY solution for SHS,

$$\eta _+=\frac{\sqrt{96\tau ^1+\tau ^2/2}+\tau ^13}{7\tau ^112},\tau _c^1=3(2+\sqrt{2})10.24,\eta _c=\frac{3}{2}\sqrt{2}20.1213.$$
(37)
The lack of a vapor branch of the spinodal line is a feature that our model inherits from the PY-SHS solution. Our model also has in common with the SHS limit, as well as with the solution of the PY equation for finite $`\lambda `$,  the existence of regions in the temperature–density plane, inside which the physical quantities cease to take real values. According to Eq. (29), this happens when $`1+\beta _1\eta +\beta _2\eta ^2+\beta _3\eta ^3+\beta _4\eta ^4<0`$. Let us call $`\eta _i(T)`$ ($`i=1,\mathrm{},4`$) the four roots of the quartic equation $`1+\beta _1\eta +\beta _2\eta ^2+\beta _3\eta ^3+\beta _4\eta ^4=0`$, with the convention that, whenever they are real, $`\eta _1\eta _2\eta _3\eta _4`$. It turns out that the roots $`\eta _1(T)`$ and $`\eta _2(T)`$ are real only if $`\tau `$ is smaller than a certain threshold value $`\tau _{\text{th}}`$ (or, equivalently, if $`T<T_{\text{th}}`$); they define a dome-shaped curve with an apex at a density $`\eta _{\text{th}}=\eta _1(T_{\text{th}})=\eta _2(T_{\text{th}})`$. Analogously, the other two roots, $`\eta _3(T)`$ and $`\eta _4(T)`$, are real only if $`\tau <\tau _{\text{th}}^{}`$ ($`T<T_{\text{th}}^{}`$) and define another dome-shaped curve with an apex at $`\eta _{\text{th}}^{}=\eta _3(T_{\text{th}}^{})=\eta _4(T_{\text{th}}^{})`$. Consequently, the parameter $`L_2`$ becomes complex inside the intervals $`\eta _1(T)<\eta <\eta _2(T)`$ (for $`T<T_{\text{th}}`$) and $`\eta _3(T)<\eta <\eta _4(T)`$ (for $`T<T_{\text{th}}^{}`$). However, the existence of the second region is a mathematical artifact since it affects unphysically high densities. For instance, for $`\lambda =1.125`$ we have $`T_{\text{th}}^{}=0.535147`$ and $`\eta _{\text{th}}^{}=0.78692`$, this density being larger than the one corresponding to the close packing value $`\eta _{\text{cp}}=\sqrt{2}\pi /60.740`$. In fact, in the SHS limit ($`\lambda 1`$) the second region collapses into the line $`\eta =1`$ and disappears. On the other hand, the region $`\eta _1(T)<\eta <\eta _2(T)`$ is inside the curve $`\eta =\eta _\pm (T)`$ and persists in the SHS limit. In the latter case, the threshold values coincide with the critical ones, i.e. $`\tau _{\text{th}}=\tau _c`$, $`\eta _{\text{th}}=\eta _c`$. For finite $`\lambda `$, $`T_{\text{th}}`$ is smaller but practically indistinguishable from $`T_c`$ and $`\eta _{\text{th}}`$ is slightly smaller than $`\eta _c`$. In the illustrative case of $`\lambda =1.125`$, we have $`T_c=0.502324`$, $`\eta _c=0.11299`$, $`T_{\text{th}}=0.502269`$, and $`\eta _{\text{th}}=0.10980`$. The existence line, $`\rho =\rho _{1,2}(T)`$, and the liquid branch of the spinodal line, $`\rho =\rho _+(T)`$, are shown in Fig. 7 for $`\lambda =1.125`$.

Recently, it has been suggested  that, while the critical point is in general very sensitive to the range of the interaction, the second virial coefficient has a fairly constant value at the critical temperature, $`B_2(T_c)\pi `$. In the case of the SW interaction, computer simulations  show that the value of the second virial coefficient $`B_2(T)=2\pi [x(\lambda ^31)1]/3`$ at $`T=T_c`$ is much less sensitive to the width $`\lambda `$ than the critical temperature. Moreover, the PY-SHS solution predicts $`B_2(T_c)=(2+3\sqrt{2})\pi /61.04\pi `$. Thus, by assuming that $`B_2(T_c)`$ is independent of $`\lambda `$ and is equal to its value in the SHS limit, the criterion of Ref.  allows one to estimate the critical temperature as

$$T_c=1/\mathrm{ln}[1+3(2+\sqrt{2})/4(\lambda ^31)].$$
(38)
Since the model presented in this paper is constructed as a simple generalization of Baxter’s solution of the PY equation for SHS, we expect its predictions for the critical values of the temperature ($`T_c`$), the density ($`\rho _c`$), and the compressibility factor ($`Z_c`$) to be more reliable for small values of $`\lambda 1`$ than for wide wells. However, to the best of our knowledge, the simulation estimates for those quantities are only available for $`\lambda 1.25`$,  and so we cannot make a comparison for smaller widths. Table I shows the values of $`T_c`$, $`\rho _c`$, and $`Z_c`$ for several values of $`\lambda `$, as estimated from computer simulations, and as predicted by the model, by Eq. (38), by perturbation theory, and by the PY integral equation.

Comparison with MC simulation data shows that the critical temperature predicted by our model, Eq. (34), tends to be smaller than the correct value, while Vliegenthart and Lekkerkerker’s criterion,  as well as second-order perturbation theory tend to overestimate it. What is indeed remarkable is the fact that Eq. (34) provides the best agreement with computer simulations for $`\lambda 1.75`$ (except at $`\lambda =1.25`$, in which case Eq. (38) is better). In the case of the critical value of the compressibility factor, the best agreement corresponds to our model, except at $`\lambda =1.75`$, where Chang and Sandler’s perturbation theory gives a better result. On the other hand, the critical densities predicted by our model are around 30–45% smaller than the simulation values, a discrepancy that can be traced back to the solution of the PY equation for SHS. Since the computer simulations and all the theories share the property that $`\rho _c`$ is a monotonically decreasing function of $`\lambda `$, we can conclude that the correct value of $`\rho _c`$ in the limit $`\lambda 1`$ is certainly larger than the simulation value $`\rho _c=0.370`$ at $`\lambda =1.25`$, while Baxter’s solution yields $`\rho _c=0.232`$. Thus, we can expect this failure to reproduce accurately the critical density is also present in the PY approximation for finite $`\lambda `$. This is confirmed by the results of a numerical solution of the PY equation for $`\lambda =1.5`$, which gives a value of $`\rho _c`$ rather close to the one obtained here.

Now let us compare the general density dependence of the compressibility factor predicted by the model with available computer simulations. We start with the smallest value of $`\lambda `$ that, to our knowledge, has been analyzed in simulations, namely $`\lambda =1.125`$. Figure 8 shows $`Z(\rho )`$ for $`\lambda =1.125`$ and $`T=0.5`$, 0.67, and 1. Strictly speaking, the curve representing the model at the lowest temperature corresponds to $`T=T_c0.502`$ rather than to $`T=0.5`$. This is because at $`T=0.5`$ there exists a small density interval around $`\rho _c0.22`$ \[cf. Fig. 7\] where the model gives complex values. The next width we consider is $`\lambda =1.4`$, for which extensive molecular dynamics simulations are available. The results for $`T=1.25`$, 1.43, 2, and 5 are plotted in Fig. 9. We observe that, except at the highest temperatures ($`T=1`$ for $`\lambda =1.125`$, $`T2`$ for $`\lambda =1.4`$), our model presents a better general agreement with the simulation data than the TL perturbation theory.

## IV Concluding remarks

The main objective of this paper has been to propose an analytical model that could be useful to describe the structural and thermodynamic properties of systems, such as colloidal dispersions, composed of particles effectively interacting through a hard-core potential with a short-range attraction. As the simplest interaction that captures both features we have considered the square-well (SW) potential (1). On the one hand, perturbation theory becomes unreliable when the range $`\lambda 1`$ of the attraction is small, especially at low temperatures, i.e. when a certain degree of “stickyness” among the particles becomes important. On the other hand, the widely used sticky-hard-sphere (SHS) interaction model combines temperature and well width in one parameter only, thus lacking the flexibility to accommodate additional changes in width and/or temperature. Our approach intends to fill the gap between these two theories.

The model presented in this paper is based on the one proposed in Ref. , where the functions $`R(t)`$ and $`\overline{R}(t)`$ defined by Eqs. (7) and (17) were approximated by rational forms, Eq. (18). The parameters in these functions are constrained to yield a finite value for the isothermal compressibility by Eqs. (21)–(24). This still leaves two free parameters, $`A`$ and $`\tau `$, the latter being defined in Eq. (27), as unknown functions of the packing fraction $`\eta `$, the temperature parameter $`xe^{1/T^{}}1`$, and the well width $`\lambda `$. Apart from their zero-density limits (25) and (28), two extra conditions are needed to fix those parameters and close the construction of the model. Thus, depending on the physical situation one is interested in and/or on the degree of simplicity one wants to keep in the model, it is possible to choose different criteria to determine $`A`$ and $`\tau `$. For instance, one could impose certain continuity conditions on the cavity function $`y(r)g(r)e^{\phi (r)/k_BT}`$ at the points where the potential is singular; alternatively, one could require thermodynamic consistency among the virial, compressibility, and energy routes. Of course, other choices are possible. In the original formulation of the model, $`A`$ was assumed to be independent of density (hence $`A=x`$) and the second condition was the continuity of $`y(r)`$ at $`r=\lambda \sigma `$, this giving rise to a transcendent equation that needed to be solved numerically. On the other hand, in this paper we have simply assumed that both $`A`$ and $`\tau `$ are independent of density, Eq. (30). The assumption for $`\tau `$ is expected to be especially adequate for narrow potentials, since in the SHS limit the role of $`\tau `$ is played by the parameter $`\tau _{\text{SHS}}=[12x(\lambda 1)]^1`$, which is indeed independent of density. With these choices for $`A`$ and $`\tau `$ the problem remains fully algebraic and all the parameters can be expressed in terms of the solution of a quadratic equation, in analogy with what happens in Baxter’s solution of the PY equation for SHS. In fact, the model includes such a solution as a limit case. Given the scarcity of simulation results for narrow wells, we have been forced to carry out a comparison for cases with $`\lambda 1.125`$. In spite of that, the results show a very good general performance of the structural and thermodynamic properties predicted by the model, correcting the inadequacy of the perturbation theory predictions in the low-temperature domain. It is interesting to remark that the model provides an explicit expression of the critical temperature as a function of the well width, Eq. (34), which is accurate even for rather wide wells.

The results of this paper could also be useful in connection with a recently proposed extension of the law of corresponding states for systems, such as colloidal suspensions, that have widely different range of attractive interactions. Given an interaction potential $`\phi (r)=\phi _{\text{rep}}(r)+\phi _{\text{att}}(r)`$, where $`\phi _{\text{rep}}`$ and $`\phi _{\text{att}}`$ are the repulsive (not necessarily hard-core) and attractive parts, respectively, one can define a (temperature-dependent) effective hard-core diameter

$$\sigma =_0^{\mathrm{}}𝑑r\left[1e^{\phi _{\text{rep}}(r)/k_BT}\right],$$
(39)
an effective well depth

$$ϵ=\phi (r_0),\frac{d\phi (r)}{dr}|_{r=r_0}=0,$$
(40)
and a (temperature-dependent) effective well width

$$\lambda =\left[1+\left(B_2^{}1\right)\left(1e^{ϵ/k_BT}\right)^1\right]^{1/3},$$
(41)
where

$$B_2^{}=\frac{3}{\sigma ^3}_0^{\mathrm{}}𝑑rr^2\left[1e^{\phi (r)/k_BT}\right]$$
(42)
is the reduced second virial coefficient. Then, according to the extended law of corresponding states, the compressibility factor for a wide range of colloidal materials is a function of only the reduced temperature $`T^{}=k_BT/ϵ`$, the reduced density $`\rho ^{}=\rho \sigma ^3`$, and the range parameter $`\lambda `$, i.e. $`Z=F(T^{},\rho ^{},\lambda )`$, where the function $`F`$ is hardly sensitive to the details of the potential. As a consequence, an accurate prescription for the function $`F`$ based on the SW interaction for variable width can be used to determine the thermodynamic properties of a wide class of colloidal suspensions. As an illustrative example, Fig. 10 shows the dependence of the reduced critical temperature $`T_c^{}=k_BT_c/ϵ`$ on the effective width $`\lambda `$ for several interaction potentials. The prediction of our model, Eqs. (31) and (34), is also plotted. The extended law of corresponding states work very well up to $`\lambda 1.3`$. For larger interaction ranges the values of $`T_c^{}`$ for the generalized Lennard-Jones and the hard-core Yukawa potentials tend to lie slightly above those corresponding to the SW interaction.

In the near future, we plan to extend the model presented in this paper in two directions. First, the case of a mixture of particles interacting through SW potentials with different values of $`\sigma `$, $`ϵ`$, and $`\lambda `$ will be analyzed. This generalization must be such that one recovers the cases of a mixture of hard spheres,  a mixture of sticky hard spheres,  and a monodisperse SW system in the appropriate limits. As a second extension, we will study the case of an interaction model made of a hard core plus a square shoulder plus a square well, for which, in addition to the conventional gas-liquid phase transition, a liquid-liquid transition in the supercooled phase appears.

###### Acknowledgements.

Partial support from the Ministerio de Ciencia y Tecnología (Spain) through grant No. BFM2001-0718 and from the Junta de Extremadura (Fondo Social Europeo) through grant No. IPR99C031 is gratefully acknowledged. A.S. is also grateful to the DGES (Spain) for a sabbatical grant No. PR2000-0117.

## A Some explicit expressions

### 1 Expression for $`S(0)`$

In this Appendix we list some of the expressions that are derived from the model. The long-wavelength value of the structure factor is obtained by using the model (20) in the general expression (14). The result is

$`S(0)`$ $`=`$ $`{\displaystyle \frac{1}{5(1+2\eta )^2}}\{520\eta [1+(2+\lambda )A^{}3(1+\lambda )L_2]`$ (A12)
$`+2\eta ^2[156(14\lambda +19\lambda ^2\lambda ^3\lambda ^4)L_2+120(1+\lambda +\lambda ^2)L_2^2`$
$`+\left(50+16\lambda +27\lambda ^22\lambda ^3\lambda ^4\right)A^{}30(5+4\lambda +3\lambda ^2)A^{}L_2`$
$`+15(3+2\lambda +\lambda ^2)A_{}^{}{}_{}{}^{2}]`$
$`2\eta ^3[10+3(753\lambda 13\lambda ^2+7\lambda ^38\lambda ^4)L_2`$
$`60(1+\lambda )(14\lambda +\lambda ^2)(1+\lambda +\lambda ^2)L_2^2`$
$`+\left(19+59\lambda +9\lambda ^2\lambda ^3+4\lambda ^4\right)A^{}`$
$`+3(1167\lambda 114\lambda ^266\lambda ^313\lambda ^4+9\lambda ^5)A^{}L_2`$
$`+3(3+2\lambda +\lambda ^2)(1+7\lambda +3\lambda ^2\lambda ^3)A_{}^{}{}_{}{}^{2}]`$
$`+\eta ^4[512(1+\lambda +\lambda ^2+11\lambda ^34\lambda ^4)L_2+240\lambda ^3(1+\lambda +\lambda ^2)L_2^2`$
$`+2\left(7+8\lambda +9\lambda ^2+10\lambda ^34\lambda ^4\right)A^{}12(1+3\lambda +6\lambda ^2+24\lambda ^3+17\lambda ^4+9\lambda ^5)A^{}L_2`$
$`+3(3+2\lambda +\lambda ^2)(1+2\lambda +3\lambda ^2+4\lambda ^3)A_{}^{}{}_{}{}^{2}]\},`$
where we have made use of Eqs. (21)–(24).

### 2 Expressions for the coefficients in Eq. (29)

Equation (27) reduces to a quadratic equation for $`L_2`$. Its physical solution is given by Eq. (29), where

$$\alpha _1=2A^{}(2+\lambda )+\frac{1}{6}\left(1+4\lambda +\lambda ^23A^{}\right)\tau ^1,$$
(A13)
$$\alpha _2=3+A^{}(7+2\lambda 3\lambda ^2)\frac{1}{12}\left[7+\lambda +16\lambda ^2+A^{}(23+15\lambda +15\lambda ^2+7\lambda ^3)\right]\tau ^1,$$
(A14)
$$\alpha _3=22A^{}(1+2\lambda +3\lambda ^2)+\frac{1}{6}\left[7+\lambda 2\lambda ^2+A^{}(7+15\lambda +21\lambda ^2+11\lambda ^3+6\lambda ^4)\right]\tau ^1$$
(A15)
$$\beta _1=44A^{}(2+\lambda )+\frac{1}{3}(5+2\lambda \lambda ^2+3A^{})\tau ^1\frac{1}{3}\tau ^2,$$
(A16)
$`\beta _2`$ $`=`$ $`6+6A^{}(3+2\lambda +\lambda ^2)+4A_{}^{}{}_{}{}^{2}(2+\lambda )^2`$ (A19)
$`{\displaystyle \frac{1}{6}}\left[9(3+\lambda )+A^{}(59+51\lambda +3\lambda ^25\lambda ^3)+12A_{}^{}{}_{}{}^{2}(2+\lambda )\right]\tau ^1`$
$`+{\displaystyle \frac{1}{36}}\left[31+8\lambda +18\lambda ^24\lambda ^3+\lambda ^4+12A^{}(5+\lambda )+9A_{}^{}{}_{}{}^{2}\right]\tau ^2,`$
$`\beta _3`$ $`=`$ $`412A^{}(1+\lambda +\lambda ^2)4A_{}^{}{}_{}{}^{2}(2+\lambda )(1+2\lambda +3\lambda ^2)`$ (A24)
$`+{\displaystyle \frac{1}{3}}[3(4+\lambda +\lambda ^2)+A^{}(29+36\lambda +30\lambda ^2+7\lambda ^33\lambda ^4)`$
$`+A_{}^{}{}_{}{}^{2}(17+43\lambda +30\lambda ^2+\lambda ^3\lambda ^4)]\tau ^1`$
$`{\displaystyle \frac{1}{36}}[29+13\lambda +27\lambda ^2+7\lambda ^34\lambda ^4+A^{}(68+80\lambda +80\lambda ^2+16\lambda ^3+7\lambda ^4+\lambda ^5)`$
$`+3A_{}^{}{}_{}{}^{2}(13+13\lambda +\lambda ^23\lambda ^3)]\tau ^2,`$
$`\beta _4`$ $`=`$ $`[1+A^{}(1+2\lambda +3\lambda ^2)]^2{\displaystyle \frac{1}{6}}[7+\lambda +4\lambda ^2`$ (A28)
$`+2A^{}(7+15\lambda +15\lambda ^2+5\lambda ^3+6\lambda ^4)+A_{}^{}{}_{}{}^{2}(1+2\lambda +3\lambda ^2)(7+15\lambda +3\lambda ^2\lambda ^3)]\tau ^1`$
$`+{\displaystyle \frac{1}{144}}[49+38\lambda +9\lambda ^2+8\lambda ^3+16\lambda ^4+2A^{}(49+88\lambda +112\lambda ^2+80\lambda ^3+35\lambda ^44\lambda ^5)`$
$`+A_{}^{}{}_{}{}^{2}(49+138\lambda +219\lambda ^2+124\lambda ^3+27\lambda ^4+18\lambda ^5+\lambda ^6)]\tau ^2,`$
$$\gamma _0=1+\lambda \frac{1}{6}\tau ^1,$$
(A29)
$$\gamma _1=2(1+\lambda \lambda ^2)\frac{1}{6}(3+\lambda ^3)\tau ^1,$$
(A30)
$$\gamma _2=\lambda \left[4\lambda +\frac{1}{3}(2+3\lambda +\lambda ^2+\lambda ^3)\tau ^1\right].$$
(A31)
The other root is incompatible with (25) and then must be discarded.

## B Low-density behavior of the model

In the low-density limit, Eqs. (21)-(24) become, respectively,

$$L_1=1\frac{3}{2}\eta [1x(\lambda ^41)]+𝒪(\eta ^2),$$
(B1)
$$S_1=𝒪(\eta ),$$
(B2)
$$S_2=\frac{1}{2}[1x(\lambda ^21)]+𝒪(\eta ),$$
(B3)
$$S_3=\frac{1}{12\eta }+\frac{1}{3}[1x(\lambda ^31)]+𝒪(\eta ).$$
(B4)
To first order in density, Eq. (29) gives

$$L_2=x\lambda (\lambda 1)+L_2^{(1)}\eta +𝒪(\eta ^2),$$
(B5)
where

$$L_2^{(1)}=\frac{3}{2}x\lambda (\lambda 1)\left\{1+x(\lambda 1)^2\left[12\lambda \lambda ^24x\lambda (1+\lambda )\right]\right\}.$$
(B6)
Substitution into Eq. (20) yields, after some algebra,

$`F(t)`$ $`=`$ $`F_{\text{exact}}(t)+\{{\displaystyle \frac{2C_2}{t^4}}+{\displaystyle \frac{C_1}{t^3}}+{\displaystyle \frac{C_0}{t^2}}[{\displaystyle \frac{2C_2}{t^4}}+{\displaystyle \frac{C_1+2C_2(\lambda 1)}{t^3}}`$ (B8)
$`+{\displaystyle \frac{C_0+C_1(\lambda 1)+C_2(\lambda 1)^2}{t^2}}]e^{(\lambda 1)t}\}\eta +𝒪(\eta ^2),`$
where we have called

$$C_23x(1+x)(\lambda ^21),$$
(B9)
$$C_12x(1+x)(\lambda 1)^2(1+2\lambda ),$$
(B10)
$$C_0\frac{L_2^{(1)}}{\lambda 1}\frac{1}{2}x(28\lambda ^3+3\lambda ^4)2x^2(\lambda 1)^2(2+4\lambda +3\lambda ^2)3x^3(\lambda ^21)^2.$$
(B11)
The expression for the exact radial distribution function, $`g_{\text{exact}}(r)`$, to first order in density was derived by Barker and Henderson for the case $`\lambda <2`$. The corresponding expression for $`F_{\text{exact}}(t)`$ can be found in Ref. . From Eq. (B8) one easily gets

$$g(r)g_{\text{exact}}(r)=\frac{1}{r}[C_2(r1)^2+C_1(r1)+C_0][\mathrm{\Theta }(r1)\mathrm{\Theta }(r\lambda )]\eta +𝒪(\eta ^2).$$
(B12)
Thus, the difference (to first order in density) is non-zero in the interval $`1>r>\lambda `$ only. In particular,

$$g(1^+)g_{\text{exact}}(1^+)=C_0\eta +𝒪(\eta ^2),$$
(B13)
where

$`g_{\text{exact}}(1^+)`$ $`=`$ $`1+x+[{\displaystyle \frac{5}{2}}+{\displaystyle \frac{x}{2}}(1516\lambda ^3+6\lambda ^4)2x^2(\lambda 1)(4+4\lambda +\lambda ^23\lambda ^3)`$ (B15)
$`+3x^3(\lambda ^21)^2]\eta +𝒪(\eta ^2).`$
Note that the relative coefficient $`C_0/(1+x)`$ vanishes in the SHS limit. As said in Sec. II, the model presented in this paper does not enforce the verification of Eq. (26). In fact, Eq. (B12) implies that

$`g(\lambda ^{})(1+x)g(\lambda ^+)`$ $`=`$ $`\lambda ^1\left[C_2(\lambda 1)^2+C_1(\lambda 1)+C_0\right]\eta +𝒪(\eta ^2)`$ (B16)
$`=`$ $`{\displaystyle \frac{x(\lambda 1)}{2\lambda }}[\lambda (11\lambda \lambda ^2)x(\lambda 1)(\lambda +2)(3+10\lambda 3\lambda ^2)`$ (B18)
$`+6x^2(\lambda 1)^2(1+3\lambda +2\lambda ^2)]\eta +𝒪(\eta ^2).`$
Finally, from Eq. (A12) or, equivalently, inserting Eq. (B12) into Eq. (3), we get

$`S(0)S_{\text{exact}}(0)`$ $`=`$ $`2x(\lambda 1)^2[7+23\lambda +30\lambda ^24\lambda ^32\lambda ^4x(\lambda 1)`$ (B20)
$`\times (25+84\lambda +78\lambda ^2+2\lambda ^39\lambda ^4)+18x^2(\lambda ^21)^2(1+2\lambda )]\eta ^2+𝒪(\eta ^3),`$
where

$`S_{\text{exact}}(0)`$ $`=`$ $`18[1x(\lambda ^31)]\eta +2[17x(\lambda 1)(19+19\lambda +19\lambda ^2+51\lambda ^33\lambda ^43\lambda ^5)`$ (B22)
$`2x^2(\lambda 1)^2(8+16\lambda 3\lambda ^238\lambda ^319\lambda ^4)+18x^3(\lambda ^21)^3]\eta ^2+𝒪(\eta ^3).`$
## C The hard-sphere and sticky-hard-sphere limits

### 1 Hard spheres

The SW potential becomes equivalent to the HS potential if $`\lambda =1`$ at any non-zero temperature $`T`$ or if $`T\mathrm{}`$ at any width $`\lambda `$. Let us first consider the latter limit. Making $`x=0`$ in Eqs. (A13)–(A31), one gets

$$\alpha _1=0,\alpha _2=3,\alpha _3=2,$$
(C1)
$$\beta _1=4,\beta _2=6,\beta _3=4,\beta _4=1,$$
(C2)
$$\gamma _0=1+\lambda ,\gamma _1=2(1+\lambda \lambda ^2),\gamma _3=4\lambda ^2.$$
(C3)
As a consequence, Eq. (29) implies that $`L_2=0`$. Thus Eq. (20) becomes

$$F(t)=\frac{1}{12\eta }\frac{1+L_1t}{1+S_1t+S_2t^2+S_3t^3},$$
(C4)
with

$$L_1=\frac{1+\frac{1}{2}\eta }{1+2\eta },$$
(C5)
$$S_1=\frac{3}{2}\frac{\eta }{1+2\eta },$$
(C6)
$$S_2=\frac{1}{2}\frac{1\eta }{1+2\eta }$$
(C7)
$$S_3=\frac{(1\eta )^2}{12\eta (1+2\eta )}.$$
(C8)
This is precisely the form adopted by $`F(t)`$ in the analytical solution of the PY equation for hard spheres. The expression for $`S(0)`$, Eq. (A12), simply reduces to

$$S(0)=\frac{(1\eta )^4}{(1+2\eta )^2}.$$
(C9)
The case $`\lambda =1`$ is not considered in this Subsection, as it is a particular case of the SHS limit.

### 2 Sticky hard spheres

Let us now take the limit $`x\mathrm{}`$, $`\lambda 1`$, with $`x(\lambda 1)=\text{finite}`$ in the model proposed in this paper. In that limit the parameter $`\tau `$ is finite, cf. Eq. (30), while $`A^{}=0`$. Equations (A13)–(A31) become

$$\alpha _1=\tau ^1,\alpha _2=32\tau ^1,\alpha _3=2+2\tau ^1,$$
(C10)
$$\begin{array}{c}\beta _1=4+2\tau ^1\frac{1}{3}\tau ^2,\beta _2=66\tau ^1+\frac{3}{2}\tau ^2,\\ \beta _3=4+6\tau ^12\tau ^2,\beta _4=12\tau ^1+\frac{5}{6}\tau ^2,\end{array}$$
(C11)
$$\gamma _0=2\frac{1}{6}\tau ^1,\gamma _1=2\frac{2}{3}\tau ^1,\gamma _3=4+\frac{7}{3}\tau ^1.$$
(C12)
Therefore, Eq. (29) reduces to

$$L_2=\frac{1\eta }{24\eta }\frac{(1+2\eta )\left[(1\eta )^2+2\eta (1\eta )\tau ^1\frac{1}{6}\eta (25\eta )\tau ^2\right]^{1/2}(1\eta )(1+2\eta \eta \tau ^1)}{(1\eta )(1+2\eta )\frac{1}{12}(1+4\eta 14\eta ^2)\tau ^1}.$$
(C13)
Taking the limit $`\lambda 1`$ in Eq. (20) we get

$$F(t)=\frac{1}{12\eta }\frac{1+L_1t+L_2t^2}{1+S_1t+S_2t^2+S_3t^3},$$
(C14)
with

$$L_1=\frac{1+\frac{1}{2}\eta }{1+2\eta }+\frac{6\eta }{1+2\eta }L_2,$$
(C15)
$$S_1=\frac{3}{2}\frac{\eta }{1+2\eta }+\frac{6\eta }{1+2\eta }L_2,$$
(C16)
$$S_2=\frac{1}{2}\frac{1\eta }{1+2\eta }+\frac{14\eta }{1+2\eta }L_2$$
(C17)
$$S_3=\frac{(1\eta )^2}{12\eta (1+2\eta )}\frac{1\eta }{1+2\eta }L_2.$$
(C18)
This coincides with the analytical solution of the PY equation for sticky hard spheres. From Eq. (A12) we have

$$S(0)=\frac{(1\eta )^2}{(1+2\eta )^2}\left[1\eta +12\eta L_2\right]^2.$$
(C19)
Of course, the results for hard spheres are recovered in the high-temperature limit ($`\tau ^10`$).


format: text

text:
# Optimal quantum estimation of the coupling between two bosonic modes

## 1 Introduction

How effectively may we estimate the strength of a simple interaction such as the transmittivity of a beam splitter or the phase-shift imposed in the internal arms of an interferometer ? The Hamiltonian describing the bilinear coupling between two bosonic modes has the form

$$H=\kappa (a^{}b+b^{}a),$$
(1)
where $`\kappa `$ depends on the specific interaction under consideration. By using the Schwinger representation of the $`SU(2)`$ Lie algebra, with generators

$$J_x=\frac{1}{2}(a^{}b+b^{}a),J_y=\frac{1}{2i}(a^{}bb^{}a),J_z=\frac{1}{2}(a^{}ab^{}b),$$
(2)
we can rewrite the Hamiltonian as $`H=2\kappa J_x`$ and the evolution operator as

$$U_\psi =\mathrm{exp}(\mathrm{i}J_x\psi ),$$
(3)
where the global coupling constant $`\psi `$ is equal to $`2\kappa \mathrm{\Delta }t`$, $`\mathrm{\Delta }t`$ being the effective interaction time. The evolution in Eq. (3) describes, for example, the interaction of two light modes in a beam splitter with transmittivity $`\tau =\mathrm{cos}^2\psi `$ , or, apart from a fixed rotation, the evolution of the arm modes in a Mach-Zehnder interferometer, with $`\psi `$ representing the phase-shift between arms . If the initial preparation of two modes is described by the density matrix $`\rho _0`$ the evolved state in the interaction picture is given by

$$\rho _\psi =\mathrm{exp}(\mathrm{i}J_x\psi )\rho _0\mathrm{exp}(\mathrm{i}J_x\psi ).$$
(4)
In this paper we devote our attention to the estimation of $`\psi `$ through measurements performed on $`\rho _\psi `$. We denote the generic POVM for such estimation process by $`\mathrm{d}\nu (\varphi )`$, such that the outcomes of the measurement are distributed according to

$$p(\varphi |\psi )\mathrm{d}\varphi =\mathrm{Tr}[\rho _\psi \mathrm{d}\nu (\varphi )],$$
(5)
where $`p(\varphi |\psi )`$ represents the conditional probability of registering the outcome $`\varphi `$ when the true value of the parameter is $`\psi `$.

Our objective is to find the best strategy, i.e. the POVM that provides the optimal estimation of the parameter $`\psi `$ . Since $`\psi `$ is manifestly a phase-shift, we can use general results from phase estimation theory, which provides the optimal POVM to estimate the phase-shift induced by a phase generator, i.e. a selfadjoint operator with discrete, equally spaced, spectrum. The optimality criterion is given in terms of the minimization of the mean value of a cost function  that assesses the quality of the strategy, i.e. it weights the errors in the estimates. Since a phase-shift is a $`2\pi `$-periodic parameter, the cost function must be a $`2\pi `$-periodic even function of $`(\psi \varphi )`$, i.e. it has Fourier expansion with cosines only. The appropriate concavity properties of the cost function are reflected by expansion coefficients which are all negative apart the irrelevant additive constant. A cost function with such Fourier expansion has been firstly considered by Holevo , and for such reason it is usually referred to as belonging to Holevo’s class. Notice that the optimal POVM for a given state $`\rho _0`$ is the same for every cost function in this class. A more general quantum estimation approach for different kinds of phase-shift and general quantum system is given in Ref. . There, it is also shown that the kind of problem we are presently dealing with is, in a certain sense, the best situation, since the spectrum of our phase generator is the whole set of integers $`\mathrm{}`$, including negative ones. In this case, there is an optimal orthogonal projective POVM, which can be regarded as the spectral resolution of a self-adjoint phase operator. However, if the estimation is performed with the constraint of bounded or fixed energy, the optimal POVM and the optimal input state $`\rho _0`$ do not correspond to a canonical quantum observable scheme . Moreover, in general the optimal POVM depends on the input state, and the optimal POVM of Ref.  holds only for pure states, whereas a generalization to a class of mixed states, the so-called phase-pure states, has been considered in Ref. .

The optimal POVM, in the sense described above, provides an unbiased estimation of $`\psi `$ for preparation $`\rho _0`$ of the two modes, in formula

$$\varphi =\psi \mathrm{with}\varphi =_0^{2\pi }\varphi \mathrm{Tr}[\rho _\psi \mathrm{d}\nu (\varphi )]\rho _\psi =U_\psi \rho _0U_\psi ^{}.$$

On the other hand, we want also to find the optimal state $`\rho _0`$ for the estimation of $`\psi `$ according to some cost function, which quantifies the noise of the estimation. The customary root mean square is not a good choice for a cost function, since the function $`(\varphi \psi )^2`$ is not $`2\pi `$-periodic. A good definition for the precision of the measurement is given by the average of the cost function $`C(\varphi \psi )=4\mathrm{sin}^2(\frac{\varphi \psi }{2})`$, i. e. a ”periodicized variance”, which obviously belongs to the Holevo’s class. If the estimates occur within a small interval around the true value of the parameter $`\psi `$, one has approximately $`C(\varphi \psi )(\varphi \psi )^2`$, whence $`\delta \psi =\sqrt{C}`$ can be assumed as a reasonable measure of the precision of the measurement.

## 2 Optimal estimation of the coupling parameter

In order to solve our estimation problem, let us consider the following unitary transformation

$$𝒰=\mathrm{exp}\left\{\frac{\pi }{4}(a^{}bb^{}a)\right\}=\mathrm{exp}\left\{\mathrm{i}\frac{\pi }{2}J_y\right\}.$$
(6)
Using Eq. (6) we may rewrite Eq. (5) in the more familiar form of rotation along the z-axis

$`p(\varphi |\psi )\mathrm{d}\varphi `$ $`=`$ $`\mathrm{Tr}\left[𝒰\rho _\psi 𝒰^{}𝒰\mathrm{d}\nu (\varphi )𝒰^{}\right]`$ (7)
$`=`$ $`\mathrm{Tr}\left[\mathrm{exp}(\mathrm{i}J_z\psi )𝒰\rho _0𝒰^{}\mathrm{exp}(\mathrm{i}J_z\psi )𝒰\mathrm{d}\nu (\varphi )𝒰^{}\right]`$
$`=`$ $`\mathrm{Tr}\left[\mathrm{exp}(\mathrm{i}J_z\psi )R_0\mathrm{exp}(\mathrm{i}J_z\psi )\mathrm{d}\mu (\varphi )\right],`$
where we used the identity $`𝒰J_x𝒰^{}=J_z`$. Equation (7) shows that the problem of estimating the shift generated by $`J_x`$ on the state $`\rho _0`$ is equivalent to that of estimating the same shift generated by $`J_z`$ on the rotated state $`R_0=𝒰\rho _0𝒰^{}`$. In particular, any POVM $`\mathrm{d}\nu (\varphi )`$ to estimate the $`J_x`$-induced shift can be written as $`d\nu (\varphi )=𝒰^{}\mathrm{d}\mu (\varphi )𝒰`$, where $`\mathrm{d}\mu (\varphi )`$ is a POVM for the $`J_z`$-induced shift estimation.

For pure states $`R_0=|\psi _0\psi _0|`$ (in the following we use double brackets for two-mode vectors) the degeneration of the spectrum of $`a^{}ab^{}b`$ can be treated using the technique introduced in , and the optimal POVM for cost functions in Holevo’s class  is proved to be of the form

$$\mathrm{d}\mu (\varphi )=\frac{\mathrm{d}\varphi }{2\pi }|E_\varphi E_\varphi |$$
(8)
with the vectors $`|E_\varphi `$ given by

$$|E_\varphi =\underset{d}{}\mathrm{e}^{id\varphi }|d.$$
(9)
The vectors $`|d`$ are certain eigenvectors of $`D=a^{}ab^{}b`$ built by picking up, in every eigenspace $`_d`$ of the eigenvalue $`d`$, the normalized vector parallel to the projection of $`|\psi _0`$ on $`_d`$. In order to be more specific, let us consider an input state of the form

$$|\psi _0=\underset{n=0}{\overset{\mathrm{}}{}}\psi _n^{(0)}|n,0+\underset{n=0}{\overset{\mathrm{}}{}}\underset{d=1}{\overset{\mathrm{}}{}}[\psi _n^{(d)}|n,d+\psi _n^{(d)}|n,d],$$
(10)
where $`|n,d`$ is given by

$`|n,d\{\begin{array}{cc}|n+d_a|n_b& \mathrm{if}d0\\ |n_a|nd_b& \mathrm{if}d<0\end{array}.`$ (13)
The projection of $`|\psi _0`$ in $`_d`$ is equal to

$$\underset{n=0}{\overset{\mathrm{}}{}}\psi _n^{(d)}|n,d,$$
(14)
such that the eigenvector $`|d`$ reads as follows

$$|d=\frac{_{n=0}^{\mathrm{}}\psi _n^{(d)}|n,d}{\sqrt{_{n=0}^{\mathrm{}}\left|\psi _n^{(d)}\right|^2}},$$
(15)
and the input state can be rewritten as

$$|\psi _0=\underset{d}{}\gamma _d|d\gamma _d=\sqrt{\underset{n=0}{\overset{\mathrm{}}{}}\left|\psi _n^{(d)}\right|^2}.$$

Notice that the dependence of the POVM on the state $`|\psi _0`$ is contained in the vectors $`|E_\varphi `$.

By adopting $`C(\varphi \psi )`$ as a cost function the average cost of the strategy corresponds to the expectation value of the cost operator $`C=2\mathrm{E}_+\mathrm{E}_{}`$, where the raising and lowering operators $`E_+`$ and $`E_{}`$ are given by

$$\mathrm{E}_+=\underset{d}{}|d+1d|\mathrm{E}_{}=\mathrm{E}_+^{}$$

(with the vectors $`|d`$ defined as above).

The optimization problem is that of minimizing the average cost of the strategy

$`\overline{C}={\displaystyle _0^{2\pi }}{\displaystyle \frac{d\psi }{2\pi }}{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\varphi }{2\pi }}C(\varphi \psi )p(\varphi |\psi )=\mathrm{Tr}\left[R_0C\right]\psi _0|2\mathrm{E}_+\mathrm{E}_{}|\psi _0`$
$`,`$ (16)
with the constraint that the solution is a normalized state. The Lagrange function is given by

$$=\overline{C}\lambda \psi _0|\psi _0,$$
(17)
with $`\lambda `$ being the Lagrange multiplier for the normalization constraint. The solution of this problem is a state with infinite mean energy $`N=\psi _0|a^{}a+b^{}b|\psi _0`$. In order to find physically realizable states, one must impose a constraint on $`N`$ too, and the Lagrange function becomes

$$=\overline{C}\mu \psi _0|a^{}a+b^{}b|\psi _0\lambda \psi _0|\psi _0,$$
(18)
$`\mu `$ being the Lagrange multiplier for the mean energy. It is useful, in order to calculate the solution of this equation, to write the generic state $`|\psi _0`$ in the following way

$$|\psi _0=\underset{d}{}\psi _d\underset{n=0}{\overset{\mathrm{}}{}}c_{n,d}|n,d.$$
(19)
The coefficients $`c_{n,d}`$ determine the normalized projection of $`|\psi `$ into the eigenspace $`_d`$, whereas the $`\psi _d`$’s are the coefficients which combine those projections.

Using Eq. (19), the Lagrange function (18) explicitly shows terms accounting for the normalization of the vectors $`_{n=0}^{\mathrm{}}c_{n,d}|n,d`$ in each $`_d`$, with $`\nu ^{(d)}`$ denoting Lagrange multipliers for the normalization of projections, and rewrites as

$``$ $`=`$ $`{\displaystyle \underset{d}{}}\{2|\psi _d|^2{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}|c_{n,d}|^2(\overline{\psi }_d{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}|c_{n,d}|^2)(\psi _{d1}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}|c_{m,d1}|^2)+`$ (20)
$`\left(\overline{\psi }_d{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}|c_{n,d}|^2\right)\left(\psi _{d+1}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}|c_{m,d+1}|^2\right)+`$
$`\mu |\psi _d|^2{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}|c_{n,d}|^2(2n+|d|)\lambda |\psi _d|^2{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}|c_{n,d}|^2\nu ^{(d)}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}|c_{n,d}|^2\}.`$
By taking derivatives of the Lagrange function with respect to $`c_{n,d}^{}`$ and $`\psi _d^{}`$ with the constraints

$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}|c_{n,d}|^2=1{\displaystyle \underset{d}{}}|\psi _d|^2=1`$
$`{\displaystyle \underset{d}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}|\psi _d|^2|c_{n,d}|^2(2n+|d|)=N,`$ (21)
and by rephasing the $`|d`$’s we arrive at the system

$`\{\begin{array}{c}(2\lambda )\psi _d\psi _{d1}\psi _{d+1}\mu _{n=0}^{\mathrm{}}(2n+|d|)|c_{n,d}|^2\psi _d=0\hfill \\ \left[(2\lambda )\psi _d^22(\psi _d\psi _{d1}+\psi _{d+1}\psi _d)\mu \psi _d^2(2n+|d|)\nu ^{(d)}\right]c_{n,d}=0\hfill \end{array}`$ (24)
The second equation in (24) implies that for a fixed $`d`$ only one coefficient $`c_{n,d}`$ can be different from zero, say for the value $`\overline{n}`$, and in this case $`|c_{\overline{n},d}|=1`$ 111Indeed, let us suppose $`c_{n,d}0`$ for two values $`n_1=m`$ and $`n_2=p`$. Then we must have $`(2\lambda )\psi _d^22(\psi _d\psi _{d1}+\psi _{d+1}\psi _d)\mu \psi _d^2(2p+|d|)\nu ^{(d)}=0`$ which implies $`2\mu (mp)\psi _d=0`$. Since the case $`\mu =0`$ is not interesting and $`\psi _d=0`$ would imply that the choice of the $`c_{n,d}`$ is completely arbitrary and irrelevant, the only possibility is $`m=p`$..

The first equation of (24) can therefore be rewritten as

$$(2\lambda )\psi _d\psi _{d1}\psi _{d+1}\mu (2n(d)+|d|)\psi _d=0,$$
(25)
which allows us to obtain from the second one

$$\nu ^{(d)}=\psi _d(\psi _{d+1}+\psi _{d1}).$$
(26)
The solutions of Eq. (26) give local minima for the average cost $`\overline{C}`$, and one should solve the equation (25) for arbitrary choices of $`n(d)`$, looking for the optimal one. In the case $`n(d)=0`$ we have

$`{\displaystyle \frac{2(\lambda ^{}+|d|)}{\frac{2}{\mu ^{}}}}\psi _d=\psi _{d1}+\psi _{d+1},`$ (27)
where $`\mu ^{}=\mu `$ and $`\lambda ^{}=\frac{2\lambda }{\mu ^{}}`$. Eq. 27 is the recursion equation for Bessel functions, with solution given by

$$\psi _d=𝒩^{1/2}(\mu ^{},\lambda ^{})J_{\lambda ^{}+|d|}\left(\frac{2}{\mu ^{}}\right),$$
(28)
with $`𝒩(\mu ^{},\lambda ^{})=_dJ_{\lambda ^{}+|d|}^2\left(\frac{2}{\mu ^{}}\right)`$ and with the boundary conditions $`J_\lambda ^{}^{}(2/\mu ^{})=0`$, i.e. $`J_{\lambda ^{}+1}=J_{\lambda ^{}1}`$. Finally, to obtain the optimal state $`\rho _0`$ one has to rotate $`R_0=|\psi _0\psi _0|`$ by the unitary transformation (6).

In order to obtain the behaviour of the average cost versus the energy we numerically solved Eq. (25) with $`n(d)=0`$. This problem can be rewritten as the eigenvalue problem $`𝐀\psi =\lambda \psi `$ for the matrix $`𝐀`$ with elements given by

$$(𝐀)_{m,n}=(2\mu |m|)\delta _{m,n}\delta _{m,n+1}\delta _{m,n1}.$$
(29)
Numerical diagonalization gives the power-law $`\overline{C}\frac{\gamma }{N^2}`$ in the range $`0N1000`$ with $`\gamma 0.1`$. This behaviour is plotted in Fig. 1, where also other solutions of equation (25) corresponding to local minima of the average cost are shown. Since the phase distribution of the optimal state is singly peaked we may also write $`\delta \psi \sqrt{\overline{C}}\gamma ^{1/2}/N`$, which means that the optimal states derived here are at the so-called Heisenberg limit of phase variance.

## 3 Conclusions

In conclusion, we dealt with the problem of estimating the coupling constant of a bilinear interaction. In our approach the coupling constant, which appears in the exponent of the time evolution operator, has been treated as a phase parameter. The optimal POVM has been derived according to the theory of quantum phase estimation . As noticed in Ref.  this resorts to an SU(2) estimation problem with the Schwinger two-mode boson realization. However, the representation is not irreducible, and a more complicated problem is faced. In this sense our results generalize those of Ref. , where the estimation of $`SU(2)`$ phase-shifts has been analyzed in irreducible subspaces: indeed we found an improved scaling of phase variance versus the total energy. The degeneracy of the spectrum of the the Hamiltonian, can be treated using the technique of Ref. , and, in this way, the problem is reduced to a nondegenerate one with spectrum $``$ for the phase-shift operator.

It is important to remark that the optimal POVM depends on the preparation state. This is true for every kind of phase estimation problem, but in the presence of degeneracy the dependence is crucial. In fact, one has to define the optimal POVM as a block diagonal operator, where the invariant subspaces are spanned by projections of the input state into the eigenspaces of the generator. From a practical point of view this means that optimal estimation of the phase-shift imposed to a state, needs a measuring device which is adapted to the shifted state. We then optimized the input state by minimizing the average cost for fixed input energy and found a power law $`\delta \psi \gamma /N`$ in a range $`0N1000`$.

Notice that the law $`N^1`$ is the same as in the optimal phase estimation with only one mode , however with a much smaller constant $`\gamma `$ ($`\gamma 0.1`$ instead of $`\gamma 1.36`$). We think that this phenomenon of improvement of phase sensitivity by increasing the number of modes is the same considered in Ref. , where an exponential improvement versus the number of modes has been estimated when increasing the number of modes and the the number of photons per mode, jointly in the same proportion.

## Acknowledgement

This work has been cosponsored by MURST under the project ”Quantum Information Transmission and Processing”. The authors thank Dipartimento di Fisica ”Alessandro Volta” for partial support.