text:
# Phase separation of penetrable core mixtures

## I Introduction

Demixing of binary or multicomponent mixtures is a very common phenomenon observed in a broad range of molecular fluids ref1 , polymer solutions and blends ref2a ; ref2b , or colloidal dispersions ref3 ; ref4 . Phase separation is generally associated with differences in the attractive interactions between particles of different chemical species. In polymer solutions these differences are usually embodied in the Flory $`\chi `$-parameter ref2a , which controls the competition between the entropy of mixing and the total interaction energy, at least within a mean-field picture. On the other hand, in multi-component colloidal systems like binary dispersions involving colloidal particles of very different sizes, or mixtures of colloidal particles and non-adsorbing polymer, phase separation can be driven by purely repulsive, excluded volume interactions. By mapping the initial multi-component system onto an effective one-component system involving only the bigger colloidal particles, the largely entropy-driven demixing can be understood in terms of attractive depletion interactions induced between the large particles by the smaller species (the “depletant”) ref4 . It should however be kept in mind that the initial bare interactions are purely repulsive, albeit strongly non-additive, as in the highly simplified Asakura-Oosawa model for colloid-polymer mixtures ref5 . In the case of fully additive hard sphere mixtures, phase separation has been predicted for sufficiently large size ratios ref6 , but it is now generally believed that the fluid-fluid demixing is metastable, and preempted by freezing ref7 . A significant degree of positive non-additivity of the core radii $`R_{\mu \nu }`$ (whereby $`R_{12}=(R_{11}+R_{22})(1+\mathrm{\Delta })/2`$, with $`\mathrm{\Delta }>0`$) is required to observe a stable demixing transition in the fluid phase ref8 .

Effective interactions between the centres of mass of fractal objects, like linear polymer coils ref9 ; ref10 ; ref11 ; ref12 , star polymers ref13  or dendrimers ref14 , obtained by averaging over individual monomer degrees of freedom, are now known to be very “soft”. More specifically the effective pair potential diverges only logarithmically for overlapping star polymers ref13 , while remaining finite, of the order of $`12k_BT`$, for linear polymers in good solvent ref9 ; ref10 ; ref11 ; ref12 . This observation has stimulated the investigation of simple models, like finite repulsive step potentials ref15 , or the gaussian core potential ref16 ; ref17 ; ref18 , which was first introduced by Stillinger, in a somewhat different context ref19 , namely

$$v(r)=ϵ\mathrm{exp}\left(r^2/R^2\right),$$
(1)
where $`ϵ`$ is the energy scale, while $`R`$ determines the range of the effective potential. It is worth stressing that the “gaussian core” model is unrelated to the “gaussian molecule” model, which was extensively studied by Michael Fisher and collaborators ref26 . In the latter model it is the Mayer $`f`$-function, rather than the pair potential, which has a gaussian shape.

Simple, penetrable particle models are also widely used in highly coarse-grained simulations of large-scale phenomena within the so-called “dissipative particle dynamics” (DPD) method ref20 ; ref21 . In DPD, effective interactions between penetrable fluid “particles” are frequently modelled by a simple parabolic potential ref21 ; ref22 :

$$v(r)=\{\begin{array}{cc}ϵ(1r/R)^2;\hfill & r<R\hfill \\ 0;\hfill & rR.\hfill \end{array}$$
(2)
It has recently been realized that binary mixtures of particles with penetrable cores, which interact via generalizations of the gaussian and parabolic potentials (1) and (2), involving different energy scales $`ϵ`$ and radii $`R`$ for the various species, may phase-separate over appropriate ranges of these parameters. Spinodal instability was first shown to occur for the gaussian core model (1) within the random phase approximation (RPA) ref17  and binodals as well as interfacial properties were then calculated within the same approximation ref18 . Similarly Gibbs ensemble Monte Carlo simulations have very recently shown that binary systems of soft particles interacting via the parabolic potential (2) phase separate beyond a critical degree of enhanced repulsion between particles of different species, in agreement with Flory-like mean field considerations ref22 .

In this paper we systematically extend our earlier results for the gaussian core model ref17  and investigate the range of validity of the RPA by detailed calculations of the pair structure, thermodynamics, and the resulting phase coexistence curve within the much more accurate hypernetted chain (HNC) approximation. The break-down of the RPA is quantified in the physically relevant regime where $`\rho R^31`$ and $`ϵk_BT`$, which would correspond to the cross-over from dilute to semi-dilute regimes of the underlying binary polymer solution. The RPA continues to provide reliable first estimates at higher densities.

## II RPA and HNC

The model under consideration is the binary gaussian core model (GCM) already introduced in references ref17  and ref18 . It consists of $`N_1`$ particles of “radius” $`R_1`$ and $`N_2`$ particles of “radius” $`R_2`$ in a volume $`V`$. The total number density is $`\rho =(N_1+N_2)/V`$, while the concentrations of the two species are $`x=N_2/N`$ and $`1x=N_1/N`$, respectively. The pair potentials are

$$v_{\mu \nu }(r)=ϵ_{\mu \nu }e^{\left(r/R_{\mu \nu }\right)^2},$$
(3)
which introduce three length and three energy parameters: $`R_{11}`$, $`R_{22}`$, $`R_{12},ϵ_{11},ϵ_{12}`$ and $`ϵ_{22}`$. If $`R_{11}`$ is chosen as unit of length, the system is entirely specified by the 5 dimensionless parameters $`R_{12}/R_{11}`$, $`R_{22}/R_{11}`$, $`ϵ_{11}^{}=\beta ϵ_{11}`$, $`ϵ_{12}^{}=\beta ϵ_{12}`$, $`ϵ_{22}^{}=\beta ϵ_{22}`$, where $`\beta =1/(k_BT)`$. For a fixed set of dimensionless parameters the reduced Helmholtz free energy per particle, $`f=F/(Nk_BT)`$, is a function only of the intensive variables $`\rho `$ and $`x`$; this may be split into the ideal gas , ideal mixing and excess (non-ideal) contributions:

$`f(\rho ,x)`$ $`=f_{\text{id}}(\rho )+f_{\text{mix}}(x)+f_{\text{ex}}(\rho ,x)`$ (4a)
$`=\mathrm{ln}(\rho \mathrm{\Lambda }^3)1+x\mathrm{ln}x+(1x)\mathrm{ln}(1x)+f_{\text{ex}}(x,\rho ),`$ (4b)
where $`\mathrm{\Lambda }`$ is an irrelevant de Broglie thermal wavelength. The equation of state $`\beta P/\rho `$, and the chemical potentials $`\mu _\nu ,(\nu =1,2)`$ follow from the standard thermodynamic relations:

$`{\displaystyle \frac{\beta P}{\rho }}`$ $`=\rho \left({\displaystyle \frac{f}{\rho }}\right)_x`$ (5a)
$`\beta \mu _1`$ $`=\left({\displaystyle \frac{\rho f}{\rho }}\right)_xx\left({\displaystyle \frac{f}{x}}\right)_\rho `$ (5b)
$`\beta \mu _2`$ $`=\left({\displaystyle \frac{\rho f}{\rho }}\right)_x+(1x)\left({\displaystyle \frac{f}{x}}\right)_\rho `$ (5c)
and may likewise be split into ideal and excess contributions. Within RPA and HNC, these excess contributions may be easily expressed in terms of the usual total and direct correlation functions $`h_{\mu \nu }(r)`$ and $`c_{\mu \nu }(r)`$, which are related by the familiar Ornstein-Zernike (OZ) relations ref23 ; these may be used to express the Fourier transforms (FT) $`\widehat{h}_{\mu \nu }(k)`$ of the total correlation functions in terms of the FT $`\widehat{c}_{\mu \nu }(k)`$ of the direct correlation functions:

$`\widehat{h}_{11}(k)`$ $`={\displaystyle \frac{1}{\mathrm{\Delta }(k)}}\left[\widehat{c}_{11}(k)(1\rho _2\widehat{c}_{22}(k))+\rho _2\widehat{c}_{12}^2(k)\right],`$ (6a)
$`\widehat{h}_{12}(k)`$ $`={\displaystyle \frac{1}{\mathrm{\Delta }(k)}}\widehat{c}_{12}(k),`$ (6b)
$`\widehat{h}_{22}(k)`$ $`={\displaystyle \frac{1}{\mathrm{\Delta }(k)}}\left[\widehat{c}_{22}(k)(1\rho _1\widehat{c}_{11}(k))+\rho _1\widehat{c}_{12}^2(k)\right],`$ (6c)
where $`\rho _1=\rho (1x),\rho _2=\rho x,`$ and

$$\mathrm{\Delta }(k)=[1\rho _1\widehat{c}_{11}(k)][1\rho _2\widehat{c}_{22}(k)]\rho _1\rho _2\widehat{c}_{12}^2(k).$$
(7)
These OZ relations must be supplemented by closure relations.

The RPA amounts to simply identifying the $`c_{\mu \nu }(r)`$ with their asymptotic (large $`r`$) behaviour:

$`c_{\mu \nu }(r)`$ $`=\beta v_{\mu \nu }(r)`$
$`=ϵ_{\mu \nu }^{}e^{\left(r/R_{\mu \nu }\right)^2}`$ (8)
Substitution of the FT’s of (II) into the relations (6) and (7) yields explicit expressions for the $`\widehat{h}_{\mu \nu }(k)`$, which may be transformed back to obtain the $`h_{\mu \nu }(r)`$. Detailed analytical expressions are given in Appendix A in the special case of mixtures with equal core radii $`R_{11}=R_{22}=R_{12}=R`$.

The HNC closure relations are ref23 :

$`g_{\mu \nu }(r)`$ $`=1+h_{\mu \nu }(r)`$ (9a)
$`=\mathrm{exp}\{\beta v_{\mu \nu }(r)+\gamma _{\mu \nu }(r)\}`$
$`\gamma _{\mu \nu }(r)`$ $`=h_{\mu \nu }(r)c_{\mu \nu }(r)`$ (9b)
An iterative procedure must be used to solve the coupled equations (6) and (9) numerically. Convergence is easily achieved for the binary GCM model over most of the potential parameter space. In fact the HNC approximation becomes exact in the high density limit ($`\rho R^3\mathrm{}`$), and is extremely accurate for densities $`\rho R^31`$, as shown earlier in the one-component case ref16 ; ref17 . Pair structure data will be examined in greater detail in section 5.

Knowledge of the pair correlation functions allows direct access to thermodynamics, via the compressibility or virial routes ref23 . Being approximate, the RPA and HNC closures are not thermodynamically consistent, i. e. the two routes lead to different answers. In the case of the GCM the two theories become strictly thermodynamically consistent only in the high density limit, ref16 ; ref17 .

Within the RPA, the compressibility route is equivalent to a simple mean-field ansatz for the free energy ref16 ; ref17 . This leads immediately to the following analytic expression for the excess part of the reduced free energy ref17 :

$`f_{\text{ex}}^\text{C}(x,\rho )`$ $`={\displaystyle \frac{1}{2}}\rho V_0(x)`$ (10a)
$`={\displaystyle \frac{1}{2}}\rho {\displaystyle \underset{\nu }{}}{\displaystyle \underset{\mu }{}}x_\mu x_\nu V_{\mu \nu }`$ (10b)
$`V_{\mu \nu }`$ $`=\beta \widehat{v}_{\mu \nu }(k=0)={\displaystyle \beta v_{\mu \nu }(r)𝑑𝐫},`$ (10c)
where $`x_1=(1x),x_2=x`$. In the special case of the GCM:

$$V_{\mu \nu }=\pi ^{3/2}ϵ_{\mu \nu }^{}R_{\mu \nu }^3$$
(11)
The pressure and chemical potentials $`\mu _\nu `$ follow directly from eqns. (5) (the superscript C refers to the compressibility route):

$`{\displaystyle \frac{\beta P^\text{C}}{\rho }}`$ $`=1+{\displaystyle \frac{\rho }{2}}{\displaystyle \underset{\mu }{}}{\displaystyle \underset{\nu }{}}x_\mu x_\nu V_{\mu \nu }`$ (12a)
$`\beta \mu _1^\text{C}`$ $`=\mathrm{ln}(\rho \mathrm{\Lambda }^3(1x))+V_{11}\rho (1x)+V_{12}\rho x`$ (12b)
$`\beta \mu _2^\text{C}`$ $`=\mathrm{ln}(\rho \mathrm{\Lambda }^3x)+V_{12}\rho (1x)+V_{22}\rho x`$ (12c)
The virial route is considerably more arduous. The pressure follows from the virial theorem:

$`{\displaystyle \frac{\beta P^\text{V}}{\rho }}`$ $`=1+{\displaystyle \frac{\rho }{2}}{\displaystyle \underset{\mu }{}}{\displaystyle \underset{\nu }{}}x_\mu x_\nu V_{\mu \nu }`$ (13)
$`{\displaystyle \frac{2\pi }{3}}\rho {\displaystyle \underset{\mu }{}}{\displaystyle \underset{\nu }{}}x_\mu x_\nu {\displaystyle _0^{\mathrm{}}}r^3{\displaystyle \frac{d\beta v_{\mu \nu }(r)}{dr}}h_{\mu \nu }(r)𝑑r`$
In the case of the symmetric mixture ($`R_{11}=R_{22}=R_{12}=R`$; $`V_{11}=V_{22}`$) the RPA pair correlation functions obtained in Appendix A may be used to yield the following expression (the superscript V referring to the virial route):

$`{\displaystyle \frac{\beta P^\text{V}}{\rho }}`$ $`={\displaystyle \frac{\beta P^\text{C}}{\rho }}{\displaystyle \frac{z_2x(1x)\mathrm{\Delta }}{\pi ^{3/2}R^3}}\mathrm{}\left({\displaystyle \frac{\rho }{z_1}}\right)`$
$`{\displaystyle \frac{z_1x(1x)\mathrm{\Delta }}{\pi ^{3/2}R^3}}\mathrm{}\left({\displaystyle \frac{\rho }{z_2}}\right),`$ (14)
where the functions $`\mathrm{}(s)`$ and the roots $`z_1`$ and $`z_2`$ are defined in appendix A (eqns. (29a) and (29b)). The resulting reduced free energy per particle is obtained by integrating the equation of state (II) with respect to density:

$$f^\text{V}(x,\rho )=f^\text{C}(x,\rho )+\frac{1}{2\pi ^{3/2}\rho R^3}\left[\mathrm{Li}_{5/2}(\rho /z_1)+\mathrm{Li}_{5/2}(\rho /z_2)+\rho V_{11}\right],$$
(15)
where the polylogarithm $`\mathrm{Li}_{5/2}`$ is defined in equation (32). Unfortunately this expression becomes singular for densities larger than the spinodal densities obtained by the compressibility route. This mathematical artifact discussed in Appendix A does not allow us to use expressions (II) and (15) to construct a phase diagram.

The lack of thermodynamic consistency of the RPA is illustrated in Figure 1, which compares the compressibility and virial equations of state as a function of density, for a symmetric mixture; the virial pressures are always lower than their compressibility counterparts, and closer to the nearly exact HNC results, to which we now turn.

The HNC approximation is very nearly thermodynamically consistent in practice, at least for temperatures and densities relevant for polymer solutions ref17 . Hence we have calculated thermodynamic properties within the more convenient virial route with the pressure given by the standard relation (13), while the chemical potentials may also be directly expressed in terms of the direct and total correlation functions, according to ref24 :

$`\beta \mu _\nu `$ $`=\mathrm{ln}(\rho _\nu \mathrm{\Lambda }^3)+{\displaystyle \underset{\mu }{}}\left\{{\displaystyle \frac{\rho _\mu }{2}}{\displaystyle 𝑑𝐫h_{\nu \mu }(r)[h_{\nu \mu }(r)c_{\nu \mu }(r)]}\rho _\mu \widehat{c}_{\nu \mu }(k=0)\right\}.`$ (16)
Note that these expressions hold only within the HNC approximation, and are consistent with the virial route (13).

## III Scaling properties

As mentioned earlier, suitably reduced equilibrium properties of the binary GCM depend on the five dimensionless combinations $`R_{12}/R_{11}`$, $`R_{22}/R_{11}`$, $`ϵ_{11}^{}`$, $`ϵ_{12}^{}`$, and $`ϵ_{22}^{}`$. The compressibility version of the RPA thermodynamics (or equivalently, the mean field approximation) allows a considerable reduction of this parameter space. From the expressions (10)-(12c) it is clear that the radii $`R_{\mu \nu }`$ and energies $`ϵ_{\mu \nu }`$ enter only in the combinations $`V_{\mu \nu }`$ defined in eqn. (11). In terms of the reduced density $`\rho ^{}=\rho V_{11}`$ and pressure $`P^{}=\beta PV_{11}`$, the thermodynamic behaviour of the mixture is uniquely determined by the dimensionless ratios $`V_{22}/V_{11}`$ and $`V_{12}/\sqrt{V_{11}V_{22}}`$. This remarkable reduction of parameter space holds only for the thermodynamics, but not for the correlation functions, which depend explicitly on the five parameters of the GCM potentials. The scale invariance of the thermodynamics is broken by the RPA virial route, as is immediately evident form the explicit expressions (II) and (15) (valid for a symmetric mixture); the same is true within the HNC approximation.

However deviations from RPA-C (mean-field) scale invariance are expected to be small in the high density, high temperature regime, where RPA is increasingly accurate, so that the considerable reduction in the number of relevant potential parameters (from 5 to 2) is expected to carry over, at least approximately, to the phase diagrams calculated with RPA-V or HNC thermodynamics, which will be presented in the next section.

## IV Phase diagrams

The phase behaviour of the binary GCM may be deduced from the knowledge of the reduced free energy per particle $`f`$ as a function of the variables $`v=1/\rho `$ and $`x`$ and from its thermodynamic derivatives (the pressure and the chemical potentials). For suitable values of the potential parameters $`ϵ_{\mu \nu }`$ and $`R_{\mu \nu }`$, the binary GCM becomes unstable against demixing at sufficiently high densities. In terms of the free energy per particle, $`f=f(v,x)`$, the standard thermodynamic stability conditions of a binary mixture read ref1 :

$`\left({\displaystyle \frac{^2f}{v^2}}\right)_x>0;\left({\displaystyle \frac{^2f}{x^2}}\right)_v`$ $`>0`$ (17a)
$`\left({\displaystyle \frac{^2f}{v^2}}\right)_x\left({\displaystyle \frac{^2f}{x^2}}\right)_v\left({\displaystyle \frac{^2f}{vx}}\right)^2`$ $`>0.`$ (17b)
The first condition ensures mechanical stability, while the second inequality guarantees stability against demixing at constant volume; the third inequality ensures stability of the mixture at constant pressure. Note that the stability at constant pressure is the more restrictive condition.

The subsequent discussion is restricted to demixing at constant pressure. The vanishing of the l. h. s. of equation (17b), corresponding to the case where the third inequality turns into an equality, signals the occurrence of spinodal instability. Within the mean-field (RPA-C) approximation (10), the equation for the spinodal is easily calculated to be ref17

$$1+\rho V_1(x)\rho ^2x(1x)\mathrm{\Delta }V=0,$$
(18)
where:

$`V_1(x)`$ $`=(1x)V_{11}+xV_{22}`$ (19a)
$`\mathrm{\Delta }V`$ $`=V_{12}^2V_{11}V_{22}.`$ (19b)
It is easily inferred from (18) that a spinodal instability occurs whenever $`V_{12}/\sqrt{V_{11}V_{22}}>1`$. This region is visualized in Figure 2 in a plot of the minimum energy ratio $`ϵ_{12}/\sqrt{ϵ_{11}ϵ_{22}}`$ required for spinodal instability, as a function of the size ratio $`R_{22}/R_{11}`$. The radius $`R_{12}`$ is taken to be given by the combination rule:

$$R_{12}^2=\frac{1}{2}(R_{11}^2+R_{22}^2),$$
(20)
as suggested by renormalization group (RG) calculations ref10  and by direct simulation ref11  of mixtures of self-avoiding polymer coils.

Figure 2 also shows the values of the ratio $`ϵ_{12}/\sqrt{ϵ_{11}ϵ_{22}}`$, as calculated for two self-avoiding walk polymers by RG techniques ref10 , plotted as a function of the ratio of their radii of gyration. The figure suggests that phase separation of real polymers would only be observed for ratios of gyration radii larger than 10, i. e. for mixtures of very long and very short polymers. However, the RG results are probably not trustworthy for such asymmetric mixtures, and the RG potentials are strictly valid only in the infinite dilution limit.

The spinodal line in the $`Px`$-plane is easily calculated from (18) to be:

$$\rho _S(x)=\frac{V_1(x)+\sqrt{[V_1(x)]^2+4x(1x)\mathrm{\Delta }V}}{2x(1x)\mathrm{\Delta }V}$$
(21)
The critical concentration is determined by the condition:

$$\frac{dP_S(x)}{dx}=0,$$
(22)
where $`P_S(x)`$ is the pressure calculated from (12a), at the spinodal density given by (21). The corresponding critical density $`\rho _C`$, calculated by substituting the critical concentration determined by (22) is plotted in Figure 3 versus the ratio $`V_{12}/\sqrt{V_{11}V_{22}}`$, for several ratios $`V_{22}/V_{11}`$.

The binodal or phase coexistence curves in the $`\rho x`$ plane are determined by the usual conditions of equality of the chemical potentials of both species and of the pressures in the two phases, using expressions (12a)-(12c). The resulting equations must in general be solved numerically to yield the binodal curves. The complete phase diagram can however be calculated analytically within the RPA-C approximation, in the symmetric case, where $`V_{11}=V_{22}`$. Due to the symmetry of the problem, all thermodynamic quantities must be invariant with respect to the transformation $`12;x(1x)`$ and $`\rho \rho `$. This simplification allows the fully analytic treatment detailed in Appendix B. The resulting spinodal and binodal are shown in Figure 4 in the case where $`V_{12}/V_{11}=2`$. The same Figure also shows the HNC binodals, for several values of $`ϵ_{11}^{}`$ (remember that within RPA-C all these binodals coincide). As expected the HNC results break the scale invariance of RPA-C and shift the binodal curves upward, to higher densities. The HNC and RPA-C critical densities differ by almost a factor of 2 for $`ϵ_{11}^{}=2`$, pointing to the limitations of the mean-field (RPA-C) description.

An example of a phase diagram in an asymmetric mixture, with potential parameters already used in earlier RPA calculations ref17 ; ref18 , is shown in Figure 5 in the $`\rho x`$ plane and in Figure 6 in the $`Px`$ plane. The spinodal and binodal curves are now asymmetric, and the difference between RPA-C and HNC results is again significant, with the HNC binodal being pushed to higher densities.

Phase separation in binary mixtures of particles interacting via the closely related parabolic core potentials (2) has been observed in the Gibbs ensemble Monte Carlo simulations of Wijmans et al. ref22 . They considered the symmetric case where all $`R_{\mu \nu }`$ are equal, $`ϵ_{11}^{}=ϵ_{22}^{}`$ and $`ϵ_{12}^{}=ϵ_{11}^{}+\mathrm{\Delta }ϵ^{}`$. We have calculated the binodal in the $`\mathrm{\Delta }ϵ^{}x`$ plane under the same conditions as the MC simulations, i. e. for $`ϵ_{11}^{}=12.5,\rho _b^{}=\rho _bR^3=3`$, in the RPA and HNC approximations. The results are compared to the MC data in Figure 7. In view of the fact that the amplitude of the parabolic repulsion is $`12.5k_BT`$, which makes it more hard-core like, the agreement may be considered to be rather satisfactory. As in the GCM case, the HNC coexistence-curve lies well above its RPA counterpart; it is closer to the simulation data.

## V Pair structure and clustering

The pair correlation functions $`h_{\mu \nu }(r)`$ of the GCM mixture can be calculated by combining the OZ relations (6) with either the RPA closure (II) or the HNC closure (9). In the symmetric case, characterized by a single gaussian range parameter $`R_{11}=R_{22}=R_{12}=R`$, analytic expressions, in the form of infinite series, can be obtained for the RPA closure, as shown in Appendix A. In the asymmetric case, solutions for the $`h_{\mu \nu }(r)`$ are readily obtained by numerical Fourier transformation of equations (6) and (7), using the straightforward Fourier transforms of the $`c_{\mu \nu }(r)`$ given by (II).

The HNC closure requires an iterative solution of the coupled OZ and closure equations. This was achieved using the standard Picard method, and well-converged solutions were generally obtained in a few iterations. Convergence was found to be slower in the vicinity of phase coexistence, and to break down rapidly inside the phase coexistence region corresponding to metastable mixtures. Concentration fluctuations build up in that region, as signalled by a $`k=0`$ peak of growing amplitude in the $`\widehat{h}_{\mu \nu }(k)`$, or equivalently in the corresponding partial structure factors

$$S_{\mu \nu }(k)=x_\mu \delta _{\mu \nu }+x_\mu x_\nu \rho \widehat{h}_{\mu \nu }(k)$$
(23)
The concentration-concentration structure factor:

$$S_{CC}(k)=x_2^2S_{11}(k)2x_1x_2S_{12}(k)+x_1^2S_{22}(k)$$
(24)
satisfies the long wavelength limit ref23 :

$$\underset{k0}{lim}S_{CC}(k)=\frac{Nk_BT}{(^2G/x^2)_{N,P,T}},$$
(25)
where $`G`$ is the Gibbs free energy. The cross-over from metastable to unstable mixture corresponds to the vanishing of $`(^2G/x^2)_{N,P,T}`$, so that $`S_{CC}(k=0)`$ is expected to diverge along a spinodal line.

This is indeed the case with the RPA closure, but the HNC closure ceases to converge before a spinodal line is reached, a well-known shortcoming of the HNC closure ref25 . Such a break-down is typical of the inadequacy of fluid integral equations to describe critical behaviour ref25a . The accuracy of HNC in the stable one-phase region for the pair correlation function $`h(r)`$, already documented in the one-component GCM ref16 ; ref17 , is tested under more severe conditions (namely $`ϵ^{}=12.5`$) in Figure 8 against recent Monte Carlo data ref22  for the parabolic core potential (2). The agreement is seen to be excellent. This suggests that the HNC closure would be very useful to generate pair structure in coarse-grained DPD fluids at a moderate computational cost.

Examples of pair correlation functions and partial structure factors for several states of GCM mixtures are shown in Figures 9-11. The most striking feature is the amount of structure observed both in $`r`$ and in $`k`$-space, even within the RPA, compared to the previously studied one-component case ref17 . The peak in the structure factors at $`k=0`$ as the spinodal is approached is expected, as explained earlier. More surprising perhaps is the appearance of peaks in the pair correlation functions as $`r0`$, or at finite $`r`$. Peaks at $`r0`$ are of course precluded in the presence of hard cores, but are quite significant in the present penetrable core model and point to a novel physical mechanism to trigger phase separation, namely the clustering of the minority species, as observed e. g. in Figure 11 for the symmetric case. The clustering can be explained by a simple “energetic” argument: the particles of the majority species tend to maximise their mutual distances. If the repulsion between unlike particles is stronger than between like ones, particles of the minority species prefer to cluster in the voids left by the majority species, rather than overlap with particles of the latter. This mechanism drives phase separation.

## VI Conclusion

The two-component extension of the gaussian core model leads to non-trivial phase behaviour, which we have investigated within RPA and HNC theories. The former is partially analytic, while the latter requires only a modest numerical effort to calculate partial pair correlation functions and the resulting thermodynamic properties. While the pair structure shows some unexpected features, the most interesting prediction is the occurrence of phase separation induced by purely repulsive pair interactions. As expected for this model of penetrable particles, phase separation is driven by an enhanced repulsion between unlike particles compared to the like-particle interaction.

A special feature of the GCM mixtures is that at high densities of both species, the mixture behaves like a “mean field fluid”, i. e. the RPA becomes asymptotically exact ref16 ; ref17 . This means that at finite densities, the RPA makes semi-quantitavely valid predictions for the phase behaviour, and allows a rapid exploration of potential parameter space to search for likely conditions for phase separation.

The nature of the link between the behaviour of the binary GCM and that of a binary polymer solution must still be worked out in detail, mainly because the effective interactions between the centres of mass of the polymer coils are state-dependent ref12 . For a one-component system at lower densities, the structure is not very sensitive to the state-dependence Bolh01a ; Bolh01b . We therefore expect the binary GCM to make qualitatively correct predictions for the structure of binary polymer solutions in the dilute and beginning of the semi-dilute regimes. The link to the phase-behaviour is less clear. Because the GCM uses a fixed potential, it is mainly relevant to behaviour in the dilute regime. While polymers in a melt are known to phase-separate rather easily ref2a , in the dilute regime for good solvents they are not expected to phase-separate ref2b . This appears to be confirmed by the results shown in Fig 2. However, it may be that for poorer solvent, or for stronger polymer incompatibilities, a regime where the polymers phase-separate at low densities could open up. If this is the case, then the binary GCM would be a useful coarse-grained model with which one could rapidly explore the qualitative phase behaviour of this regime.

Another open question relates to the critical behaviour of the binary GCM in the immediate vicinity of the critical consolute point. The RPA clearly predicts mean field exponents.As regards the exact exponents, it is not obvious whether the correct critical behaviour of the model necessarily belongs to the Ising universality class, since there is no clear correspondence between GCM mixtures and standard lattice models, due to the penetrability of the particles. The question of the correct universality class, a subject dear to the heart of Michael Fisher, must be considered as open.

## VII Acknowledgements

R. Finken is indebted to the Oppenheimer Fund for their support. A. A. Louis would like to thank the Isaac Newton Trust, Cambridge, for funding.

## Appendix A Solution of the virial route in RPA

We consider the GCM where the potentials have equal range $`R_{11}=R_{12}=R_{22}=R`$, and the heights of the repulsion potentials of like particles are the same for both species $`ϵ_{11}=ϵ_{22}`$. Without loss of generality we can set $`V_{11}=1`$. Then we have

$$\widehat{c}_{\mu \nu }(k)=V_{\mu \nu }e^{k^2R^2/4},$$
(26)
with

$$V_{\mu \nu }=ϵ_{\mu \nu }^{}R^3.$$
(27)
If we introduce the abbreviation $`\widehat{z}(k)\rho V_{11}e^{k^2R^2/4}`$ and set $`\mathrm{\Delta }V=V_{12}^2V_{11}^2=V_{12}^21`$, the OZ equations (6) reduce to (using equation (7)):

$`\rho \widehat{h}_{11}(k)`$ $`={\displaystyle \frac{x_2\mathrm{\Delta }V\widehat{z}(k)^2\widehat{z}(k)}{\mathrm{\Delta }(k)}}`$ (28a)
$`\rho \widehat{h}_{12}(k)`$ $`={\displaystyle \frac{V_{12}\widehat{z}(k)}{\mathrm{\Delta }(k)}}`$ (28b)
$`\rho \widehat{h}_{22}(k)`$ $`={\displaystyle \frac{x_1\mathrm{\Delta }V\widehat{z}(k)^2\widehat{z}(k)}{\mathrm{\Delta }(k)}}`$ (28c)
$`\mathrm{\Delta }(k)`$ $`=1+\widehat{z}(k)x(1x)\mathrm{\Delta }V\widehat{z}(k)^2`$ (28d)
The denominator $`\mathrm{\Delta }(k)`$ in equation (6) has the zeros

$`z_1`$ $`={\displaystyle \frac{1+\sqrt{1+4x(1x)\mathrm{\Delta }V}}{2x(1x)\mathrm{\Delta }V}}>0`$ (29a)
$`z_2`$ $`={\displaystyle \frac{1\sqrt{1+4x(1x)\mathrm{\Delta }V}}{2x(1x)\mathrm{\Delta }V}}<0`$ (29b)
The positive root $`z_1`$ is the spinodal density at any given composition.

The virial equation (13) gives

$`{\displaystyle \frac{\beta P^\text{V}}{\rho }}{\displaystyle \frac{\beta P^\text{C}}{\rho }}`$ $`={\displaystyle \frac{2\pi }{3}}\rho {\displaystyle \underset{\mu }{}}{\displaystyle \underset{\nu }{}}x_\mu x_\nu {\displaystyle _0^{\mathrm{}}}r^3{\displaystyle \frac{d\beta v_{\mu \nu }(r)}{dr}}h_{\mu \nu }(r)𝑑r`$ (30a)
$`={\displaystyle \frac{1}{3}}\rho {\displaystyle \underset{\mu }{}}{\displaystyle \underset{\nu }{}}x_\mu x_\nu {\displaystyle 𝑑𝐫ϵ_{\mu \nu }(r^2/R^2)e^{r^2/R^2}h_{\mu \nu }(r)}.`$ (30b)
Evaluating this integral in Fourier space and making use of the Parseval theorem leads to

$`{\displaystyle \frac{\beta P^\text{V}}{\rho }}{\displaystyle \frac{\beta P^\text{C}}{\rho }}`$ $`={\displaystyle \frac{1}{3}}\rho {\displaystyle \underset{\mu }{}}{\displaystyle \underset{\nu }{}}x_\mu x_\nu {\displaystyle 𝑑𝐤V_{\mu \nu }\left[\frac{3}{2}\frac{k^2R^2}{4}\right]e^{r^2/R^2}\widehat{h}_{\mu \nu }(k)}`$
$`={\displaystyle \frac{1}{12}}{\displaystyle 𝑑𝐤\frac{\widehat{z}(k)^2}{\mathrm{\Delta }(k)}[x_1x_2\widehat{z}(k)\mathrm{\Delta }V2x_1x_2\mathrm{\Delta }V1][6k^2R^2]}`$ (31)
This integral cannot be evaluated analytically in general. However, if the density is below the spinodal density for a given composition, one can expand the denominator in the integrand in a power series with respect to $`\widehat{z}(k)`$. The integral can then be evaluated term by term, interchanging the integration and the summation. In order to make the symmetry between the two species explicit, we write $`x=\frac{1+\xi }{2}`$. The shifted concentration $`\xi `$ now varies between $`1`$ and 1, with $`\xi =0`$ corresponding to $`x=0.5`$.

With the definition of the polylogarithmic function

$$\mathrm{Li}_\alpha (x)=\underset{n=1}{\overset{\mathrm{}}{}}\frac{x^n}{n^\alpha },$$
(32)
and the auxiliary function

$$\mathrm{}(\alpha )=\frac{1}{2\alpha }[\mathrm{Li}_{3/2}(\alpha )\mathrm{Li}_{5/2}(\alpha )]$$
(33)
evaluation of (A) leads to

$$\begin{array}{cc}\hfill \frac{\beta P}{\rho }& =1+\frac{1}{4}\rho [1+V_{12}\xi ^2\mathrm{\Delta }V]\hfill \\ & \frac{z_2(1\xi ^2)\mathrm{\Delta }V}{4\pi ^{3/2}R^3}\mathrm{}\left(\frac{\rho }{z_1}\right)\hfill \\ & \frac{z_1(1\xi ^2)\mathrm{\Delta }V}{4\pi ^{3/2}R^3}\mathrm{}\left(\frac{\rho }{z_2}\right).\hfill \end{array}$$
(34)
This can be integrated to give the free energy per particle

$$\begin{array}{cc}\hfill \beta f(x,\rho )& =f^{\text{id}}+f^{\text{mix}}(x)\hfill \\ & +\frac{1}{4}\rho [1+V_{12}\xi ^2\mathrm{\Delta }V]\hfill \\ & +\frac{1}{2\pi ^{3/2}\rho R^3}[\mathrm{Li}_{5/2}(\rho /z_1)+\mathrm{Li}_{5/2}(\rho /z_2)+\rho ]\hfill \end{array}$$
(35)
From that we get the chemical potentials as

$`\beta \mu _1`$ $`=\mathrm{log}(\mathrm{\Lambda }^3\rho )+\mathrm{log}\left({\displaystyle \frac{1\xi }{2}}\right)+{\displaystyle \frac{1\xi }{2}}\rho +V_{12}{\displaystyle \frac{1+\xi }{2}}\rho +{\displaystyle \frac{1}{2\pi ^{3/2}R^3}}`$
$`+{\displaystyle \frac{1}{2\pi ^{3/2}R^3(1\xi )\rho }}\mathrm{Li}_{3/2}\left({\displaystyle \frac{\rho }{z_1}}\right)\left[1+{\displaystyle \frac{\xi }{\sqrt{1+\mathrm{\Delta }V(1\xi ^2)}}}\right]`$
$`+{\displaystyle \frac{1}{2\pi ^{3/2}R^3(1\xi )\rho }}\mathrm{Li}_{3/2}\left({\displaystyle \frac{\rho }{z_2}}\right)\left[1{\displaystyle \frac{\xi }{\sqrt{1+\mathrm{\Delta }V(1\xi ^2)}}}\right],`$ (36)
$`\beta \mu _2`$ $`=\mathrm{log}(\mathrm{\Lambda }^3\rho )+\mathrm{log}\left({\displaystyle \frac{1+\xi }{2}}\right)+{\displaystyle \frac{1+\xi }{2}}\rho +V_{12}{\displaystyle \frac{1\xi }{2}}\rho +{\displaystyle \frac{1}{2\pi ^{3/2}R^3}}`$
$`+{\displaystyle \frac{1}{2\pi ^{3/2}R^3(1+\xi )\rho }}\mathrm{Li}_{3/2}\left({\displaystyle \frac{\rho }{z_1}}\right)\left[1{\displaystyle \frac{\xi }{\sqrt{1+\mathrm{\Delta }V(1\xi ^2)}}}\right]`$
$`+{\displaystyle \frac{1}{2\pi ^{3/2}R^3(1+\xi )\rho }}\mathrm{Li}_{3/2}\left({\displaystyle \frac{\rho }{z_2}}\right)\left[1+{\displaystyle \frac{\xi }{\sqrt{1+\mathrm{\Delta }V(1\xi ^2)}}}\right].`$ (37)
It is easily checked that $`\mu _1(\xi ,\rho )=\mu _2(\xi ,\rho )`$, due to the symmetry of the model. Note that these expressions become singular when $`\rho z_1\rho _s`$. This is a mathematical artefact originating in the interchange of summation and integration in (A), which renders the expressions for $`P`$ and $`\mu `$ unusable for binodal calculations.

## Appendix B Analytic solution of the RPA-C binodal in the symmetric case

The phase diagram simplifies considerably in the symmetric case where $`V_{11}=V_{22}`$. We express $`V_{12}`$ as $`V_{12}=V_{11}(1+\gamma )`$ and note that demixing occurs when the dimensionless parameter $`\gamma `$ is positive. As above we introduce the reduced quantities $`\rho ^{}=\rho V_{11}`$ and $`P^{}=\beta PV_{11}`$. Because of the symmetry of the problem, all quantities must be invariant with respect to the transformation $`\rho \rho ,x1x,12`$, which corresponds to a relabeling of species 1 and 2. As in appendix A, we make use of the shifted concentration $`\xi `$ to make this symmetry explicit. The pressure and the chemical potentials can now be expressed in terms of the variables $`\rho ^{}`$ and $`\xi `$:

$`P^{}`$ $`=\rho ^{}+{\displaystyle \frac{1}{4}}\rho _{}^{}{}_{}{}^{2}(2+\gamma ){\displaystyle \frac{1}{4}}\rho _{}^{}{}_{}{}^{2}\gamma \xi ^2`$ (38a)
$`\beta \mu _1`$ $`=\mathrm{ln}(\rho ^{})+\mathrm{ln}\left({\displaystyle \frac{1\xi }{2}}\right)+{\displaystyle \frac{1}{2}}\rho ^{}(2+\gamma )+{\displaystyle \frac{1}{2}}\rho ^{}\xi \gamma `$ (38b)
$`\beta \mu _2`$ $`=\mathrm{ln}(\rho ^{})+\mathrm{ln}\left({\displaystyle \frac{1+\xi }{2}}\right)+{\displaystyle \frac{1}{2}}\rho ^{}(2+\gamma ){\displaystyle \frac{1}{2}}\rho ^{}\xi \gamma `$ (38c)
In these equations we have set $`\mathrm{\Lambda }=V_{11}`$, which merely shifts the free energy per particle by a constant amount and thus does not change the phase diagram.

Two coexisting phases $`\alpha ,\beta `$ must have the same pressure and chemical potentials, respectively:

$$P^\alpha =P^\beta ,\mu _1^\alpha =\mu _1^\beta ,\mu _2^\alpha =\mu _2^\beta $$
(39)
Because of the symmetry in this model, coexisting phases must have the same density $`\rho ^\alpha =\rho ^\beta `$ and the corresponding concentrations must fulfill $`x^\alpha =1x^\beta `$. Considering this, equal pressure in both phases is trivially fulfilled, and the two conditions of equal chemical potentials become equivalent. The relation between the binodal density $`\rho _b`$ and concentration is

$$\rho _{}^{}{}_{b}{}^{}=\frac{1}{\gamma \xi }\mathrm{ln}\left(\frac{1+\xi }{1\xi }\right).$$
(40)
Expanding the binodal densities around the critical point gives

$$\rho _{}^{}{}_{b}{}^{}=\frac{2}{\gamma }\left[1+\frac{1}{3}\xi ^2+𝒪(\xi ^4)\right].$$
(41)
From the series expansions we see that for all $`\gamma >0`$ the spinodal lies inside the binodal. Within the RPA approximation we thus get an analytically exact solution for a fluid-fluid demixing transition of a liquid consisting of purely repulsive particles. As expected in this treatment, the resulting critical exponents are mean-field exponents.


format: text

text:
# Luttinger-liquid-like behavior in bulk crystals of the quasi-one-dimensional conductor NbSe3

## I Introduction

Transport properties of electron systems strongly depend on dimensionality. In one dimension the Coulomb interaction between carriers modifies the ground state and dramatically changes the transport properties. In a one-dimensional (1D) metal, instead of a Fermi liquid of electrons, a so-called Luttinger liquid (LL) is developed LLReview . Elementary excitations in a LL are collective charge and spin boson modes, rather than fermion quasiparticles as in 3D metals. The ground state of a LL is characterized by a power-law suppression of the tunneling density of states. Impurities in a LL act as tunneling barriers. Therefore, the conduction $`G`$ of impure LLs is a power-law function of temperature, $`GT^\alpha `$, and voltage, $`GV^{\beta 1}`$ ($`IV^\beta `$ in current-voltage curves) KF . These dependencies are considered as a fingerprint of 1D behavior.

The experimental study of 1D effects in long nanowires remains a challenge. The main efforts are focused on developing a technology to fabricate metal or semiconductor quantum wires with transverse sizes on the order of 10 nm or less and attaching them to an external circuit. At the present, both problems were successfully solved for a number of physical systems including single-wall carbon nanotubes carbon , Mo6O6 molecular wires Mo6O6 , Si whiskers Si , InSb nanowires in asbestos matrices InSb , Bi nanowires in various matrixes (see e.g. Bi\_Al2O3  and references therein), etc. Experimental observation of the power-law-supporting LL theories has been reported for single-wall carbon nanotubes carbon ; Sonin , Mo6O6 nanowires Mo6O6nano  and InSb nanowires InSb . In all cases the nanowire diameters do not exceed 10 nm.

There is a wide class of materials suitable for experimental search and study of 1D effects: quasi-one dimensional (quasi-1D) conductors with a charge-density wave (CDW) ground state CDWReview . Such materials consist of metallic chains connected in a crystal through the relatively weak Van der Waals interaction. At temperatures below the Peierls transition temperature, $`T_p`$, the Peierls instability leads to the development of the CDW, which is a 3D-ordered periodic modulation of the electron density accompanied by a periodic distortion of the crystalline lattice.

An example of these materials is NbSe3. It has a monoclinic unit cell with lattice parameters $`a=10.01`$ Å, $`b=3.48`$ Å, $`c=15.63`$ Å, and $`\beta =109.5^o`$ CDWReview . The conducting chains are along the $`b`$ direction. NbSe3 exhibits two Peierls transitions leading to the partial dielectrization of the electron spectrum, at $`T_{P1}=145`$ K and $`T_{P2}=59`$ K. The remaining non-condensed electrons provide metallic conduction down to the lowest temperatures.

In some sense, the CDW state results from an instability of LLs with respect to 3D coupling at low enough temperatures. It is also widely accepted that in quasi-1D conductors LL exists at temperatures above $`T_P`$ (see e.g. Voit ). As $`T_p`$ decreases with the sample cross-section area BZZN , one could expect the area of existence of LL to grow towards the low-temperature region for sufficiently thin crystals.

Transitions from quasi-1D to 1D behavior with reduction of the transverse sample sizes has recently been observed in the quasi-1D conductors NbSe3 and TaS3 ZZPM . Here we present the results of our observation of even more dramatic changes in the conduction of quasi-1D conductors found in two novel (for quasi-1D conductors) types of structures: CDW/normal-metal/CDW (C/N/C) junctions and nanoconstrictions. Surprisingly, 1D behavior has been observed in samples with transverse sizes in the $`\mu `$m-range scale.

## II Technology

Structures were prepared at the Delft Institute for Microfabrication and Submicron technology (DIMES). We used a standard Fei Inc. 200xP focus-ion-beam (FIB) system with a gallium ion source. This FIB system was used for imaging, etching and platinum deposition in situ. In all cases the FIB was used at 30 kV. When imaging, the lowest beam currents were used ($`I=1`$ \- 4 pA) to minimize possible damage of incoming gallium ions.

For the fabrication of the CDW devices, a long, narrow and flat quasi-1D crystal is placed on top of a sapphire substrate with preliminary deposited gold contacts. For nanoconstriction fabrication, the FIB is used to remove part of a crystal between two gold terminals (Fig. 1(a)).

C/N/C junctions were fabricated in two steps. In the first step the FIB is used to make a gap in the crystal (right part of Fig. 1(b)). In the second step the gap is filled in situ with platinum by FIB deposition (left part of Fig. 1(b)). In some cases additional platinum voltage terminals were attached to the crystal using FIB platinum deposition.

## III Results

### III.1 Nanoconstrictions

Damage during etching may have important consequences for the interpretation of the data. However, up to now we have not found any indication for serious damage to our samples due to FIB milling. For example, Shapiro step measurements on samples etched with the FIB still show complete mode-locking indicating their good quality. In addition, for NbSe3 the resistance ratios, $`RRR=R(300\mathrm{K})/R(4.2\mathrm{K})`$, of etched sample parts are the same as the ones of unetched parts.

Reduction of the transverse sample sizes reveals interesting finite-size effect (Fig. 2). First of all, when the resistance per unit length, $`\rho =R(300\mathrm{K})/L`$, of thin cleaved (none-etched) crystals of NbSe3 exceeds 0.1 k$`\mathrm{\Omega }/\mu `$m, the resistance ratio, RRR, starts to decrease. When $`\rho `$ approaches to 1 k$`\mathrm{\Omega }/\mu `$m, the low-temperature metallic conduction ($`dR/dT<0`$) converts into nonmetallic behavior ($`dR/dT<0`$). Further reduction of the transverse sample sizes done by the FIB leads to further changes in the entire $`R(T)`$ dependence: disappearance of both Peierls transitions and appearance of a nonmetallic behavior which can be approximated by the power law, $`RT^\alpha `$, with $`\alpha =1`$ \- 3. Finally, the $`R(T)`$ curves of the thinnest crystals have nothing in common with the initial $`R(T)`$ of thick samples.

The inset in Fig. 2 shows a typical low-temperature I-V curve of a nanoconstriction. It is nonlinear and the nonlinear part can be approximated by the the power law, $`IV^\beta `$ with $`\beta 3`$.

### III.2 C/N/C junctions

Transport measurements of C/N/C junctions have been done in a two-terminal configuration. Before a junction is fabricated into a crystal, the $`R(T)`$ is measured and demonstrates the expected behavior for NbSe3. After fabrication of a metallic junction, the $`R(T)`$ measurement exhibits a steep rise in resistance below 160 K, as illustrated in Fig. 3. At low temperatures the resistance increases dramatically and at $`T10`$ K it has increased over 6 orders of magnitude when compared to the room temperature resistance. Below 100 K, the logarithm of conductance is a linear function of the logarithm of temperature (Fig. 3), i.e. $`G(T)T^\alpha `$. Among 10 studied C/N/C junctions, three exponents lies in the region $`2.9\alpha 3.4`$, and seven belong to the interval $`5.5\alpha 6.5`$.

Fig. 4 shows a typical temperature set of current-voltage characteristics (I-V curves) of a C/N/C junction. At temperatures above 100 K the I-V curves are linear. Upon cooling, the I-V curves become more and more nonlinear. Finally, at the lowest temperature studied, $`T=10`$ K, nonlinear growth of the conduction, $`G=I/V`$, by four orders of magnitude is observed.

Four-terminal samples have been fabricated to investigate to what extend the interface between NbSe3 and platinum is important for the observation of behavior described above. For this purpose two platinum probes are deposited on top of the crystal surface by FIB deposition 40 $`\mu `$m apart of each other. At this stage temperature variations of the sample resistance have been measured to verify its character before junction fabrication. Then a junction is fabricated on some distance $`L`$ outside this 40 $`\mu `$m region. The platinum probes are used as voltage probes in a four-terminal setup. It has been found that for $`L20`$ $`\mu `$m the junction fabrication does not affect the sample properties, whereas for $`L=5`$ and 10 $`\mu `$m a rapid increase in resistance is observed when going down in temperature. These measurements rule out the possibility that the resistance growth is solely an interface phenomena. It was also found that both etching alone and deposition of platinum on the top of the crystal alone do not dramatically change the CDW properties of crystals as measured by electric transport. Thus, the resistance increase outside the junction region only occurs when the crystal has been etched through and platinum has subsequently been deposited.

## IV Discussion

Our results show that quasi-one dimensional conductors with a CDW can demonstrate the temperature and electric-field dependent conduction expected for 1D systems. The electric properties of nanoconstrictions in NbSe3 reproduce quantitatively the results obtained for long thin crystals of NbSe3 ZZPM . The electric properties of C/N/C junctions reveal qualitatively the same type of behavior, but for much bigger transverse sizes of samples: 1D behavior is observed even in samples with a width of the order of 10 $`\mu `$m.

A quasi-1D conductor can be considered as a set of metallic chains with LL coupled by an interchain hopping $`t_{}`$. At $`T<t_{}`$ the LLs are unstable towards 3D coupling Voit , and the 3D ordered CDW may be formed.

Reduction of the transverse sizes of quasi-1D conductors decreases the stiffness of the CDW. This results in growth of dynamic CDW fluctuations. In NbSe3 a growth of fluctuations is noticeable in thin samples with submicrometer transverse sizes Thorne . At a critical transverse size, when the amplitude of fluctuation of a CDW phase difference between neighboring chains reaches a value of the order of $`\pi `$, the stiffness of the shear deformation modes vanishes. This leads to destruction of the 3D order of the CDW, and the system starts to behave as a set of LLs, like it is expected in a sliding LL phase at sufficiently small $`t_{}`$ 2DcoupledLL . An additional contribution may come from static CDW disorder. Growth of static CDW fluctuations may result from an increase of the surface contribution, as well as from introduction of defects and impurities. In this case one can also expect frustration of the formation of a 3D-ordered CDW and recovery of LL properties. It has been argued that the power laws expected for LLs, survive for the sliding LL phase formed in 2D and 3D ordered chains of coupled LLs, at least in some range of $`t_{}`$ 2DcoupledLL ; SlidingLL .

If the scenario described above is correct, one could expect that the measured temperature dependencies of the linear resistance, as well as I-V curves, can be fitted by LL-like equations. According to Ref. MasterCurve , for a single-channel LL the temperature set of I-V curves collapses into a master curve described by:

$$\frac{I}{T^{\alpha +1}}=C\mathrm{sinh}\left(\gamma \frac{eV}{kT}\right)\left|\mathrm{\Gamma }\left(1+\frac{\beta }{2}+i\frac{\gamma }{\pi }\frac{eV}{k_BT}\right)\right|^2,$$
(1)
where $`\mathrm{\Gamma }`$ is the gamma-function, and $`C`$ and $`\gamma `$ are the fitting parameters. One could expect that a similar equation, but with modified $`\alpha `$ and $`\beta `$, describes stacks of LLs. Fig. 5 shows the best fit of the data of Fig. 4 by Eq. 1 with $`C=5\times 10^{20}`$, $`\gamma =0.108`$, $`\alpha =6`$, $`\beta =4`$. Apparently, Eq. 1 describes the data well over the entire studied temperature and voltage range.

The conduction of a set of coupled LLs is provided by both tunneling along LLs through impurity barriers and tunneling between the chains with LLs. The conduction exponents are $`1/g1`$ for impurity barriers, and $`(g+1/g2)/2`$ for interchain tunneling TL\_parameters . Both conduction channels are connected in parallel. The relative contribution of each conduction channel depends also on the strength of impurities, temperature and the electric field. As the resistance along the chains grows faster upon cooling than the resistance across the chains, interchain tunneling dominates at low temperatures.

It is interesting to compare the observed exponents with the predictions of LL theories. The biggest observed exponent is $`6`$. In the case of intrachain tunneling we get 1/g=7. Taking $`1/g\sqrt{1+U/2E_F}`$ LLReview ; KF , where $`U`$ denotes the potential energy of electron-electron interaction and $`E_F`$ is the Fermi energy, one gets $`U100E_F`$. As $`E_F1`$ eV in NbSe3 CDWReview , the respective $`U10^2`$ eV is too big (In the case of interchain tunneling, the discrepancy is even more dramatic). Similar discrepancies are observed if we use a more accurate way of estimating $`g`$ TL\_parameters . In accordance with Ref. TL\_parameters , to get $`\alpha 6`$, we need $`k_Fd0.1`$, where $`d`$ is the effective wire diameter, and $`k_F`$ is the Fermi wave vector. As $`k_F=Q_{CDW}/2=0.45`$ Å, where $`Q_{CDW}`$ is the CDW wave vector, the diameter $`d0.2`$ Å. Such an effective wire diameter is too small. Also for this estimation, an even more serious disagreement arises if we assume interchain tunneling as the dominant mechanism of conduction.

The observed disagreement with predictions of LL theory for a single chain may be caused by a renormalization of exponents by interchain interactions 2DcoupledLL ; SlidingLL . We would like to note that the CDW-like component of the electron density SabGind  —which is expected to be large in our quasi-1D crystals — can increase dielectrization. Additional dielectrization (especially for C/N/C junctions) may come from Friedel oscillations around impurity sites FrOscComment ; RomEgg . Further study is necessary to clarify the origin of the observed behavior.

## V Conclusion

Our results show that quasi-1D conductors with a CDW ground state can demonstrate the temperature and electric-field dependent conduction expected for 1D systems. For FIB-etched thin crystals we reproduced results obtained earlier on thin un-etched crystals ZZPM . For C/N/C structures, the LL features are more pronounced: we observed power laws up to 6 orders of magnitude in conduction. The observed behavior corresponds to the destruction of 3D CDW ordering by fluctuations and the recovery of the Luttinger liquid features of metallic chains.

## VI Acknowlegments

We are grateful to Yu.A. Firsov, V.A. Sablikov and Yasha Gindikin for valuable discussions, and to R.E. Thorne for providing batches of high-quality crystals of NbSe3. This work is supported by the Dutch Foundation for Fundamental Research on Matter (FOM), the Dutch Organization for Scientific Research (NWO), the Russian Foundation for basic research (01-02-17771), and by MNTP “Physics of Solid State Nanostructures” 97-1052. HSJvdZ is supported by the Dutch Royal Academy of Arts and Sciences (KNAW).


text:
# The evolution of the cluster environments of radio sources at 𝑧<1.8

## 1. Introduction

The discovery and study of clusters at redshifts beyond one is of great importance for many aspects of structure formation and cosmology, as discussed throughout this volume. Unfortunately, identification of clusters at these redshifts is very challenging observationally: at X-ray wavelengths existing surveys suffer sensitivity limits and the field of view of the current generation of X-ray telescopes make them inefficient for wide-area surveys; at optical and near–infrared wavelengths distant clusters show only low contrast above the high background counts at faint magnitudes. Fortunately an alternative method of locating distant clusters does exist, using powerful radio sources as probes.

Radio sources are hosted by giant elliptical galaxies (e.g. McLure & Dunlop 2000, Best et al 1998) and are amongst the most massive galaxies known in the early Universe. At low redshifts they are typically found in galaxy groups, or poor clusters, but environmental richness appears to increase with redshift: at $`z0.5`$, the large amplitude of the galaxy cross–correlation function around radio galaxies and an Abell clustering classification indicate that about 40% of radio sources are located in clusters of Abell richness class 0 or greater (e.g. Hill & Lilly 1991). At $`z1`$ there is overwhelming evidence from observations at a wide variety of wavelengths that at least some powerful radio sources are located at the centres of clusters (e.g. Best 2000 and references therein). Best (2000) analysed UKIRT K–band data of a sample of 3CR radio galaxies at $`z1`$ and found a sharp peak in the cross–correlation function surrounding the radio galaxies, corresponding to a mean environmental richness of between Abell class 0 and 1. Colour–magnitude (C-M) relations showed red sequences around some (but not all) radio sources (see also Barr et al, this volume), but a detailed analysis of the C-M relations and the source–to–source variations was prohibited by the small field of view (70 by 70 arcsec) of the UKIRT frames.

At still higher redshifts, strong evidence is found for (proto-)clusters around radio galaxies at $`z>2`$ (e.g. Pentericci et al, Kurk et al, this volume). At these redshifts, however, studies have been limited to objects selected by line emission: the sequence of red galaxies characteristic of low redshift clusters is not seen.

In this article the first results are presented of a very deep optical–IR imaging study of the environments of radio sources at $`1.44<z<1.7`$, in order to investigate their environments, to determine to what redshift the centres of clusters are dominated by a population of red galaxies, and ultimately to use this red sequence to investigate the evolution of the cluster ellipticals.

## 2. The fields of radio galaxies at $`z1.6`$

We have defined a complete subsample of 9 radio sources with redshifts $`1.44<z<1.7`$ and galactic latitudes $`|b|>20^{}`$ from the equatorial sample of Best et al (1999; $`S_{408\mathrm{M}\mathrm{H}\mathrm{z}}>5`$ Jy; $`30<\delta <+10^{}`$). Deep 5 by 5 arcminute R, J, K images of the fields of five of these sources were obtained in Sept 2000 and Aug 2001 using the NTT, reaching typical depths of $`K20.5`$, $`J22.3`$ and $`R26`$. Although final results await the completion of the imaging program, these multi–colour data provide some promising early results, including:

$``$ An excess of K–band counts. The number counts in the infrared K–band have been determined as a function of magnitude in each field (see Table 1). These results show a clear excess of number counts fainter than $`K18.5`$. Note that the magnitude of the radio source host galaxy is typically $`K18`$.

$``$ A peak in the angular cross–correlation function around the radio sources. For each field the angular cross–correlation function has been determined following the method outlined by Best (2000); cross–correlation statistics have been calculated both for all galaxy–galaxy pairs, and for only radio galaxy–galaxy pairs. Combining the data from the different fields, there is a clear excess of clustering around the radio host galaxies, as compared to the field in general (Figure 1), the mean cross–correlation amplitude indicating cluster environments of Abell richness class 0 to 1. To investigate how the cross–correlation amplitude varies with galaxy colour, the data were split into different colour bins. The clustering excess is predominantly associated with red galaxies; those galaxies with colours consistent with being cluster ellipticals are very strongly clustered around the radio galaxy, whilst blue galaxies show no excess clustering. There is some field-to-field variation; this will be investigated when the dataset is completed.

$``$ Overdensities of red galaxies. In many of the fields there is a large overdensity of extremely red galaxies (see also Hall et al 2001), with colours consistent with passively evolving elliptical galaxies at the radio source redshift. There is especially strong clumping of red galaxies on a $`100`$kpc scale. This can be seen, for example, in colour cuts as a function of radius (Figure 2).

## 3. Conclusions

We have presented some of the first results of a program of optical–IR imaging of the fields of radio galaxies at redshifts $`z1.6`$. Three different analyses all provide strong evidence for significant galaxy overdensities around some of these radio sources. Work to spectroscopically confirm these redshift and to study the properties of the cluster galaxies is continuing. The presence of the extremely red galaxies is particularly exciting, as it indicates that a substantial population of elliptical galaxies have already formed, and must have ceased star formation a significant time ($`\mathrm{}>1`$ Gyr) before redshift 1.5. This red sequence also makes these amongst the most distant ‘normal’ clusters known, and thus offers unique prospects for advancing studies of galaxy formation and evolution.

## References

Best P.N., Longair M.S., Röttgering H.J.A., 1998, MNRAS, 295, 549

Best P.N., Röttgering H.J.A., Lehnert M.D., 1999, MNRAS, 310, 223

Best P.N., 2000, MNRAS, 317, 720

Hall P.B., et al., AJ, 121, 1840

Hill G.J. & Lilly S.J., 1991, ApJ, 367, 1

McLure R.J. & Dunlop J.S., MNRAS, 317, 249


text:
# Lattice Model for Approximate Self-Affine Soil Profiles

## 1 Introduction

The concepts of the fractal geometry have been widely used to describe and quantify irregularity in nature, and statistical self-affine properties have recently been identified in various Earth’s terrain  and profiles . Several works in soil science incorporates the fractal geometry to describe and modeling soil physical properties, soil physical processes and quantify the soil spatial variability .

A very important property closely related to the fractal geometry is the soil surface roughness, defined as the configuration of the soil microrelief. The soil surface roughness exerts great influence on water infiltration, erosion and run-off effects. Its quantification is important for understanding the soil behavior during degradation processes like rainfall erosion or abrupt changes such as those induced by tillage . In the last years, a considerable effort was done to simulate the soil structure; several models propose the simulation of the particle size distributions , the soil surface roughness , the morphology of the pore-solid structure , etc. The major properties considered in these models are the fractal dimension of the soil surface, $`D`$, the particle size distribution (PSD), the pore size distribution and the surface roughness.

In this work we present a simple lattice model to simulate the soil structure and reproduce the surface soil roughness. We basically use the ideas of the fractal growth models  to generate an approximate self-affine profile, that exhibit, as much as possible, a similitude with the fractal properties of real soils. It is approximate in the sense that the scaling properties of these profiles are valid only in a limited range of scales. Therefore, the model was validated by means of a comparison between the fractal dimension estimated for simulated profiles and for natural soil surfaces. The major improvements in our model are the power law distribution of the sizes of the deposited particles or aggregates and a random selection in the allowed morphologies of the aggregates. This two features, not improved before by any model, are responsible for the soil structure in our model. In the section II we present a detailed description of the model, showing the particle size distribution and explaining the algorithm for the random choice of the morphologies. We also present a brief summary of the theory of growth surfaces. In the section III is presented the results for several maximum aggregates sizes and different fragmentation dimensions. Finally, in the section IV we present our conclusions and perspectives.

## 2 The Model

The motivation in elaborate this model lies in aptitude to simulate some of the majors properties of agricultural soil, like self-affine profiles and porous medium. The model have to allow the possibility to change some parameters, as the aggregate size distribution, and the maximum mass of the aggregates. These two parameters try to cover the structural variability found in natural soils. The aggregate size distribution is assumed to have power law behavior , and their exponent, $`D_f`$, is one of the parameters of the model. The other parameter, $`M_{max}`$, is the maximum mass of particles utilized to generate the aggregate, which is approximately related with the square of maximum aggregate radius. In spite of others classical lattice models, where the particles generated always have the same morphology, this model simulate the variability observed in natural soil, generating aggregates by a bond percolation algorithm. All this characteristics have to be integrated in a very simple model to allow the simulation of relatively large systems ($`L=1000`$).

In this way, the major improvement to choose a lattice model is the simplicity and the velocity to run the code, what possibilities the simulation of large systems with great quantity of deposited particles (hundred millions). This massive number of deposited particles are needed in order to avoid the effect of the constant attachment of micro random variability, represented by the random morphology of the aggregates. Another advantage to work in the lattice, is the theoretical background furnished by the fractal growth theory, presented by Barabási and Stanley .

We consider the soil structure composed by a set of particles and aggregates whose radius are power law distributed according to the following relation, known as Turcotte’s empirical law ,

$$N(r>R)R^{D_f},$$
(1)
where $`N(r>R)`$ is the cumulative number of particles (or aggregates) with radius $`r`$ greater than $`R`$ and $`D_f`$ is the fragmentation dimension of the particle size distribution (PSD). There is some controversy in the literature about the allowed range of the values of the $`D_f`$. Tyler & Wheatcraft  argued that $`D_f<3`$ because, under usual hypothesis (constant density, spherical particles, etc), the mass distribution scales with $`M(r>R)R^DR^{3D_f}`$, where $`D`$ is the fractal dimension of the soil. Thus, considering a fractal distribution, only the values $`D_f<3`$ have physical meaning. Martin & Taguas  present several mathematical arguments supporting this conjecture, and show some PSD simulations. However, there are several experimental studies, summarized by Perfect & Kay  where the range of values of $`D_f`$ is $`2.6<D_f<3.5`$. Gimenez et al  also affirm that there is not any experimental relation between the fractal dimension and the fragmentation dimension. We consider that the usual hypothesis of constant density of the particles is not valid when the radius of the aggregates grows up, due the presence of pores in the structure. Thus, fragmentation dimensions greater than 3 are, in principle, allowed.

According to the USDA classification of soil texture , the basic particle size classification is:
$``$ sand $`50<r<2000\mu m`$
$``$ silt $`2<r<50\mu m`$
$``$ clay $`r<2\mu m`$.
The real soils can vary widely in the percentile of each range of sizes. The study of Nemes et al proposes a standardisation of the classification of the European soils and shows several experimental data . In the figure 1, we present a typical PSD generated by the algorithm of the model. This PSD is constructed considering that the particles are approximately circles; so, there is a direct relationship between the mass number and the particle radius, that is used to build the PSD.

The model has two parameters, the fragmentation dimension, $`D_f`$, and the maximum mass of the aggregates, $`M_{max}`$. To every particle is selected a mass number, $`M`$; when $`M>1`$, the model choose a random configuration for it; this choice is one of the possible bond percolation clusters with a given size. In the figure 2, we show some of the possible aggregate morphologies. These two parameters try to cover, in a simple way, the structural variability found in natural soils.

The deposit of each particle follows the model of the solid-on-solid deposition with surface relaxation . The initial position for the deposition is random chose and the particle follows a straight line until touch the soil surface. Then, the algorithm simulate a surface relaxation (without rotation) to the particles, which only are adhered to the bulk when they reach the site with the local minimum energy .

Considering the theory of growth surfaces, the solid-on-solid deposition belongs to the universality class of the Edwards-Wilkinson equation

$$\frac{h(x,t)}{t}=\nu ^2h+\eta (\stackrel{}{x},t),$$
(2)
where $`h(\stackrel{}{x},t)`$ is the height of the site $`\stackrel{}{x}`$ at time $`t`$, and $`\eta (\stackrel{}{x},t)`$ represents a random noise, that expresses the fluctuations in the arrive of the particles and had the following properties, $`<\eta (\stackrel{}{x},t)>=0`$ and $`<\eta (\stackrel{}{x},t)\eta (\stackrel{}{x^{}},t^{})>=2C\delta ^d(\stackrel{}{x}\stackrel{}{x^{}})\delta (tt^{})`$.The $`\delta `$ is the usual Kronecker delta and $`C`$ is a diffusion constant.

The width of the growth surfaces, $`w(L,t)`$ then obeys a scaling relation

$$w(L,t)L^\alpha f(t/L^z),$$
(3)
where $`f(u)`$ is a scaling function: $`f(u)u^\beta `$ for $`u<<1`$ and $`f(u)`$ is a constant for $`u>>1`$; $`\alpha `$, $`\beta `$ and $`z`$ are the roughness, the growth and the dynamical exponents, respectively. They are known as the scaling exponents and are related by the scaling law

$$z=\frac{\alpha }{\beta }.$$

The EW scaling exponents are known exactly: $`\alpha =1/2`$, $`z=2`$.

Another equation related to the deposition is the Kardar-Parisi-Zhang equation, that has a nonlinear term,

$$\frac{h(x,t)}{t}=\nu ^2h+\lambda /2(h)^2+\eta (\stackrel{}{x},t).$$
(4)
To this equation the calculeted exponents are $`\alpha =1/2`$, $`z=3/2`$.

So, we perform simulations to test the influence of the model parameters $`D_f`$ and $`M_{max}`$ in the scaling exponents and fractal dimension of the soil profiles, and compare the results with the experimental data disposable.

## 3 Results

We perform simulations in one dimensional lattice with length $`L=1000`$ sites. The results for the scaling exponents represent the average value in a set of 20 samples with the same parameters in each simulation, and different deposition sequence of the random aggregates. The number of particles/aggregates deposited in each sample is $`N=100000000`$. The behavior of the local width, $`w(L,T)`$, defined by

$$w^2(L,t)=\frac{1}{L}\underset{i=1}{\overset{L}{}}\left(h_i(t)\overline{h}(t)\right)^2,$$
(5)
with the model parameters are shown in the figure 3. Note the roughness saturation at $`N10000000`$. The fractal dimension of the final profile is estimated using the Hurst exponent $`H`$. The Hurst exponent is associated to the fractal dimension via the relation $`D=2H`$. There are several ways to measure $`H`$. At this work, we use the detrended fluctuation analysis (DFA), improved by the first time by Moreira et al , and the semivariogram method, that uses a height-height correlation function.

The DFA consists in measure the roughness around the mean square straight line . The roughness $`W(L,ϵ,t)`$ at the scale $`ϵ`$, is given by

$$W(L,ϵ,t)=\frac{1}{L}\underset{i=1}{\overset{L}{}}w_i(ϵ,t)$$
(6)
and the local roughness $`w_i(ϵ,t)`$ is defined by

$$w_i^2(ϵ,t)=\frac{1}{2ϵ+1}\underset{j=iϵ}{\overset{j=i+ϵ}{}}\{h_j(t)[a_i(ϵ)x_j+b_j(ϵ)]\}^2$$
(7)
where $`a_i(ϵ)`$ and $`b_i(ϵ)`$ are the linear fitting coefficients to the displacement data on the interval $`[iϵ,i+ϵ]`$ centered at the site $`i`$. Self affine profiles satisfy the scaling law

$$W(ϵ)ϵ^H$$
(8)
that is used to measure $`H`$. The semivariogram method estimates the spatial variability, through the calculus of the semivariance as a function of the distance between points. This function is can be estimated by:

$$\gamma (ϵ,t)=\frac{1}{2n(ϵ)}\underset{i=1}{\overset{n(ϵ)}{}}[h_i(t)h_{i+ϵ}(t)]^2,$$
(9)
where $`h_i(t)`$ the height in the location $`i`$ at time $`t`$, and $`n(ϵ)`$ represents the number of pairs of points which are separated by $`ϵ`$. In the case of self-affine profiles $`\gamma `$ exhibit a power law behavior,

$$\gamma (ϵ)ϵ^{2H},$$

and, in the same manner as in the DFA method, the $`H`$ exponent can be utilized to calculate the fractal dimension.

The semivariogram method is especially useful in cases where the data are not regular spaced, and for this reason, is widely used in estimating fractal dimension of soil surfaces .

The range of the $`H`$ values is shown in the figure 4A and 4B. In the figure 4C and 4D, we show the dependence of the growth exponent with the model parameters. We conclude that the increase of the maximum particle mass $`M_{max}`$ have similar effects to the decrease of the fragmentation dimension: both introduces a nonlinear correlation in the system, expressed by the increase in the value of the growth exponent $`\beta `$, changing the universality class of the deposition. Nevertheless, the fractal dimension of the surface do not alters considerably with the model parameters.

## 4 Conclusions and Perspectives

In this work, we present a new model to simulate the soil structure and its surface roughness. Two parameters are used to control the simulated profiles: the maximum mass of the particles and the fragmentation dimension. The major improvements of this model are: the random configuration allowed to the particles or aggregates and the power law distribution of its radius. These two features are not present in any model discussed in the literature until today, and permits the reproduction of the variability observed in natural soils. The results obtained shows a good agreement for the exponent $`H`$ calculated from the simulations and measured experimentally . We also observe a dependence of the growth exponent $`\beta `$ with the maximum mass allowed to the particles: as the $`M_{max}`$ grows, the value of the $`\beta `$ exponent grows from the EW value ($`\beta =1/4`$) to the KPZ value ($`\beta =1/3`$). So, the increase of the averaged particle radius corresponds to the introduction of nonlinear correlations into the system.

Next, we expect increase the system size and control better the shape of the deposited particles, in order to perform simulations closer to the real soils. We also intend simulate the rainfall effect over the simulated structure, to verify the dependence of the Hurst exponent with the rain .

This work was sponsored by the brazilian agencies CNPq, FAPEMIG and FINEP, and by the spanish agency AECI.