text:
# Sedimentation profiles of systems with reentrant melting behavior

## I Introduction

Colloidal particles in a suspension under gravitational influence show spatial inhomogeneities due to the symmetry breaking induced by the gravitational field. The problem of sedimentation of particles in the presence of gravity has been of long scientific interest. The simplest approximation is the one of non-interacting particles, valid in the limit of dilute solutions. This approach leads to an exponential sedimentation density profile, which was observed by Perrin for a calculation of Boltzmann’s constant in 1910 Perrin (1910). Taking into account particle interactions at higher concentrations will yield corrections to the exponential density profile. For very small gravitational strength, a local-density-approximation (LDA) of density functional theory (DFT) is justified Biben et al. (1993); Löwen (1994). In this case, there is a one-to-one correspondence between the sedimentation density profile and the isothermal equation of state. This fact was exploited to extract the hard sphere equation of state experimentally by investigating sterically stabilized colloids Piazza et al. (1993). Furthermore, within the LDA, a change in the height $`z`$ corresponds to a local change of the chemical potential $`\mu `$ of the bulk system. This implies that, in the limit of small gravity, the phase behavior becomes visible as a function of height $`z`$, a feature which was also been exploited to estimate the hard sphere freezing transition Piazza et al. (1993). Surprisingly, comparison with Monte-Carlo (MC) simulations show that the LDA is even reliable for relatively strong inhomogeneities or gravitational strengths Biben et al. (1993). This was further confirmed by comparing LDA against the exactly soluble hard rod model in one spatial dimension. While the LDA yields a monotonic decaying density profile $`\rho (z)`$, a layering shows up near the hard wall of the container bottom. Even crystallization can be induced by the bottom wall Biben et al. (1994). As shown recently Heni and Löwen (2000), details of this surface-induced crystallization may be significantly influenced by a periodic wall pattern. Indeed, pure colloidal crystals can be grown from sedimentation on a patterned substrate van Blaaderen et al. (1997); van Blaaderen (1997); van Blaaderen and Wiltzius (1997). In this case the gravitational field acts as an external force enforcing and accelerating heterogenous nucleation and growth. Other fascinating phenomena in a gravitational field relevant for colloidal suspensions are phase transitions such as wetting de Gennes (1981), surface melting Löwen and Beier (1990), as well as dynamical effects as shock like fronts Mills and Snabre (1994), metastable phase formation Bocquet and Löwen (1994), long-range velocity correlations Segré et al. (1997), stratification Mueth et al. (1996), and crystal growth Ackerson et al. (1999).

While the equilibrium sedimentation of hard sphere suspensions is well-understood Biben et al. (1993, 1994); Piazza et al. (1993); Hong (1999); Levin (2000), charged suspensions are much more subtle as they reveal an apparent mass which is smaller than the bare mass at least for intermediate heights Piazza et al. (1993); Biben and Hansen (1994); Löwen (1998); Simonin (1995). In this paper, we study a third kind of effective interaction between colloids, namely a very soft core as realized for star polymer solutions Likos (2001). The qualitative new feature of those solutions as compared to the traditional hard-sphere and Yukawa interactions is that their phase diagram exhibits a reentrant melting behavior for increasing density Watzlawek et al. (1999). In fact, our analysis holds for any system with a reentrant melting behavior but we will mainly focus explicitly on star polymers. Star polymers consist of $`f`$ linear polymer arms attached to a central common core. The complete bulk phase diagram for star polymers in a good solvent was calculated in Watzlawek et al. (1999) and exhibits several unusual solid lattices as well as reentrant melting. As will be discussed in detail in the following sections, due to the reentrant melting behavior, unusual density profiles, featuring interesting effects, arise and a wealth of scaling laws can be established.

The paper is organized as follows: In Sec. II results of computer simulations of a system of star polymers, interacting by means of an ultrasoft pair potential Likos et al. (1998) are presented. In Sec. III, we present a phenomenological theory giving account of the sedimentation profiles observed in the computer simulations. Scaling laws are put forth. Also in Sec. III, density functional theory in a simplified hybrid weighted density approximation (HWDA) is used to reproduce density oscillations at the fluid-gas interface found in the simulation data. Concluding remarks are contained in Sec. IV.

## II Computer Simulation

We performed canonical MC computer simulations keeping particle number $`N`$, volume $`V`$, and temperature $`T`$ constant. We used a simulation box with squared periodic boundary conditions in $`x`$, $`y`$-direction and semi-infinite geometry in $`z`$-direction where the particles were confined only by the gravitational field for $`z>0`$. The bottom wall at $`z=0`$ was hard and interacting with the star polymers by means of an effective star polymer-wall potential which is derived from the effective star polymer-hard sphere interaction in the limit of a sphere with zero curvature. The calculation was performed in Jusufi et al. (2000). It is of the following form:

$`\beta V_{\mathrm{sw}}(z)`$ $`=`$ $`\mathrm{\Lambda }f^{3/2}`$
$`\times `$ $`\{\begin{array}{cc}\mathrm{}\hfill & z<0\hfill \\ \xi _2\mathrm{ln}(\frac{2z}{\sigma })(\frac{4z^2}{\sigma ^2}1)(\xi _1\frac{1}{2})\hfill & z<\frac{\sigma }{2}\hfill \\ \xi _2(1\text{erf}(2\kappa z))/(1\text{erf}(\kappa \sigma ))\hfill & \text{else.}\hfill \end{array}`$
With $`z`$ we denote the distance from the center of one star polymer to the surface of the flat wall. $`\sigma `$ defines the so-called corona diameter of a star polymer, which is related to its diameter of gyration $`\sigma _\mathrm{g}`$ through $`\sigma 0.66\sigma _\mathrm{g}`$, see Jusufi et al. (2000). The constants are $`\mathrm{\Lambda }=0.24`$, $`\kappa \sigma =0.84`$, $`\xi _1=1/(1+2\kappa ^2\sigma ^2)`$, $`\xi _2=\frac{\sqrt{\pi }\xi _1}{\kappa \sigma }\mathrm{exp}(\kappa ^2\sigma ^2)(1\text{erf}(\kappa \sigma ))`$ and the inverse thermal energy $`\beta =1/k_BT`$. We emphasize that the range of the star-wall interaction is of the order of one or two corona diameters, so that the behavior of the sedimentation profiles for larger distances is not influenced. The star polymer pair potential is ultrasoft and is described by the following equation Likos et al. (1998):

$`\beta V_{\mathrm{ss}}(r)`$ $`=`$ $`{\displaystyle \frac{5}{18}}f^{3/2}\{\begin{array}{cc}\mathrm{ln}(\frac{r}{\sigma })+\frac{1}{1+\sqrt{f}/2}\hfill & r<\sigma \hfill \\ \frac{\sigma /r}{1+\sqrt{f}/2}\mathrm{exp}(\frac{\sqrt{f}}{2\sigma }(r\sigma ))\hfill & \text{else,}\hfill \end{array}`$
with center-to-center distance $`r`$. Both interactions are purely entropic, hence they scale linearly with temperature. Previous work Watzlawek et al. (1999) showed that a system of star polymers interacting by means of the potential (II) possess a very rich and interesting bulk phase diagram, see Fig. 1, exhibiting reentrant melting and reentrant freezing transitions for arm numbers $`f_\mathrm{c}f54`$, with the critical arm number $`f_\mathrm{c}=34`$. As we will discuss in more detail below, it is the reentrant melting that makes this type of system appropriate for the analysis presented in this article. The suspending liquid is assumed to be incompressible. Furthermore we treat the solvent to be continuous, neglecting possible effects of the discreteness of the solvent particles. Given the size of the colloidal particles under observation, the star polymers, this is a reasonable assumption. In the simulation, the initial configuration of the system was chosen to be a body-centered cubic (bcc) solid to facilitate equilibration. Its lattice constant $`a`$ was determined from a bulk system with a packing fraction $`\eta =\frac{\pi }{6}\rho \sigma ^30.5`$ lying in the bcc-regime in the bulk phase diagram, see Fig. 1. The lateral box-dimensions were chosen to be multiples of the lattice constant $`a`$.

The total number of particles was then fixed by prescribing a certain value of the thermodynamic variable $`\tau `$, giving the number density per unit surface. The density profile $`\rho (z)`$ is normalized as

$$\tau =_0^{\mathrm{}}\rho (z)𝑑z.$$
(3)
$`\tau \sigma ^2`$ is the number of particles piled up over the area $`\sigma ^2`$ of the bottom wall. Typical system sizes were $`N=2000`$ particles and the Monte Carlo runs were extended over $`N_{\mathrm{MC}}\mathrm{500\hspace{0.17em}000}`$ cycles, each cycle comprising one trial move for each of the $`N`$ particles. Besides the aforementioned thermodynamic variable $`\tau `$, two further parameters characterize the state of the system: First, the arm number $`f`$ of the star polymers, being the number of polymer chains grafted on the central core. Second, the dimensionless gravitational strength (or Peclet number)

$$\alpha =\frac{mg\sigma }{k_BT},$$
(4)
which describes the ratio of the potential energy gain to the thermal energy $`k_BT`$ for a particle of mass $`m`$, displaced by $`\sigma `$ in height in an external field with acceleration $`g`$. The three parameters $`f`$, $`\tau `$ and $`\alpha `$ were varied over a broad range of values. The particles were moved by employing the standard Metropolis algorithm.

In Fig. 2 we show results for different gravitational strengths $`\alpha `$, while $`\tau \sigma ^2=48.87`$ and $`f=39`$ are fixed. The $`f=39`$ star polymer system displays reentrant melting in the fluid $``$ bcc $``$ fluid sequence, as seen along the arrow in Fig. 1. The gravitational field forces the local density $`\sigma ^3\rho (z)`$ to take values that scan the range from $`\sigma ^3\rho (z)=0`$ up to high values, $`\sigma ^3\rho (z)3`$. Thus, the local density ‘crosses through’ the range of the phase diagram where the system displays a bulk bcc phase. It it intuitively expected that the system will then feature a solid regime (for intermediate densities) intercalated between two fluid regimes, at low and high densities. We have found that this is indeed what happens but provided that the gravitational strength does not exceed a critical value $`\alpha ^{}`$, as we discuss below.

Let us start from the case where no solid phase appears. For $`\alpha >\alpha ^{}`$ \[Figs. 2(a) and (b)\], we obtain density profiles $`\rho (z)`$ that show three distinct features: First, there is layering on the wall due to packing effects, typically extending over several layers. As $`z`$ increases a fluid regime with density decaying as a linear function of altitude $`z`$ can be distinguished. At some height ($`z25\sigma `$ in (a)) the density rapidly decays to zero. At this strong inhomogeneity oscillations in density with wavelength $`\sigma `$ can be distinguished in the sedimentation profile, which is smooth elsewhere in the linear regime. The linear dependence of the density profile on $`z`$, can be understood in terms of a local density functional mean-field theory, as will be shown in Sec. III.2; the corresponding results from this theory are shown in Figs. 2(a) and (b) with dotted lines. The density oscillations observed in the simulations were reproducible in the framework of density functional theory using a simplified form of the HWDA, as will be discussed in further detail in Sec. III.3.

By lowering the gravitational strength $`\alpha `$ further, a critical strength $`\alpha ^{}`$ in the range $`16.0<\alpha ^{}<17.0`$ is discovered. Below $`\alpha ^{}`$ the density profiles qualitatively change and exhibit a new feature. Strong density oscillations appear, a clear indication for a crystalline phase. They extend over $`10`$ to $`20`$ star diameters, equivalent to several crystalline layers. The length of the crystal grows, as $`\alpha `$ decreases. A typical simulation snapshot is shown in Fig. 3 next to the corresponding equilibrium density profile. Here, the well-ordered crystal phase in the middle of the simulation box ($`20\sigma z30\sigma `$) is clearly visible.

As an additional check for crystalline order, we calculate the local order parameter, $`\mathrm{\Psi }_4`$, that checks for fourfold symmetry in two dimensions around a given particle. It is defined by

$$\mathrm{\Psi }_4(z)=\left|\frac{1}{4N_\mathrm{l}}\underset{j=1}{\overset{N_\mathrm{l}}{}}\underset{<k>}{}e^{4i\varphi _{jk}}\right|,$$
(5)
where the $`k`$-sum includes the four nearest neighbors of the given particle and the $`j`$-sum extends over $`N_\mathrm{l}`$ particles in the corresponding layer. A layer is defined by a slab of thickness $`\delta 0.2a`$, centered around the given particle at elevation $`z`$, which is motivated by the ‘Lindemann melting rule’, assuming a maximum particle displacement of approximately $`10\%`$ around the equilibrium position in a possible crystal regime. The angular brackets indicate a canonical ensemble average. $`\varphi _{jk}`$ is the polar angle of the interparticle distance vector with respect to a fixed reference frame. For ideal fourfold symmetry, i.e., for a particle contained in a bcc-solid layer, $`\mathrm{\Psi }_4=1`$. Due to thermal motion, small defects of the perfect crystalline symmetry arise and usually values of $`\mathrm{\Psi }_4>0.8`$ Heni and Löwen (2000) are taken to be conclusive evidence for a crystalline phase with fourfold-in-layer-symmetry. As can be seen in Figs. 2(c)-(f) our simulation data do indeed show values up to $`\mathrm{\Psi }_40.95`$ in the region of the density profile $`\rho (z)`$ which we already identified to be solid due to the pronounced density oscillations.

Comparing the interval of the packing fraction in which crystallization occurs to the bulk phase diagram in Fig. 1, we may thus conclude that the intercalated solid regime is a manifestation of the reentrant melting in the bulk phase diagram, mapped onto the $`z`$-axis in a system under gravitational influence. The absence of freezing for strong gravitational fields ($`\alpha >\alpha ^{}`$) can now be at least qualitatively understood: for high values of $`\alpha `$, the density profiles grow too fast as $`z`$ approaches the wall, so that the mapping onto the $`z`$-axis results into a domain which is too narrow to sustain crystalline order. In fact, as we will show in detail in Sec. III.2, a minimal, nonvanishing thickness of the crystalline layer is necessary so that the latter can be stably ‘nested’ between the two fluid phases.

## III Theory

### III.1 Density Functional Theory in Local Density Approximation (LDA)

In order to predict scaling relations characterizing crystallization in sedimentation profiles of star polymer solution we apply density functional theory within the framework of the local density approximation (LDA). The latter is a reliable theoretical tool in cases where the density profile of the system varies slowly with $`z`$, so that it can be considered as staying essentially constant at length scales set by the microscopic natural length of the system ($`\sigma `$ in this case). As can be seen in Figs. 2(a) and (b), this is indeed the case if we discard the strong oscillations close to the wall (the layering effect). As the range of these oscillations is much shorter than the range of the density profile itself, the bulk of the free energy of the system resides in the smooth “ramp-like” part of the density profile and the use of the LDA is justified. Accordingly, we will omit the star-wall potential from our considerations in this subsection and consider only the effects of the external gravitational field $`\mathrm{\Phi }_{\mathrm{ext}}(z)=mgz`$.

We work in the grand canonical ensemble and introduce the chemical potential $`\mu `$ and a variational grand potential per unit area, $`\stackrel{~}{\mathrm{\Sigma }}(T,\mu ;[\rho (z)])`$ which is a functional of the density profile. Introducing the ideal and excess per unit area contributions to the intrinsic Helmholtz free energy of the system, $`F_{\mathrm{id}}[\rho (z)]`$ and $`F_{\mathrm{ex}}[\rho (z)]`$ respectively, we find that in the LDA, the expression for $`\stackrel{~}{\mathrm{\Sigma }}(T,\mu ;[\rho (z)])`$ reads as:

$`\stackrel{~}{\mathrm{\Sigma }}(T,\mu ,[\rho (z)])`$ $`=`$ $`F_{\mathrm{id}}[\rho (z)]+F_{\mathrm{ex}}[\rho (z)]`$
$`+{\displaystyle dz\mathrm{\Phi }_{\mathrm{ext}}(z)\rho (z)}\mu {\displaystyle dz\rho (z)}`$
$`=`$ $`k_BT{\displaystyle _0^{\mathrm{}}}dz\rho (z)\left[\mathrm{ln}(\rho (z)\lambda ^3)1\right]`$
$`+`$ $`{\displaystyle _0^{\mathrm{}}}dz\left[f(\rho (z))+(mgz\mu )\rho (z)\right],`$
where $`\lambda =\sqrt{h^2/2\pi mk_BT}`$ is the thermal de Broglie wavelength and $`f(\rho (z))`$ is the Helmholtz free energy per unit volume of the bulk fluid. The minimization of $`\stackrel{~}{\mathrm{\Sigma }}`$ with respect to $`\rho (z)`$ yields the equilibrium profile $`\rho _0(z)`$; the value of the functional at equilibrium, $`\stackrel{~}{\mathrm{\Sigma }}(T,\mu ,[\rho _0(z)])`$ is then the Gibbs free energy per unit area, $`\mathrm{\Sigma }(T,\mu )`$ of the system. Setting $`\delta \stackrel{~}{\mathrm{\Sigma }}(T,\mu ,[\rho (z)])/\delta \rho (z)|_{\rho _0(z)}=0`$ in Eq. (LABEL:dft.lda), leads to:

$$k_BT\mathrm{ln}(\rho _0(z)\sigma ^3)+f^{}(\rho _0(z))=\mu ^{}mgz,$$
(7)
where $`f^{}(x)`$ denotes the derivative of $`f(x)`$ and $`\mu ^{}=\mu 3\mathrm{ln}(\lambda /\sigma )`$ is a shifted chemical potential.

Due to the ultrasoft character of the logarithmic-Yukawa star-star interaction $`V_{\mathrm{ss}}(r)`$, the star polymer system belongs to the class of mean-field fluids Watzlawek et al. (1998); Lang et al. (2000); Likos et al. (2001); Louis et al. (2000); Louis (2001), for which the excess free energy density is a quadratic function of $`\rho `$, namely:

$$f(\rho )\frac{\rho ^2}{2}4\pi _0^{\mathrm{}}\mathrm{d}rr^2V_{\mathrm{ss}}(r)=:\frac{\widehat{V}_{\mathrm{ss}}(0)\rho ^2}{2},$$
(8)
with the Fourier transform $`\widehat{V}_{\mathrm{ss}}(k)`$ of the pair potential. This property is valid for high density fluids provided their pair potential $`V(r)`$ is only slowly diverging at the origin and decays fast enough to zero as $`r\mathrm{}`$, so that it is integrable. For the more restrictive case of a nondiverging potential at $`r=0`$, the stronger condition $`c(r)=\beta V(r)`$ holds approximately Lang et al. (2000); Likos et al. (2001); Louis et al. (2000); Louis (2001), with $`c(r)`$ denoting the direct correlation function of the fluid Hansen and McDonald (1986). This gives again rise to Eq. (8) above through the compressibility equation of state Hansen and McDonald (1986).

Using the dimensionless variables $`xz/\sigma `$, $`\overline{\rho }(x)\rho (z)\sigma ^3`$, $`B\beta \widehat{V}_{\mathrm{ss}}(0)/\sigma ^3`$ and $`\overline{\mu }\beta \mu ^{}`$ and introducing Eq. (8) into Eq. (7), we obtain the equilibrium profile through the equation:

$$\mathrm{ln}(\overline{\rho }_0(x))+B\overline{\rho }_0(x)=\overline{\mu }\alpha x.$$
(9)
For star functionality $`f=39`$ we obtain $`B250`$ and, for $`f=32`$, $`B204`$ \[see Eq. (II)\]. Hence, the second term in the lhs of Eq. (9) above dominates over the logarithmic term for densities $`\overline{\rho }(x)0.10`$. As almost the entire simulation density profile fulfills this condition, we finally omit the logarithmic term from Eq. (9) above and obtain thereby a linear density profile:

$$\overline{\rho }_0(x)=\{\begin{array}{cc}0\hfill & \text{for }x<0,\hfill \\ \frac{\overline{\mu }\alpha x}{B}\hfill & \text{for }0<x<\overline{\mu }/\alpha ,\hfill \\ 0\hfill & \text{for }\overline{\mu }/\alpha <x.\hfill \end{array}$$
(10)
The chemical potential $`\overline{\mu }`$ is now determined through the normalization condition $`_0^{\overline{\mu }/\alpha }dx\overline{\rho }_0(x)=\tau \sigma ^2\overline{\tau }`$, yielding:

$$\overline{\mu }=\sqrt{2\alpha B\overline{\tau }},$$
(11)
and from Eq. (10) the final expression for the density profile:

$$\overline{\rho }_0(x)=\{\begin{array}{cc}0\hfill & \text{for }x<0,\hfill \\ \sqrt{\frac{2\alpha \overline{\tau }}{B}}\frac{\alpha }{B}x\hfill & \text{for }0<x<\sqrt{\frac{2B\overline{\tau }}{\alpha }},\hfill \\ 0\hfill & \text{for }\sqrt{\frac{2B\overline{\tau }}{\alpha }}<x.\hfill \end{array}$$
(12)
The prediction (12) is compared against the MC simulation results in Figs. 2(a) and (b); theory and simulation are in excellent agreement. This linear dependence of the density profile on $`z`$ is the first scaling prediction we make for such systems. Moreover, by introducing Eq. (10) into Eq. (LABEL:dft.lda), and once more ignoring the logarithmic term, we find that the Gibbs free energy per unit area $`\mathrm{\Sigma }(T,\mu )`$ is a power-law of the chemical potential, namely:

$$\beta \sigma ^2\mathrm{\Sigma }(T,\overline{\mu })=\frac{\overline{\mu }^3}{6\alpha B}.$$
(13)
Accordingly, the Helmholtz free energy per unit area, $`\beta \sigma ^2F(T,\overline{\tau })=\beta \sigma ^2\mathrm{\Sigma }(T,\overline{\mu })+\overline{\mu }\overline{\tau }`$ obeys the scaling law:

$$\beta \sigma ^2F(T,\overline{\tau })=\frac{2}{3}\sqrt{2\alpha B}\overline{\tau }^{3/2}.$$
(14)
The thermodynamic relation $`\overline{\mu }=(\beta \sigma ^2F)/\overline{\tau }`$ returns Eq. (11).

We now examine whether the density oscillations occurring at high $`z`$-elevations, which are clearly visible in Figs. 2(a), (b), can be obtained in the framework of the full LDA, with the logarithmic term included, Eq. (9). Though the latter is an implicit equation for $`\overline{\rho }(x)`$, we do not need to solve it in order to answer the question at hand. The key observation is that Eq. (9) delivers an explicit functional form for the inverse function:

$$x(\overline{\rho })=\alpha ^1(\mathrm{ln}\overline{\rho }+B\overline{\rho }\overline{\mu }).$$
(15)
If the LDA profile displayed oscillations, then $`\overline{\rho }(x)`$ would go through various maxima and minima and there should be several points $`x_\mathrm{m}`$ where the derivative $`\overline{\rho }^{}(x_\mathrm{m})`$ would vanish, with the implication that the derivative of the inverse function, $`x^{}(\overline{\rho }_\mathrm{m})`$, would diverge at the corresponding density values $`\overline{\rho }_\mathrm{m}`$. From Eq. (15) above, we obtain $`x^{}(\overline{\rho })=\alpha ^1(\overline{\rho }^1+B)<0`$ for all $`0<\overline{\rho }<\mathrm{}`$. The only divergence of $`x^{}(\overline{\rho })`$ occurs for the trivial limit $`\overline{\rho }0`$ and corresponds to the exponential decay $`\overline{\rho }(x)e^{\alpha x}`$, valid for high elevations. The LDA is incapable to reproduce this effect, a feat that, in fact, could have been anticipated: these oscillations occur at length scales $`\sigma `$, whereas the LDA is applicable when the spatial inhomogeneity of the profile has a characteristic length much larger than the latter. In Sec. III.3, we resort to a more powerful density functional approximation in order to reproduce this feature of the density profile.

### III.2 Phenomenological Landau Theory

If one focuses close to the reentrant melting transition point and studies the case for small gravity $`(\alpha 1)`$, a phenomenological Landau-like approach can be adopted to explore further scaling predictions for the crystallization transition. We study the situation sketched in Fig. 4

of a crystalline sheet of width $`l`$ intervening between two fluid parts of the sedimentation profile. Let us define the excess grand canonical free energy per unit area, $`\mathrm{\Sigma }_{\mathrm{ex}}(l)`$, in such a situation with respect to a situation where no crystallization takes place. This quantity is given by:

$$\mathrm{\Sigma }_{\mathrm{ex}}(l)=\mathrm{\Sigma }(l)\mathrm{\Sigma }(l=0),$$
(16)
with $`\mathrm{\Sigma }(l)`$ being the grand canonical free energy per unit area when a solid of thickness $`l`$ is present. Evidently, $`\mathrm{\Sigma }(l=0)`$ is the quantity given by Eq. (13) above.

The excess grand canonical free energy $`\mathrm{\Sigma }_{\mathrm{ex}}`$ per unit area comprises of three parts:

1. The equilibrium surface tensions $`\gamma _1`$ and $`\gamma _2`$: These describe the additional free energy cost in creating the two solid-fluid interfaces at $`z=z_1`$ and $`z=z_2`$.
2. A thermodynamic contribution which essentially depends on the arm number. If $`f>f_\mathrm{c}`$ this contribution favors a solid sheet.
3. A free energy penalty due to an elastic distortion of the solid in the external field.

Hence:

$$\mathrm{\Sigma }_{\mathrm{ex}}=\gamma _1+\gamma _2+\mathrm{\Sigma }_{\mathrm{TD}}+\mathrm{\Sigma }_{\mathrm{elast}}.$$
(17)
Let us discuss the different contributions in more detail:

The surface tension will mainly control the relative orientation of the solid with respect to the $`z`$-direction. One expects that a close-packed surface of the bcc-solid (i.e. a (100) orientation) will have smallest surface tension and will hence be the realized orientation. In fact, this is what we found in our simulation data presented in Figures 2 and 3. For hard-sphere fcc-solids the interfacial fluid-solid free energy has been calculated recently in equilibrium by computer simulation Davidchak and Laird (2000). Its order of magnitude is

$$\gamma _i\frac{k_BT}{\sigma ^2},(i=1,2),$$
(18)
where $`\sigma `$ is a microscopic length scale.

The thermodynamic contribution could be calculated within the LDA with $`f(\rho ,T)`$ taken from liquid state theories for the fluid and solid cell theory for the crystal. Here, we will simply focus on a Landau-type theory close to the reentrant melting point characterized by a critical star number density $`\rho _\mathrm{c}=\frac{6}{\pi \sigma ^3}\eta _\mathrm{c}`$ and the critical arm number $`f_\mathrm{c}`$, see Fig. 1. Performing a Landau expansion and dropping the temperature dependence one gets

$`f_\mathrm{s}(\rho )`$ $`=`$ $`f_\mathrm{s}(\rho _\mathrm{c})+A_\mathrm{s}(f_\mathrm{c}f)`$ (19)
$`+f_\mathrm{s}^{}(\rho _\mathrm{c})(\rho \rho _\mathrm{c})`$
$`+{\displaystyle \frac{1}{2}}f_\mathrm{s}^{\prime \prime }(\rho _\mathrm{c})(\rho \rho _\mathrm{c})^2+\mathrm{}.`$
Here, $`f_\mathrm{s}(\rho )`$ is the free energy per unit volume of the solid phase and $`A_\mathrm{s}`$ is a constant governing the first leading term in an expansion around $`f=f_\mathrm{c}`$. Likewise in the fluid phase one has

$`f_\mathrm{f}(\rho )`$ $`=`$ $`f_\mathrm{f}(\rho _\mathrm{c})+A_\mathrm{f}(f_\mathrm{c}f)`$ (20)
$`+f_\mathrm{f}^{}(\rho _\mathrm{c})(\rho \rho _\mathrm{c})`$
$`+{\displaystyle \frac{1}{2}}f_\mathrm{f}^{\prime \prime }(\rho _\mathrm{c})(\rho \rho _\mathrm{c})^2+\mathrm{},`$
with $`f_\mathrm{s}(\rho _\mathrm{c})=f_\mathrm{f}(\rho _\mathrm{c})`$, $`f_\mathrm{s}^{}(\rho _\mathrm{c})=f_\mathrm{f}^{}(\rho _\mathrm{c})`$, but $`f_\mathrm{s}^{\prime \prime }(\rho _\mathrm{c})>f_\mathrm{f}^{\prime \prime }(\rho _\mathrm{c})`$ in general. Performing the inversion of $`f(\rho )`$ in order to get the density profile leads to a piecewise linear profile for the averaged density with two density jumps at $`z=z_1`$ and $`z=z_2`$, as determined by the Maxwell construction, see Fig. 4.

$$\rho (z)=\{\begin{array}{cc}\frac{\mu f_\mathrm{s}^{}(\rho _\mathrm{c})mgz}{f_\mathrm{s}^{\prime \prime }(\rho _\mathrm{c})}+\rho _\mathrm{c}\hfill & \text{for }z_1<z<z_2\hfill \\ \frac{\mu f_\mathrm{f}^{}(\rho _\mathrm{c})mgz}{f_\mathrm{f}^{\prime \prime }(\rho _\mathrm{c})}+\rho _\mathrm{c}\hfill & \text{else.}\hfill \end{array}$$
(21)
Consequently, by inserting this into the free energy function one gets

$$\mathrm{\Sigma }_{\mathrm{TD}}=a(ff_\mathrm{c})l+\left(\frac{1}{f_\mathrm{f}^{\prime \prime }(\rho _\mathrm{c})}\frac{1}{f_\mathrm{s}^{\prime \prime }(\rho _\mathrm{c})}\right)\frac{m^2g^2l^3}{12}.$$
(22)
Note that $`a>0`$ in order to stabilize the solid for $`f>f_\mathrm{c}`$.

Third, the elastic part can be calculated by continuum elastic distortion theory of the solid. For a different situation of a solid in an external field this has been formulated by Gittes and Schick Gittes and Schick (1984). Following these ideas, we assume a $`z`$-independent lateral strain $`ϵ_{}`$ but consider a $`z`$-dependent vertical strain $`ϵ_{}`$. By symmetry, $`ϵ_{}`$ has to be zero for the crystal being stable at $`z=z_0`$, i.e., for the crystal at the reentrant melting point. Elasticity theory predicts for $`\mathrm{\Sigma }_{\mathrm{elast}}`$

$$\mathrm{\Sigma }_{\mathrm{elast}}\frac{1}{2}_{\frac{l}{2}}^{\frac{l}{2}}dzCϵ_{}^2(z),$$
(23)
where $`C>0`$ is related to the elastic constants of the solid. As

$$ϵ_{}\rho \rho _\mathrm{c}mgz,$$
(24)
we obtain

$$\mathrm{\Sigma }_{\mathrm{elast}}=C^{}l^3m^2g^2,$$
(25)
with another constant $`C^{}`$. Eq. (25) has a similar form as the second term of Eq. (22).

In summary, the total grand canonical excess free energy is

$$\mathrm{\Sigma }_{\mathrm{ex}}(l)=a(ff_\mathrm{c})l+b\alpha ^2l^3+\gamma _1+\gamma _2,$$
(26)
where $`a,b>0`$. Moreover we will use $`\gamma =\gamma _1+\gamma _2`$ from now on. A first order phase transition takes place if the minimum of $`\mathrm{\Sigma }_{\mathrm{ex}}`$ with respect to $`l`$ yields $`\mathrm{\Sigma }_{\mathrm{ex}}=0`$.

This determines the following resulting scaling relations:

1. The realized crystalline thickness $`l`$ as obtained by minimizing $`\mathrm{\Sigma }_{\mathrm{ex}}`$ with respect to $`l`$ for fixed $`\alpha `$ and $`f`$ scales as

$$l\frac{\sqrt{ff_\mathrm{c}}}{\alpha }.$$
(27)
2. The phase transition to a crystalline sheet is first order. It happens beyond a critical $`\alpha `$-dependent arm number $`f_{\mathrm{crit}}`$ where

$$f_{\mathrm{crit}}f_\mathrm{c}\gamma ^{2/3}\alpha ^{2/3},$$
(28)
with a scaling exponent of $`\frac{2}{3}`$.
3. The width $`l_0`$ corresponding to the transition scales as

$$l_0\gamma ^{1/3}\alpha ^{2/3}\frac{\gamma }{ff_\mathrm{c}}.$$
(29)

The analysis presented here is general, the scaling predictions derived are valid for any reentrant melting behavior in equilibrium (e.g., laser-induced freezing Chakrabarti et al. (1995); Frey et al. (1999) or polydisperse systems Bartlett and Warren (1999)). Furthermore all these relations can, in principle, be checked by simulation. Relations (28) and (29), however, require high computational efforts. In order to check on scaling relation (27), we measured the crystal length $`l`$ in MC simulations varying $`\alpha `$ or $`f`$, while keeping the density $`\tau `$ fixed. The crystal length is determined by the range $`\mathrm{\Delta }z`$, where the order parameter $`\mathrm{\Psi }_4(z)`$ has values larger than $`0.8`$. The results are plotted in Fig. 5, showing excellent agreement with the scaling predictions.

### III.3 Weighted Density Approximation of the Density Functional

In order to verify the density oscillations close to the fluid-gas interface of the sedimentation profiles for large values of $`\alpha `$, we apply a simplified form of the HWDA (hybrid weighted density approximation). The full HWDA was constructed by Leidl and Wagner in Leidl and Wagner (1993). Given an external field $`\mathrm{\Phi }_{\mathrm{ext}}(z)`$ the free energy is a unique functional of the density profile $`\rho (z)`$. Thus, the excess free energy per unit surface in the HWDA framework is given by

$`_{\mathrm{exc}}[\rho ]={\displaystyle _0^{\mathrm{}}}\rho (z)f_0(\overline{\rho }(z))dz,`$ (30)
where $`f_0(\rho )`$ denotes the excess free energy per particle of a homogeneous liquid of density $`\rho `$. The weighted density $`\overline{\rho }(z)`$ follows from a convolution with the weighting function $`\omega (r;\rho )`$

$`\overline{\rho }(z)={\displaystyle \rho (z^{})\omega (|𝐫𝐫^{}|;\widehat{\rho })d𝐫^{}}`$ (31)
with a global density $`\widehat{\rho }`$. The weighting function $`\omega (r;\rho )`$ is fixed by a simple quadratic equation in Fourier space Leidl and Wagner (1993):

$`2f_0^{}(\rho _0)\stackrel{~}{\omega }(k;\rho _0)`$ $`+`$ $`\rho _0f_0^{\prime \prime }(\rho _0)\stackrel{~}{\omega }^2(k;\rho _0)`$ (32)
$`=`$ $`\beta ^1\stackrel{~}{c}^{(2)}(k;\rho _0)`$
The primes denote differentiations with respect to the density $`\rho `$ and $`\stackrel{~}{c}^{(2)}(k;\rho _0)`$ is the Fourier transform of the direct correlation function of the homogeneous fluid. A unique solution of $`\stackrel{~}{\omega }(k;\rho _0)`$ is determined by the normalization $`\stackrel{~}{\omega }(k=0;\rho _0)=1`$, also ensuring the compressibility rule to be satisfied. We have solved the homogeneous problem with Ornstein-Zernike fluid integral equations using the Rogers-Young closure Rogers and Young (1984). Resulting correlation functions and structure factors are in very good agreement with MC simulations of the bulk system Watzlawek et al. (1998). In difference to the complete HWDA, where the global density $`\widehat{\rho }`$ is chosen to be a functional of $`\rho (z)`$, we keep $`\widehat{\rho }`$ fixed. This simplification is sufficient to verify the observed oscillations, accompanied by the advantage that the numerical effort is enormously reduced. Best agreement with simulation results could be achieved when choosing the global density $`\widehat{\rho }`$ to be of the order of the averaged density near the bottom wall $`z=0`$. The tails of the density profiles are nearly unaffected by the choice of $`\widehat{\rho }`$. A similar approach is used in the SWDA Kim and Suh (1996) for an inhomogeneous fluid in contact with a bulk fluid of density $`\rho _b`$; there $`\widehat{\rho }`$ was chosen to be $`\rho _b`$.

Applying the usual Euler-Lagrange minimization for the Helmholtz free energy per unit area $`[\rho ]`$ with chemical potential $`\mu `$

$`{\displaystyle \frac{\delta }{\delta \rho (z)}}=\mu \mathrm{\Phi }_{\mathrm{ext}}(z),`$ (33)
and using $`\mathrm{\Phi }_{\mathrm{ext}}(z)=V_{\mathrm{sw}}(z)+\alpha z/\beta \sigma `$, we obtain for the density profile $`\rho (z)`$

$`\rho (z)=\{\begin{array}{cc}\xi \mathrm{exp}\{c^{(1)}(z;[\rho ])\alpha z/\sigma \beta V_{\mathrm{sw}}(z)\}\hfill & z>0\hfill \\ & \\ 0\hfill & \text{else}.\hfill \end{array}`$ (34)
The fugacity $`\xi `$ is determined by the normalization condition $`\xi =\tau /_0^{\mathrm{}}dz\mathrm{exp}\{c^{(1)}(z;[\rho ])\alpha z/\sigma \beta V_{\mathrm{sw}}(z)\}`$. $`c^{(1)}(z;[\rho ])`$ is the one-particle correlation function:

$`\beta ^1c^{(1)}(z;[\rho ])`$ $`=`$ $`{\displaystyle \frac{\delta _{\mathrm{ex}}[\rho ]}{\delta \rho (z)}}`$
$`=`$ $`f_0(\overline{\rho }(z))`$
$`+{\displaystyle dz^{}\rho (z^{})f_0^{}(\overline{\rho }(z^{}))\omega (|zz^{}|;\widehat{\rho })}.`$
Equation (34) was solved for the profile by standard iterative techniques, see, e.g., Ref. Denton and Ashcroft (1991). The results for an arm number $`f=32`$, $`\tau \sigma ^2=21.8`$ and three different values of $`\alpha `$ are shown in Fig. 6 together with MC simulation data.

The global density $`\widehat{\rho }`$ for all three profiles is fixed at $`\widehat{\rho }\sigma ^3=1.8`$. The DFT results are in very good agreement with the simulation profiles. In particular, the interface oscillations with wavelength $`\sigma `$ also occur in the DFT. For $`\alpha =10.0`$, the lowest $`\alpha `$ that is shown, the profile is nearly indistinguishable from MC data far from the wall, while for increased $`\alpha `$, $`\alpha =20.0`$ and $`\alpha =30.0`$ the interface oscillations are underestimated in our DFT approach.

These oscillations are not a specific feature of the star polymers; we have also performed MC simulations using a repulsive Yukawa interaction of the form

$`V(r)\mathrm{exp}(\kappa r)/r,`$ (36)
corresponding to the part of the star polymer pair potential (II) valid for distances $`r\sigma `$. Here exactly the same behavior could be found at the tails of the density profiles. For hard sphere systems, on the other hand, no such density oscillations are present. This suprising fact might be attributed to the long ranged tail in the interaction potential.

## IV Conclusions

Concluding, we have presented results for systems exhibiting reentrant melting in the bulk phase diagram, under gravitational influence. It was shown that a phase transition occurs when the gravitational strength $`\alpha `$ is varied: Below a critical $`\alpha ^{}(f,\tau )`$, intercalated crystallization occurs in the sedimentation profiles of the observed star-polymer solutions, whereas for $`\alpha >\alpha ^{}(f,\tau )`$ we find monotonic sedimentation profiles $`\rho (z)`$. In MC computer simulations scaling relations for the crystallization, predicted in the framework of a phenomenological theory, valid for all systems exhibiting reentrant melting in the bulk phase diagram, could be verified. Using density functional theory, density oscillations at the fluid-gas boundary, observed in the MC simulations, could be reproduced.

In principle, our results can be verified in surface-sensitive scattering experiments or real-space imaging methods for colloidal suspensions. Unlike non-monotonicities on the liquid side of the gas-liquid equilibrium interface (see e.g. Ohnesorge et al. (1994)), the density oscillations on top of the sedimentation profile are not affected by capillary wave fluctuations and may thus be verified in real samples. The intervening solid sheet should be signalled by a Bragg-like peak in surface reflection measurements.

###### Acknowledgements.

The authors wish to thank M. Heni, M. Schmidt, and A. Esztermann for stimulating discussions and A. Esztermann for a critical reading of the manuscript. This work has been supported by the Deutsche Forschungsgemeinschaft within the SFB 237.