text:
# Sensor Alignment by Tracks

## I INTRODUCTION

For full exploitation of high resolution position sensitive detectors, it is crucial to determine the detector location and orientation to a precision better than their intrinsic resolution. It is a very demanding task to assemble a large number of detector units in a large and complex detector system to this high precision. Also, after assembly, the position determination of the modules by optical survey has its limitations because of detectors obscuring each other. Therefore the final tuning of detector and sensor positions is made by using reconstructed tracks.

In this paper we present an effective method by which individual sensors in a detector setup can be aligned to a high precision with respect to each other. The basic idea is illustrated in Figure 1. Using a large number of tracks, an optimum of each sensor position and orientation is determined such that the track fit residuals are minimized.

The outline of this paper is as follows: In Section II we briefly review published alignment methods. In Section III we introduce the basic notations and coordinate systems involved in our method. In Section IV we present the detailed formulation of the method. In Sections V and VI we demonstrate the performance of the method applied to a test beam setup and to a simulated pixel vertex detector, respectively. The CMS CMS-tracker  Pixel detector is is used as a model in the simulation.

## II BRIEF REVIEW OF ALIGNMENT METHODS

Most HEP experiments equipped with precise tracking detectors have to deal with misalignment issues, and several different approaches for alignment by tracks have been used and reported. Most methods are iterative with 5-6 parameters solved at a time.

Several papers concerning different aspects of alignment in the DELPHI experiment can be found in the literature. For instance, $`Z^0\mu ^+\mu ^{}`$ and cosmic rays are used for the global alignment between sub-detectors VD, OD and TPC litDelphiGlob . The most detailed DELPHI alignment paper deals with the alignment of the Microvertex detector litDelphiPed .

In the ALEPH experiment, alignment was carried out wafer by wafer, and with 20 iterations and 20000 $`Z^0q\overline{q}`$ and 4000 $`Z^0\mu ^+\mu ^{}`$ events an accuracy of a few $`\mu m`$ can be achieved litAleph .

A different, computationally challenging approach is chosen in the SLD experiment, where the algorithm requires simultaneous solution of 576 parameters leading to a 576 by 576 matrix inversion litSLD576 . In the SLD vertex detector, a recently developed matrix singular value decomposition technique is also used for internal alignment litSLDsvd .

## III COORDINATE SYSTEMS AND TRANSFORMATIONS

Our method is applicable to detector setups which consist of planar sensors like silicon pixel or strip detectors. For track reconstruction one conventionally uses the local (sensor) coordinate system and the global detector system. The local system $`(u,v,w)`$ is defined with respect to a detector module (sensor) as follows: The origin is at the center of the sensor, the $`w`$-axis is normal to the sensor, the $`u`$-axis is along the precise coordinate and the $`v`$-axis along the coarse coordinate. The global coordinates are denoted as $`(x,y,z)`$.

The transformation from the global to the local system goes as:

$$𝐪=𝐑(𝐫𝐫_0)$$
(1)
where $`𝐫=(x,y,z)`$, $`𝐪=(u,v,w)`$, $`𝐑`$ is a rotation and $`𝐫_0=(x_0,y_0,z_0)`$ is the position of the detector center in global coordinates.

In the very beginning of the experiment the rotation $`𝐑`$ and the position $`𝐫_0`$ are determined by detector assembly and survey information. In the course of experiment this information will be corrected by an incremental rotation $`\mathrm{\Delta }𝐑`$ and translation $`\mathrm{\Delta }𝐫`$ so that the new rotation and translation become:

$`𝐑`$ $``$ $`\mathrm{\Delta }\mathrm{𝐑𝐑}`$ (2)
$`𝐫_0`$ $``$ $`𝐫_0+\mathrm{\Delta }𝐫.`$ (3)
The correction matrix $`\mathrm{\Delta }𝐑`$ is expressed as:

$$\mathrm{\Delta }𝐑=𝐑_\gamma 𝐑_\beta 𝐑_\alpha $$
(4)
where $`𝐑_\alpha ,𝐑_\beta `$ and $`𝐑_\gamma `$ are small rotations by $`\mathrm{\Delta }\alpha ,\mathrm{\Delta }\beta ,\mathrm{\Delta }\gamma `$ around the $`u`$-axis, the (new) $`v`$-axis and the (new) $`w`$-axis, respectively. The position correction $`\mathrm{\Delta }𝐫`$ transforms to the local system as:

$$\mathrm{\Delta }𝐪=\mathrm{\Delta }\mathrm{𝐑𝐑}\mathrm{\Delta }𝐫$$
(5)
with $`\mathrm{\Delta }𝐪=(\mathrm{\Delta }u,\mathrm{\Delta }v,\mathrm{\Delta }w)`$. Using (1-5) we find the corrected transformation from global to local system as:

$$𝐪^c=\mathrm{\Delta }\mathrm{𝐑𝐑}(𝐫𝐫_0)\mathrm{\Delta }𝐪.$$
(6)
where the superscript $`c`$ stands for ’corrected’. The task of the alignment procedure by tracks is to determine the corrective rotation $`\mathrm{\Delta }𝐑`$ and translation $`\mathrm{\Delta }𝐫`$ or $`\mathrm{\Delta }𝐪`$ as precisely as possible for each individual detector element.

## IV DESCRIPTION OF THE ALIGNMENT ALGORITHM

### IV.1 Basic Formulation

Since the alignment corrections are small, the fitted trajectories can be approximated with a straight line in a vicinity of the detector plane. The size of this small region is determined by the alignment uncertainty which is expected to be at most a few hundred microns so that the straight line approximation is perfectly valid.

The equation of a straight line in global coordinates, approximating the trajectory in a vicinity of the detector, can be written as:

$$𝐫_s(h)=𝐫_x+h\widehat{𝐬}$$
(7)
where $`𝐫_x`$ is the trajectory impact point on the detector in question, $`\widehat{𝐬}`$ is a unit vector parallel to the line and $`h`$ is a parameter. Equation (7) is for uncorrected detector positions.

Using Eq. (6) the corrected straight line equation in the local system reads:

$$𝐪_s(h)=𝐑_c(𝐫_x+h\widehat{𝐬}𝐫_0)\mathrm{\Delta }𝐪$$
(8)
where $`𝐑_c=\mathrm{\Delta }\mathrm{𝐑𝐑}`$. A point $`𝐪_s=𝐪_s(h_x)`$ which lies in the detector plane must fulfill the condition $`𝐪_s\widehat{𝐰}=0`$, where $`\widehat{𝐰}=(0,0,1)`$ is normal to the detector. From this condition we can solve the parameter $`h_x`$ which gives the corrected impact or x-ing point on the detector:

$$h_x=\frac{[\mathrm{\Delta }𝐪𝐑_c(𝐫_x𝐫_0)]\widehat{𝐰}}{𝐑_c\widehat{𝐬}\widehat{𝐰}}.$$
(9)
The corrected impact point coordinates $`𝐪_x^c`$ in the local system are then:

$$𝐪_x^c=𝐑_c(𝐫_x𝐫_0)+\frac{[\mathrm{\Delta }𝐪𝐑_c(𝐫_x𝐫_0)]\widehat{𝐰}}{𝐑_c\widehat{𝐬}\widehat{𝐰}}𝐑_c\widehat{𝐬}\mathrm{\Delta }𝐪.$$
(10)
Since the uncorrected impact point is $`𝐪_x=𝐑(𝐫_x𝐫_0)`$, Eq. (10) can be written as:

$$𝐪_x^c=\mathrm{\Delta }𝐑𝐪_x+\frac{(\mathrm{\Delta }𝐪\mathrm{\Delta }𝐑𝐪_x)\widehat{𝐰}}{\mathrm{\Delta }𝐑\widehat{𝐭}\widehat{𝐰}}\mathrm{\Delta }𝐑\widehat{𝐭}\mathrm{\Delta }𝐪.$$
(11)
where $`\widehat{𝐭}=𝐑\widehat{𝐬}`$ is the uncorrected trajectory direction in the detectors local frame of reference. Eq. (11) evaluates to:

$$𝐪_x^c=\mathrm{\Delta }𝐑𝐪_x+(\mathrm{\Delta }w[\mathrm{\Delta }𝐑𝐪_x]_3)\frac{\mathrm{\Delta }𝐑\widehat{𝐭}}{[\mathrm{\Delta }𝐑\widehat{𝐭}]_3}\mathrm{\Delta }𝐪.$$
(12)
This expression provides us with a ’handle’ by which the unknowns $`\mathrm{\Delta }𝐪`$ and $`\mathrm{\Delta }𝐑`$ can be estimated by minimizing a respective $`\chi ^2`$ function using a large number of tracks.

### IV.2 General $`\chi ^2`$ Solution

We denote a measured point in local coordinates as $`𝐪_m=(u_m,v_m,0)`$. The corresponding trajectory impact point is $`𝐪_x^c=(u_x,v_x,0)`$. For simplicity we omit the superscripts $`c`$ in the coordinates $`u_x`$ and $`v_x`$. In stereo and pixel detectors we have two measurements, $`u_m`$ and $`v_m`$, and in non-stereo strip detectors only one, $`u_m`$. In the latter case the coarse coordinate $`v_m`$ is redundant. The residual is either a 2-vector:

$$\epsilon =\left(\begin{array}{c}\epsilon _u\\ \epsilon _v\end{array}\right)=\left(\begin{array}{c}u_xu_m\\ v_xv_m\end{array}\right)$$
(13)
or a scalar $`\epsilon =\epsilon _u=u_xu_m`$. In the following we treat the more general 2-vector case. The scalar case is a straightforward specification of the 2-vector formalism.

The $`\chi ^2`$ function to be minimized for a given detector is:

$$\chi ^2=\underset{j}{}\epsilon _j^T𝐕_j^1\epsilon _j$$
(14)
where the sum is taken over the tracks $`j`$. $`𝐕_j`$ is the covariance matrix of the measurements $`(u_m,v_m)`$ associated with the track $`j`$. The alignment correction coefficients, i.e. the three position parameters $`(\mathrm{\Delta }u,\mathrm{\Delta }v,\mathrm{\Delta }w)`$ and the three orientation parameters $`(\mathrm{\Delta }\alpha ,\mathrm{\Delta }\beta ,\mathrm{\Delta }\gamma )`$ are found iteratively by a general $`\chi ^2`$ minimization procedure. At each step of the iteration one uses the so far best estimate of the alignment parameters in the track fit.

Let us denote these parameters as $`𝐩=(\mathrm{\Delta }u,\mathrm{\Delta }v,\mathrm{\Delta }w,\mathrm{\Delta }\alpha ,\mathrm{\Delta }\beta ,\mathrm{\Delta }\gamma )`$. Then, according to the general $`\chi ^2`$ solution, the iterative correction to $`𝐩`$ has the following expression:

$$\delta 𝐩=\left[\underset{j}{}𝐉_j^T𝐕_j^1𝐉_j\right]^1\left[\underset{j}{}𝐉_j^T𝐕_j^1\epsilon _j\right]$$
(15)
where $`𝐉_j`$ is a Jacobian matrix of $`\epsilon _j(𝐩)`$:

$$𝐉_j=_𝐩\epsilon _j(𝐩).$$
(16)
An adequate starting point for the iteration is a null correction vector $`𝐩`$=0.

In the general case of two measurements $`(u_m,v_m)`$, $`𝐉_j`$ is a $`6\times 2`$ matrix. In case of scalar $`\epsilon `$, for single sided strip detectors, $`𝐉_j`$ is a vector of 5 elements, because $`\mathrm{\Delta }v`$ is redundant and cannot be fitted. It will also be foreseen that only a sub-set of the 6 alignment parameters would be fitted and the others kept fixed. In this case the dimension of the Jacobian matrix reduces accordingly.

The derivatives of the Jacobian matrix can be computed to a good precision in the small correction angle approximation (see below). The elements of the matrix $`𝐉`$ for a given track are then:

$$𝐉=\left(\begin{array}{cc}1& 0\\ 0& 1\\ \mathrm{tan}\psi & \mathrm{tan}\vartheta \\ v_x\mathrm{tan}\psi & v_x\mathrm{tan}\vartheta \\ u_x\mathrm{tan}\psi & u_x\mathrm{tan}\vartheta \\ v_x& u_x\end{array}\right)$$
(17)
The quantities $`\mathrm{tan}\psi `$ and $`\mathrm{tan}\vartheta `$ are defined in the next section.

### IV.3 Linearized Solution with the Tilt Formalism

We call ”tilts” the angle corrections $`x`$ which are small enough to justify the approximations $`\mathrm{cos}x1`$ and $`\mathrm{sin}xx`$. In this approximation the correction matrix $`\mathrm{\Delta }𝐑`$ reads:

$$\mathrm{\Delta }𝐑=\left(\begin{array}{ccc}1& \mathrm{\Delta }\gamma & \mathrm{\Delta }\beta \\ \mathrm{\Delta }\gamma & 1& \mathrm{\Delta }\alpha \\ \mathrm{\Delta }\beta & \mathrm{\Delta }\alpha & 1\end{array}\right)$$
(18)
Using Eq. (18) we linearize Eq. (12) and get the following expressions for the corrections of the impact point coordinates as a function of the alignment correction parameters:

$`\mathrm{\Delta }u_x`$ $`=`$ $`\mathrm{\Delta }u+\delta \mathrm{tan}\psi +\mathrm{\Delta }\gamma v_x`$ (19)
$`\mathrm{\Delta }v_x`$ $`=`$ $`\mathrm{\Delta }v+\delta \mathrm{tan}\vartheta \mathrm{\Delta }\gamma u_x`$ (20)
where $`\delta =\mathrm{\Delta }w+\mathrm{\Delta }\beta u_x+\mathrm{\Delta }\alpha v_x`$. The quantity $`\psi `$ is the angle between the track and the $`vw`$plane and $`\vartheta `$ is the angle between the track and the $`uw`$plane: $`\mathrm{tan}\psi =\widehat{t}_1/\widehat{t}_3`$, $`\mathrm{tan}\vartheta =\widehat{t}_2/\widehat{t}_3`$.

With this approximation the residuals (13) depend linearly on all 6 parameters. Hence the $`\chi ^2`$ minimization problem is linear and can be solved by standard techniques without iteration.

## V ALIGNMENT OF A TEST BEAM SETUP

From Eqs. (19) and (20) we can estimate the contributions of various misalignments to the hit measurement errors. For example the contribution of a misalignment $`\mathrm{\Delta }\alpha `$ around the $`u`$-axis to the $`v`$-coordinate is:

$$\mathrm{\Delta }vv\mathrm{\Delta }\alpha \mathrm{tan}\vartheta .$$
(21)
The error is small near normal incident angles, but grows rapidly as a function of $`\vartheta `$. At $`\vartheta =45^o`$ and near the edge of the sensor ($`v=3`$ cm) the error goes as $`30000\mu \text{m}\mathrm{\Delta }\alpha `$ so that for only $`1\text{mrad}`$ error in $`\mathrm{\Delta }\alpha `$ the systematic error in the $`v`$-coordinate is $`30\mu \text{m}`$.

The silicon detector team of Helsinki Institute of Physics made a precision survey of detector resolution as a function of the angle of incidence of the tracks calib\_paper . The study was made in the CERN H2 particle beam with a setup described in Figure 2. One of the silicon strip detectors was fixed on a rotative support which allowed the tracks to enter between $`0`$ and $`30`$ degrees of incident angle. The angular dispersion of the beam was about 10 mrad and the hits covered the full area of the test detector.

In order to obtain reliable results it was extremely important to calibrate the tilt angle to a very high precision. Our algorithm was used in the alignment calibration. In Table 1 we show the result of the alignment demonstrating the precision obtained by about 3000 beam tracks.

With the precise alignment we have been able to determine the optimal track incident angle which minimizes the detector resolution calib\_paper .

## VI MONTE CARLO SIMULATION

### VI.1 Simulated Detector

A Monte Carlo simulation code was written to test the alignment algorithm. High momentum tracks were simulated and driven through a set of detector planes. The simulated hits were fluctuated randomly to simulate measurement errors. Gaussian multiple scattering was added quadratically using the Highland Highland  approximation. The algorithm involves misalignment of a detector setup in order to simulate a realistic detector.

The experimenters’ imperfect knowledge of the true position of the detector planes is simulated by reconstructing the trajectories in the ideal (not misaligned) detector. This means that in the transformation from local to global coordinate system one uses the ideal positions of the detector planes. The full algorithm in brief is as follows:

1. Creation of an ideal detector setup with no misalignments
2. Creation of a misaligned, realistic detector
3. Generation of the particles and hits in the misaligned detector simulating the real detector
4. Reconstruction of the particle trajectories in the nominal (ideal) detector thus using slightly wrong hit positions. This simulates the realistic situation in which the detector alignment is not yet performed.

For the simulated detector type we choose a vertex detector which is a simplification of the CMS Pixel barrel detector CMS-tracker ; pixel-det  with two layers. The setup is illustrated in Figure 3. There are 144 sensors in layer 1 and 240 sensors in layer 2. The distance of the layer 1 from the beam line is about 4 cm and the layer 2 about 8 cm.

In the simulation we used the following conditions:

1. Misalignment of chosen sensors: The shifts $`\mathrm{\Delta }u,\mathrm{\Delta }v,\mathrm{\Delta }w`$ were chosen at random, each in the range $`\pm 100\mu \text{m}`$ and the tilts $`\mathrm{\Delta }\alpha ,\mathrm{\Delta }\beta ,\mathrm{\Delta }\gamma `$ were chosen at random each in the range $`\pm 20\text{mrad}`$.
2. Beam and vertex constraints: The vertex positions were Gaussian fluctuated around the center of the beam diamond with $`\sigma _x=\sigma _y=20\mu \text{m}`$ and $`\sigma _z=7`$cm and the tracks were fitted with the constraint to start from one point, i.e. from the primary vertex.

In the following we consider two different cases of misaligned detectors:

1. All sensors in layer 2 fixed, all sensors in layer 1 misaligned.
2. Only one sensor in layer 2 fixed, all remaining 383 sensors misaligned.

In case I the total number of fitted parameters is $`6\times 144=864`$ and we used about $`2\times 10^5`$ tracks. The case I appears to be an ’easy’ one with which the algorithm copes very well, as we see below. The second case we call ’extreme’ since the alignment is based on one reference sensor which covers only about 0.26 % of the detector setup area. The total number of fitted parameters in this case was $`6\times 383=2298`$. In the following sectios we show perfomance results of the algorithm in these two cases.

### VI.2 Convergence of the Algorithm

The convergence rate of the alignment procedure as a function of the iteration cycle is shown in Figure 4. It appears that the convergence is fast in the ’easy’ case (the 6 plots on the left) where more than 60 % of the sensors provide the reference. The convergence takes place after a couple of iterations.

In the case where only one sensor is taken as a reference (plots on the right of the figure), the situation is different. It appears that the number of iterations needed varies between 20 and 100 from parameter to parameter. It is also seen that the converged parameter values are somewhat off from the true values, but the precision is reasonable.

### VI.3 Comparison of Fitted and True Parameters

The precision of the fitted parameters in comparison with the true values is shown in Figure 5 on the left for the case I. The correlations are very strong. The typical deviation of the fitted parameters from the true value is less than $`1\mu \text{m}`$ for the offsets and a fraction of a milliradian for the tilts. The precision appears to be better than actually needed in this case, indicating that a smaller statistics would give a satisfactory result.

In case II (the plots on the right of the figure) a good correlation is observed, but the precision is somewhat more modest. For example the error in $`\mathrm{\Delta }w`$ (the shift normal to the sensor plane) is still in most cases below $`10\mu \text{m}`$.

## VII CONCLUSIONS

We have developed a sensor alignment algorithm which is mathematically and computationally simple. It is based on repeated track fitting and residuals optimization by $`\chi ^2`$ minimization. The computation is simple, because the solution involves matrices whose dimension is at most $`6\times 6`$. The method is capable of solving simultaneously all six alignment parameters per sensor for a detector setup with a large number of sensors.

We have successfully applied the method in a precision survey of silicon strip detector resolution as a function of the tracks incident angle. Furthermore, we have demonstrated the performance of the algorithm in case of a simulated two-layer pixel barrel vertex detector. The method performs very well in the case where the outer layer is taken as a reference and all inner sensors are to be aligned. The algorithm performs reasonably well also in the extreme case where only one sensor, representing some 0.26 % of the total area, is taken as a reference for the alignment.

###### Acknowledgements.

The authors wish to thank K. Gabathuler, R. Horisberger and D. Kotlinski for inspiring discussions. Work supported by Ella and Georg Ehrnrooth foundation, Magnus Ehrnrooth foundation, Arvid and Greta Olins fund at Svenska kulturfonden and Graduate School for Particle and Nuclear Physics, Finland.


format: text

text:
# Entropy puzzle in small exploding systems

## Abstract

We use a simple hard-core gas model to study the dynamics of small exploding systems. The system is initially prepared in a thermalized state in a spherical container and then allowed to expand freely into the vacuum. We follow the expansion dynamics by recording the coordinates and velocities of all particles until their last collision points (freeze-out). We have found that the entropy per particle calculated for the ensemble of freeze-out points is very close to the initial value. This is in apparent contradiction with the Joule experiment in which the entropy grows when the gas expands irreversibly into a larger volume.

Entropy puzzle in small exploding systems

J.P. Bondorfa, I.N. Mishustina,b,c, G. Neergaarda

aNiels Bohr Institute, DK-2100 Copenhagen Ø, Denmark
bKurchatov Institute, Russian Research Centre, 123182 Moscow, Russia
cInstitute for Theoretical Physics, Goethe University, D-60054 Frankfurt am Main, Germany

Energetic nucleus-nucleus collisions open a unique possibility to study explosive dynamics of strongly interacting many-body systems in the laboratory. Highly excited matter produced in such collisions expands into vacuum until its constituents decouple (the freeze-out stage). There exist many models for describing this process which range from simple macroscopic to fully microscopic ones. Within thermal and fluid dynamical models it is usually assumed that the matter expansion is isentropic, i. e. proceeds at constant entropy. On the other hand, as well known from statistical physics , only slow reversible processes conserve entropy. It is known from the Joule experiment  that the entropy grows if the state of the system changes too fast. In this letter we examine the entropy conservation hypothesis on the basis of a microscopic model.

For this study we employ a simple gas model where the constituent particles collide like billiard balls and follow the classical Newtonian dynamics. This model, first introduced for simulating heavy-ion collisions in ref. , was recently applied  for investigating de-equilibration dynamics in expanding matter. We consider a gas of identical balls of radius $`r_c`$, which perform classical non-relativistic elastic scatterings at impact parameters $`b<2r_c`$ with conservation of energy and momentum. Rotational degrees of freedom of the balls are ignored. The initial system consists of $`N`$ such particles placed randomly within a sphere of radius $`R`$, rejecting configurations where particles overlap within the hard-core distance. The particle velocities are generated from a Gaussian distribution with variance $`T/m`$ where $`T`$ is interpreted as temperature. Then the particles are allowed to collide for a certain time (“cooking” stage) in order to fully equilibrate the system. For our simple interaction the total energy of the gas is obviously given by the ideal-gas relation $`E=3NT/2`$ independent of gas density.

In our simulations we arbitrarily choose conventional nuclear scales: The mass of the constituent particle is $`m=940`$ $`MeV/c^2`$, the hard core radius is $`r_c=0.5`$ fm, and the initial radius of the gas sphere is $`R=r_0N^{1/3}`$. Most simulations are performed with $`r_0=1.2`$ fm corresponding to the normal nuclear density $`\rho _00.14`$ fm-3. To investigate the role played by the finite-size effects we have performed simulations for 4 systems, N=50, 100, 200 and 400. This covers the range of baryon numbers actually achievable in heavy-ion collisions. The initial average energy $`E_i118`$ MeV per particle and the corresponding temperature $`T_i78`$ MeV were chosen to safely ignore quantum and relativistic effects. The characteristic sound velocity for an ideal gas at this temparature is $`c_s\sqrt{T/m}0,186c`$.

When the gas is confined in a container the particles collide not only with each other but also with the container wall. When the container expands the gas particles loose energy and momentum while colliding with the moving walls. In the case of slow expansion these losses are rapidly redistributed over all particles and the gas remains in thermal equilibrium. This case corresponds to the reversible process when the temperature decreases with volume according to the adiabatic relation

$$TV^{\gamma 1}=\mathrm{const},$$
(1)
where $`\gamma 5/3`$ is the adiabatic index. However, when the expansion is fast fewer particles reach the wall and the energy losses are smaller than needed for the adiabatic expansion. In the case of a very fast expansion of the container no particles can collide with the wall and therefore the energy of the gas remains constant. If the wall stops at a larger radius, the gas will eventually relax to a new equilibrium state in the larger volume. The relaxation time can be estimated as $`\mathrm{\Delta }R/c_s`$. Since the total energy and accordingly the temperature is practically unchanged, the entropy of the equilibrated gas increases due to the larger volume, as expected for a fast irreversible process. This simple physics is behind the Joule experiment. Although the traditional Joule experiment was performed with nonspherical containers, the general principles are obviously valid for the spherical geometry considered in this paper. Our simulations show that the transition from the slow to fast expansion is rather sharp and takes place at wall velocities of approximately 0.5$`c_s`$.

More difficult problems arise when gas or fluid expand into the vacuum without any container. The question which we want to address is whether the expanding matter itself generates a sort of wall effect which may simulate an isentropic process. This question is closely related to the problem of freeze-out and collective flow in exploding systems. The simplest scenario which is often used in the literature is to say that the system expands isentropically until all interactions between the constituents cease. Then the change in the internal energy of the matter is transferred into a collective flow energy. But the problem is that the freely expanding system has no well defined volume. In our previous study  we have defined the instantaneous volume by taking a high, 20th, moment of the particle spatial distribution. For each time step we have defined the entropy as $`S=_kp_k\mathrm{ln}p_k`$ where $`p_k`$ is the occupation probability of the phase space cell $`k`$ in the comoving grid. This entropy was compared with a reference entropy $`S_{ref}`$ defined for the equilibrated system of the same volume. From simulations at different initial conditions for a system of 50 particles we have found that the equilibration measure $`\mathrm{\Sigma }=\mathrm{exp}(SS_{ref})`$ at late times was not equal to 1 but rather close to 0.6.

Below we adopt a slightly different strategy using the microscopic information on the freeze-out field. The initial state is prepared in the same way as before but now after a “cooking” stage the container wall is completely removed and the gas is allowed to expand freely into the vacuum. It is important to stress that this free expansion starts from a state with well defined temperature and density. All particles are followed until their last collision when their coordinates and momenta are recorded. For each system we generate many such events and define the freeze-out field as the set of all such coordinates and momenta. As demonstrated in ref. , these fields are nonlocal in space and time, in contrast to a simplified Cooper-Frye picture  assuming a sharp freeze-out hypersurface. We point out also that the number of freeze-out points per event is generally less than the number of particles because some particles leave the system without any collision (see Table 1 below).

After obtaining the freeze-out field we calculate the average characteristics of the phase space occupation. Utilising spherical symmetry of the system we divide it into a number of spherical shells of radii $`R_k`$. For each shell we calculate the average density of freeze-out points $`\rho (r)`$, collective velocity $`u(r)`$ and temperature $`T(r)`$. The collective velocity is defined simply as the mean radial velocity of frozen-out particles in a given shell i. e. between $`R_k`$ and $`R_{k+1}`$,

$$u(r)\overline{v}(r)=\frac{1}{N_k}\underset{i=1}{\overset{N_k}{}}v_r(𝐫_i),$$
(2)
where $`N_i`$ is the number of freeze-out points in this shell, $`_kN_k=N`$. The temperature is determined from the variance of velocities in the shell, assuming the ideal-gas relation,

$$T(r)=\frac{m}{3}\left(\overline{v^2}(r)\overline{v}^2(r)\right),$$
(3)
and the mean-square velocity is defined in the standard way,

$$\overline{v^2}(r)=\frac{1}{N_k}\underset{i=1}{\overset{N_k}{}}v^2(𝐫_i).$$
(4)
With this information in hand we can calculate the final entropy of the gas. Since at this late stage of expansion the gas is very dilute one can use the ideal-gas formulae. The freeze-out entropy in a given shell is defined as

$$S(r)=N_k\mathrm{ln}\left[\frac{V_ke^{5/2}}{N_k\lambda _T^3}\right],\lambda _T=\left(\frac{2\pi \mathrm{}^2}{mT(r)}\right)^{1/2},$$
(5)
and the total entropy is obviously given by the sum over the shells. We believe that this definition of entropy is valid despite the fact that it is applied not to the real gas but to the ensemble of freeze-out points in the phase space. Here one can use an analogy with the microwave background radiation in the Universe which keeps its entropy constant despite the fact that the photons have decoupled from the matter at the recombination stage a long time ago.

To make statistical errors similar for different systems, the number of generated events is chosen to be inversely proportional to the system’s particle number (see Table 1). This guarantees that the total number of freeze-out points is approximately the same for all considered systems. The dynamical simulations were performed with the time step of 0.5 fm/c which was sufficient to resolve practically all collisions.

The results of the simulations are presented in Figs. 1-4 and Table 1 (see below). Because of the limited statistics the spatial distributions shown in the figures are sensitive to the binning of data. Most calculations were done by sampling freeze-out points in spherical shells of equal volume (r3 binning). This guarantees uniform statistical errors for the bulk parts of distributions but leads to enhanced fluctuations on their tails. Moreover, for unambiguous separation of flow and thermal components the radial bin size should be sufficiently small.

It is necessary to emphasise that the spatial characteristics presented in Figs. 1-4 correspond to sampling of freeze-out points irrespective of the times when particles have actually decoupled from the system. They represent the whole freeze-out history and in this respect differ from the time evolution of the gas characteristics usually presented in gasdynamical calculations. We believe that this representation is more adequate for calculating observable characteristics of small exploding systems. This is especially true for the interpretation of experimental data on energetic nucleus-nucleus collisions.

Fig. 1 shows the spatial density of freeze-out points averaged over all events. In all cases it has a bulk part and a tail. The bulk density is about 0.1 particles/fm3 and almost independent of the system. This should be compared with an initial density of 0.14 particles/fm3. In the tail region the density rapidly drops to zero over a radial distance of about 2-3 fm. Such behaviour should be anticipated from the general consideration of the freeze-out process . The tail is formed by particles emitted through the surface at early times and the inner part contains particles from the bulk freeze-out.

Fig. 2 presents the collective velocity profiles calculated on the freeze-out fields. With a certain degree of imagination one can recognize a Hubble-like behaviour. As expected the collective velocity grows to the outer edge of the distribution. The peak value of about 0.5-0.6 is reached somewhere in the tail region. This value is in good agreement with gasdynamical calculations  predicting for leading particles a velocity of about 3$`c_s`$.

The temperature profiles presented in Fig. 3 are in a certain sense complementary to the flow profiles. One can see that the temperature reaches maximum values at the edge of the bulk region and these values decrease progressively with the system’s size. This can be explained by the fact that the freeze-out process in larger systems develops at later stages of expansion leading to lower freeze-out temperatures. As well known (see e.g. ref. ), in a macroscopic system the freeze-out temperature approaches zero and the whole thermal energy is finally transformed into collective flow. We clearly see the transition from “small” to “large” systems by analyzing the average number of collisions per particle. As we see from our simulations, this number scales roughly as $`\sqrt{N}`$. While it is only about 1 for $`N=50`$ (small system), it is already 2-3 for $`N=400`$ (mesoscopic system), and will be about 10 for $`N=5000`$ (large system).

Finally we come to the most interesting quantity, i.e. the entropy per particle as defined by eq. (5). Fig. 4 shows the corresponding profiles. One can see two clear features. First, the entropy per particle is rather constant over the bulk region and its value varies very little with the system size. Second, there is a significant rise in the entropy per particle in the tail region well above the bulk value 2-3. Moreover, the smaller the system the stronger the rise. This latter trend can be explained by the bigger volume per particle in the outer tail region.

Now we can go back to our discussion of entropy conservation. For this analysis we use the total entropy per particle calculated at freeze-out and compare it with the initial entropy. The latter is calculated by applying the same eq. (5) for the whole gas in the initial volume. The results are presented in Table 1 together with the number of particles, the number of events, the total number of freeze-out points for each size ($`N`$) and the freeze-out fraction $`F_{freeze}=n_{freeze}/(Nn_{event})`$.

From Table 1 one can see that surprisingly enough the initial and final entropies per particle are rather close to each other in all cases. The increase of about 0.5 units is largest for the smallest system considered ($`N=50`$). Formally this is a 20% effect which is quite significant. However, one should bear in mind that in classical statistics the absolute value of entropy is defined up to a constant. This increase in entropy for small systems is an effect of the surface of the system: As seen from Fig. 4, the entropy per particle increases significantly with radius, and since there is relatively more surface in a small system also the total entropy per particle is bigger.

One may wonder, if entropy per particle is conserved, what happens to the total entropy of the system? Indeed, a significant fraction of particles $`(1F_{freeze})`$ leaves the system without scatterings. The simulations show that these ”missing” particles are emitted early in the expansion and come predominantly from the surface region. We have not analyzed their characteristics in detail but we know that these particles were initially in thermal equilibrium with the rest of the system. Thus they should carry away approximnately the same amount entropy per particle as in the initial state. Therefore, we expect that the total entropy of the system is also approximately conserved.

In conclusion, we have used a simple model for a repulsive gas to study the explosive dynamics of small systems. We have demonstrated that in the course of free expansion the temperature drops and collective flow develops in the gas. In contrast to our expectations we have found that the entropy per particle defined on the freeze-out field is almost conserved even in systems with a few hundred particles. This justifies the application of thermal and hydrodynamic models for describing matter evolution in energetic collisions of medium and heavy nuclei. Based on these results we put forward a new interpretation of the old Joule experiment. The gas expansion in this case is approximately isentropic at freeze-out before the particles hit the wall. Then the entropy is produced while the system equilibrates in the larger volume.

The authors thank H. Feldmeier and T, Dossing for productive discussions. This work was supported in part by the Danish National Research Council, Duetsche Forschung Gemeinschaft, DFG grant 436 RUS 113/711/0-1, and Russian Fund of Fundamental Research, RFFI grant 03-02-04007.


text:
# Clarens Client and Server Applications

## I Introduction

The Clarens web services platform, described in a companion article, architecture  acts as a go-between for distributed clients on the wide area network to access services using the widely implemented XML-RPC and SOAP data serialization standards on top of the lower-level HTTP protocol. Its usefulness must ultimately be measured by the usefulness of the services and clients themselves, however.

This paper is divided into two parts describing the currently implemented services as well as clients taking advantage of these services.

## II Services

A rather terse overview of the current Clarens services are given below. The reader is referred to the Clarens web page clarens\_web  for full documentation of all methods and modules.

### II.1 File access

Accessing files on a remote machine remains the most useful service any middleware product can provide. This is evident from the more than 40 million web servers that are deployed worldwidenetcraft . Although the percentage of static files served is hard to gauge, the SPECweb99 specbench  benchmark puts this number at 70% of requests.

Clarens serves files in two different ways: in response to standard HTTP GET requests, as well as via a file.read() service method. A virtual server root directory can be defined for each of the above via the server configuration file which may be any directory on the server system.

The file.read() method takes a filename, an offset and the number of bytes to return to the client. Error message are returned as serialized RPC responses. Network I/O is handed off to the web server, which uses the zero-copy sendfile() system call where available to minimize CPU usage and increase throughput.

Other file access methods include file.ls() to obtain directory listings, file.stat() to obtain file or directory information, and file.md5() to obtain a hash value for checking file integrity.

### II.2 Proxy escrow

Despite the use of asymmetric key cryptography, key management remains a problem, with private keys and certificates (credentials) having to be present on the client system. In the case where the same credentials must be used from different places, e.g. a person’s desktop, laptop and other computer systems, this is inconvenient, as well as degrading the security of the credentials.

In analogy to the MyProxy myproxy  project, Clarens offers the ability to store and retrieve short-lived, self-signed certificates (called proxy certificates in Grid-oriented literature).

The RPC API provides the methods proxy.store(), proxy.retrieve(), and proxy.delete() for managing stored proxy certificates. Combinations of private key/certificate pairs or proxy/certificate pairs may be stored in this way, with the certificate distinguished name (DN) acting as a unique identifier. The credential information is stored in encrypted form using a symmetric cipher using a password provided when invoking the store method. The password itself is not stored, for obvious reasons.

From the above it should be obvious that this approach presents a chicken-and-egg problem of having to present credentials in order to obtain credentials. This may be solved in two slightly different ways. Firstly, a web portal may be constructed that takes input from the user and acts as an intermediary to log into the Clarens server in question and retrieve the credentials. This can be done with the portal residing on an arbitrary machine. Secondly, the ability of Clarens to respond to HTTP GET requests (i.e. act as as simple web server) may be used to construct such a portal on the Clarens server itself, which has access to the server methods as an unprivileged user. These web portals may also be accessed programmatically from within programs or scripts.

Finally, it should be noted that it is up to the user whether to store any credentials in this way, since it may be perceived as less secure than keeping the credentials on possibly multiple systems’ as files or in web browsers.

### II.3 Virtual organization management

The Clarens authentication architecture architecture  is built around the concept of a hierarchical virtual organization (VO) of groups and subgroups of individuals identified by unique distinguished names (DNs). These individuals may be both people or servers.

To ease administration of the VO, a set of methods is provided in the group module to create, delete, and list groups and their members and administrators.

The most important of these methods are described in Table 1.

In addition to managing the VO structure, the group module also provide methods to store, retrieve and search for certificates. Certificate Authority certificates may similarly be searched and retrieved, but not stored, since these certificates are used to ensure the uniqueness of client distinguished names. Instead, CA certificates are managed by the system administrator in the form of files stored in a designated directory.

### II.4 Access control management

The Clarens authorization architecture is built around the concept of hierarchical access control lists for RPC methods architecture . To administer these ACLs, methods are provided to create and modify lists of users and groups allowed access to a particular method, or module.

### II.5 Storage Resource Broker access

The SDSC Storage Resource Broker (SRB) provides a uniform interface to external applications to access various storage media, including local and remote file systems, and tape storage. It is also a client-server system, with Clarens acting as an SRB client on behalf of its own clients. The Clarens SRB API provides methods to initialize a connection to the SRB server and browse, store, and retrieve files.

This interface is not currently deployed, however, since it exposed a critical impedance mismatch between these two client-server systems: Clarens uses an entirely stateless connection protocol, while SRB uses a stateful protocol. This means that in the current implementation a new SRB connection must be initiated upon each method invocation by the Clarens client, which results in very poor performance.

Work is underway to remedy this mismatch by utilizing Clarens agents that can hold persistent connections to SRB on behalf of clients. Another approach being considered is to implement a stateful protocol interface to Clarens itself.

### II.6 Job execution

Next to file access, the ability to execute jobs on a remote machine remains one of the cornerstones of distributed computing. Using the Clarens distinguished name to system user mapping architecture , system commands can be executed by remote users.

As is common practice, Clarens makes use of a small compiled program, called suexec to change the permissions of the resulting process. This program can be more easily audited for security than the entire dependency chain of the main Clarens code. Upon receipt of a shell command, a directory owned by the remote user is created where the command’s output and error messages are stored. The working directory of the resultant process is also changed to this directory, and any newly created files may be accessed remotely using a job ID which is managed by Clarens.

It should be pointed out that this interface is not designed to schedule jobs on Clarens servers, but is most useful for handing jobs to the schedulers (e.g. condor , pbs ) themselves, and retrieving the resultant output files.

This interface is currently used to develop interactive remote analysis using the CMS ORCA analysis package, where analysis jobs themselves become Clarens servers (albeit less featureful ones) that can act as personal remote application servers to allow interactivity with a long-running analysis process on a cluster.

### II.7 Relational database analysis interface

An interface to the Stl based Object Caching And Transport System (SOCATS) being developed by Caltech, allows remote users to query large Physics datasets using standard relational database (SQL) queries.

SOCATS results are returned in the form of a object file formatted as a ROOT root  tree, which can be retrieved by the above-mentioned file access methods. A ROOT interface to Clarens is also described below, which was successfully demonstrated at the 2002 Supercomputing conference.

## III Clients

One of the express aims of the Clarens architecture is to use widely deployed interfaces to lower the implementation barrier for new client applications. The availability of HTTP, SSL and SOAP/XML-RPC implementations on most platforms and programmings languages helps us achieve this goal with the minimum of new code.

### III.1 Python

Python python  is a weakly typed, object oriented scripting language. Its programs are compiled to a platform-independent byte-code, similar in many ways to Java. It’s built-in support for both HTTP, SSL, and XML-RPC, combined with the rapid prototyping abilities inherent in a scripting language made it natural to be used as the default client-side development language.

The Python Clarens client is implemented as a pure Python class, called clarens\_client that takes an argument for the server URL, and optionally the certificate and private key files to be used in the connection to its constructor method. The constructor method initiates a connection with the server an authenticated the user using the credentials stored in the standard places by the Globus toolkit, including proxy credentials if they exist, or those provided to the constructor method.

The clarens\_client object maps all non-local method calls to remote procedure calls, handling serializing/deserializing of the method arguments and return values transparently. The following example of the echo.echo() method demonstrates this:

```
>>> obj=clarens_client("http://server.org/")
>>> obj.echo.echo("Hello")
["Hello"]

```

I.e. the echo.echo() remote method is invoked directly from the command line, and the result is returned as a native Python string to the caller. The result is always returned as a list, indicated by the square brackets.

### III.2 ROOT

The object-oriented modular architecture of the ROOT root  analysis package, combined with its rich set of built-in objects and wide adoption in the high energy physics community makes it a very useful client for remote analysis functionality.

In analogy to the Python client ,Clarens ROOT client handles authentication transparently, and provides the user with a Clarens object to communicate wit the remote server. This client does not do automatic serialization/deserialization of arguments and return values to native ROOT objects, a lower level interface must be used instead for general remote procedure calls.

It does, however, provide a convenient interface to read remote files using the Clarens server, via a TCWebFile object derived from the native TFile object. This object provides all the functionality of the local version, allowing transparent analysis of remote files from the ROOT command line, scripts or compiled code by changing only the object type. Using the interactive ROOT object browser it is possible to browse the structure and content of remote ROOT files quickly and easily, with the ability to display histograms and other types of plots contained in the remote files interactively or programmatically.

The TFile base class supports a dynamic local caching mechanism that is used to optimize file transfers, so that the extremes of transferring the whole remote file to the local client, as well as making a remote procedure call for small parts of the file can be avoided, striking a balance between bandwidth and latency constraints.

Another convenience class is the TCSystemDirectory, derived from the TSystemDirectory base class. This class provides an interface to the directory browsing functionality in the Clarens file module, allowing the client to interactively traverse the remote directory structure using the ROOT object browser. Any remote ROOT files may be opened by clicking on their icons in the browser. The TCSystemDirectory class may also be used programmatically from ROOT scripts or compiled C++ programs.

### III.3 IGUANA

IGUANA is an interactive visualization toolkit, used for detector and event visualization, as well for interactive data analysis and presentation. It contains support for accessing remote OpenInventor-formatted files via Clarens for local display and manipulation.

### III.4 Browser interface

The web browser is currently the most widely used distributed computing tool, bar none. Modern browsers have native support for SSL-encrypted connections and client-side certificate authentication, making it an ideal platform for a Clarens interface.

After initial experiments with client-side Java applets for communicating with the server, it was decided to use the Javascript language embedded in most browsers to handle this task. Since Clarens is able to serve web pages in response to HTTP GET requests, the browser interface is implemented as a series of static web pages that embed Javascript scripts to handle communication and interface display using dynamic HTML. This implementation eliminates the need for clients to install any additional software apart from a web browser, which most people already have.

The browser interface uses XML-RPC for data serialization since it is by far simpler than the more complex SOAP protocol. As with the Python interface, argument and return value serialization/deserialization is handled transparently by the provided Javascript libraries, made easy by the object oriented, loosely typed nature of the language.

Functionality currently provided include browsing a remote file repository, with the ability to download files, and virtual organization management.

## IV Conclusion

Clarens provides a growing list of services and useful client implementations for doing distributed computing in a Grid-based environment.

###### Acknowledgements.

This work supported by Department of Energy contract DE-FC02-01ER25459, as part of the Particle Physics DataGrid project ppdg , and under National Science Foundation Grant No. 0218937. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors, and do not necessarily reflect the views of the National Science Foundation. Clarens development is hosted by SourceForge.net sourceforge .


text:
# Boundary Slip as a Result of a Prewetting Transition

## I Introduction

It is accepted in hydrodynamics that the velocity of a liquid immediately adjacent to a solid is equal to that of the solid batchelor.gk:1967.a . Such an absence of a jump in the velocity of a simple liquid at a surface seems to be a confirmed fact in macroscopic experiments. However, it is difficult to obtain the same conclusion using microscopic models. It has been noticed that, even in case of simple liquids, the no-slip boundary condition is not justified on a microscopic level.

Therefore, the conventional condition of continuity of the velocity, or the no-slip boundary condition, is not an exact law but a statement of what may be expected to happen in normal circumstances. While the normal component of the liquid velocity must vanish at an impermeable wall for kinematic reasons, the requirement of no-slip can be relaxed. In other words, instead of imposing a zero tangential component of the liquid velocity at the solid, it is possible to allow for an amount of slippage, described by a slip length $`b`$. The slip length for a simple shear flow is the distance behind the interface at which the liquid velocity extrapolates to zero

$$v_s=b\left[_zv(z)\right]_{\mathrm{wall}},$$
(1)
where $`v_s`$ is the tangential velocity at the wall, and the $`z`$ axis is perpendicular to the surface. The definition of $`b`$ is explained in Fig. 1. It is clear that boundary slip is important only when the length-scale over which the liquid velocity changes approaches the slip length. Therefore, it is not surprising that the slippage effect has not been detected in macroscopic experiments. In microfluidic devices, however, where the liquid is highly confined, the boundary slip is important vinogradova.oi:1999.a .

Indeed, water flow in capillaries of small diameter with smooth hydrophobic walls has been investigated and slip at the wall had to be postulated to recover the viscosity of water churaev.nv:1984.a ; watanabe.k:1999.a . These results were confirmed by directly probing the fluid velocity at a solid surface using total internal reflection-fluorescence spectroscopy pit.r:2000.a ; tretheway.dc:2002.a  as well as double focus confocal fluorescence cross-correlation lumma.d:2003.a . Several indirect methods were also used: the quartz-crystal microbalance (QCM) krim.j:1996.a , the surface force apparatus (SFA) horn.rg:2000.a ; zhu.yx:2001.a ; zhu.yx:2002.a ; baudry.j:2001.a , and the atomic force microscope (AFM) craig.vsj:2001.a ; vinogradova.oi:2003.a . The magnitude of the slip length $`b`$ was sometimes greater than $`100\mathrm{n}\mathrm{m}`$ for partially wetted walls pit.r:2000.a ; zhu.yx:2001.a . In some cases the shear rate did not affect the amount of slip in the observed range pit.r:2000.a ; vinogradova.oi:2003.a ; in others a strong dependence on the velocity was found horn.rg:2000.a ; zhu.yx:2001.a . It was also shown that both surface roughness and strength of the fluid-surface interactions affect the wall slip churaev.nv:1984.a ; zhu.yx:2002.a ; vinogradova.oi:2003.a .

The no-slip condition can also be violated in more complex systems. Boundary slip has been suggested for polymer melts brochard.f:1992.a ; ajdari.a:1994.a ; brochardwyart.f:1996.a  and liquid crystals francescangeli.o:1999.a  (for the rotational motion of molecules).

The origin of such large slippages remains unclear despite considerable theoretical effort. On the theoretical side, molecular dynamics simulations have shown that the molecules can slip directly over the solid due to the fact that the strength of attraction between the liquid molecules is greater than the competing solid-liquid interaction sun.m:1992.a ; sun.m:1992.b ; bocquet.l:1993.a ; thompson.pa:1997.a ; barrat.jl:1999.a ; barrat.jl:1999.b . In general, wall slip was found on non-wetted surfaces, i.e. when the contact angle is large. However, the simulation results were not entirely consistent with the experimental data, by predicting a much lower slip length cieplak.m:2001.a ; sokhan.vp:2001.a  and substantial slippage only at large contact angles barrat.jl:1999.b . It has also been demonstrated that the surface roughness may both reduce cottin\_bizonne.c:2003.a  and increase hocking.lm:1976.a ; ponomarev.iv:2003.a  the friction of the fluid past the boundaries.

Other ideas invoke the formation of a new phase at the wall. The possible source of the surface layer could be a gas (lubricant) dissolved in the liquid, forming bubbles nucleating at the liquid/solid interface ruckenstein.e:1983.a ; degennes.pg:2002.a ; vinogradova.oi:1995.a . Experimental evidence for the formation of a wetting layer under flow has been found by Tanaka tanaka.h:2001.a ; tanaka.h.1993.c . The boundary layer can have a lower viscosity than the bulk value vinogradova.oi:1995.b . Tuning the size and the properties of this layer one can obtain large values of the slip length (see Fig. 2).

Indeed, in a sharp interface limit, when the width of the interface is much smaller than the width of the slab and the thickness of the boundary layer, we can neglect the structure of the interface. Then the problem is reduced to the shear flow of two phases ($`s`$ for “surface”, $`b`$ for “bulk”) with viscosities $`\eta _s`$ and $`\eta _b`$, respectively, and thicknesses $`\delta `$ and $`l`$, respectively (see Fig. 2). Denoting the velocity profiles in the surface layer and the bulk by $`v_s(z)`$ and $`v_b(z)`$, respectively, and using stick boundary conditions at $`z=0`$ and $`z=l+\delta `$, we have

$`v_s(0)`$ $`=`$ $`0,`$
$`v_s(\delta )`$ $`=`$ $`v_b(\delta ),`$ (2)
$`v_b(l+\delta )`$ $`=`$ $`v_0,`$
where the second equation is the condition for continuity of the flow field at the interface between the phases, and $`v_0`$ is the velocity at the top of the slab. Furthermore, we use the condition of zero divergence of the shear stress tensor at the interface,

$$_z\left(\eta (z)_zv(z)\right)=0,$$
(3)
or

$$\eta _sv_s^{}(\delta )=\eta _bv_b^{}(\delta ).$$
(4)
The solution in the bulk

$$v_b(z)=\kappa v_0\left(\frac{\eta _b}{\eta _s}1+\frac{z}{\delta }\right),$$
(5)
with $`\kappa =\left(\eta _b/\eta _s+l/\delta \right)^1`$, results in a slip length vinogradova.oi:1995.b

$$b=\delta \left(\frac{\eta _b}{\eta _s}1\right).$$
(6)
According to Eq. (6) the boundary slip can be observed if the viscosity depends on the composition ($`\eta _b\eta _s`$) and the less viscous fraction of the liquid wets the walls better than the more viscous one ($`\eta _b>\eta _s`$). It is also clear that there are two ways to obtain a large slip length. First, by having a macroscopically thick boundary layer, since the slip length has the same order of magnitude as the thickness of this layer. Second, by providing large values of the viscosity contrast $`\eta _b/\eta _s`$, e.g. for a gas layer 111when the gas is in the Knudsen regime, the slip length does not depend on the thickness of the boundary layer degennes.pg:2002.a ..

A more realistic description should allow for a prewetting transition binder.k:1983.a ; bonn.d:2001.a  for the liquid/gas or liquid/lubricant mixture and take into account the structure and the finite width of the interface region.

The aim of this paper is to include these effects and relate the slip length to the wettability of the walls, composition of the mixture, and thermodynamic parameters of the system. We show that the prewetting transition can give rise to a large boundary slip of the fluid by generating a thick film of a phase-separated ‘lubricant’ at the wall which has a lower than the bulk fluid viscosity. Indeed, if we choose typical values, thickness of the wetting layer (thick film) $`\delta 10\mathrm{n}\mathrm{m}`$, viscosity contrast $`\eta _b:\eta _s=3:1`$, we find $`b20\mathrm{n}\mathrm{m}`$, i.e. the prewetting film can indeed give a large slip length. This value can be further increased. Indeed, if the phase separation occurs upon cooling (heating), the thickness of the wetting layer increases with the increase (decrease) of the temperature: for a system with short-range forces in the vicinity of the wetting temperature $`T_w`$, the thickness diverges as bonn.d:2001.a  $`\delta \mathrm{ln}\left(\left|T_wT\right|\right)`$.

The paper is organized as follows. In Sec. II we outline the approach that allows us to calculate the phase diagram, order parameter profiles, and slip length. Section III gives a summary of the results. Finally, in Sec. IV, we present some brief conclusions.

## II Theoretical Model

Phase separation phenomena in binary and polymer mixtures have been intensively studied by theory, experiment, and simulation. While the most detailed information is available for the static bulk behavior, much more interesting phenomena occur when studying the dynamics gunton.jd:1983.a ; bray.aj:1994.a ; wagner.aj:1998.a ; kendon.vm:1999.a , or the behavior near surfaces and in confined geometries binder.k:1983.a ; dietrich.s:1988.a ; binder.k:2003.a ; duenweg.b:2003.a . The theoretical understanding of phenomena which combine both aspects (i. e. dynamics near surfaces) is still at its infancy puri.s:1994.a , while there are many experiments tanaka.h:1993.a ; tanaka.h.1993.c ; vinogradova.oi:1995.b .

Our theoretical approach splits the problem of shear flow of a binary mixture near surfaces into several independent tasks. First, we calculate the equilibrium order parameter profiles, completely disregarding the flow. Restriction to equilibrium thermodynamics allows us to introduce a change of ensemble: we fix the chemical potential difference (semi-grand-canonical ensemble) instead of fixing the composition, which is conceptually and computationally easier. The order parameter profile then results in a viscosity profile, which in turn allows us to calculate the stationary velocity profile by solving Eq. (3), again using stick boundary conditions at both surfaces. It is not immediately obvious that this split–up is justified. We discuss the restrictions and underlying assumptions of this approach in Appendix A. Finally, the slip length is calculated from the stationary velocity profile.

### II.1 Free Energy of a Binary Mixture

To describe the bulk phase as well as the interface structure, ‘phase-field’ models binder.k:1983.a  are often used. In this approach the order parameter $`\varphi `$ is introduced. For a binary mixture $`\varphi `$ is a composition variable, defined as $`\varphi =(n_1n_2)/(n_1+n_2)`$, where the $`n_i`$ are the number densities of the two species. This order parameter varies slowly in the bulk regions and rapidly on length scales of the interfacial width. The unmixing thermodynamics is described via a free energy functional.

Since the material is confined in a container in any experiment, phase separation is always affected by surface effects evans.r:1986.a ; evans.r:1987.a ; gelb.ld:1997.a ; gelb.ld:1997.b ; gelb.ld:1999.a . To include them, appropriate surface terms responsible for the interaction of the liquid with the container walls are added to the free energy cahn.jw:1977.a ; binder.k:1983.a ; degennes.pg:1985.a ; dietrich.s:1988.a ; bonn.d:2001.a ; schick.m:1990.a .

In the phase-field approach the semi-grand potential of a binary mixture is written as bonn.d:2001.a

$$\mathrm{\Omega }\{\varphi \}=\frac{1}{a^3}𝑑V\left[\frac{k}{2}a^2\left(\varphi \right)^2+f(\varphi )\mu \varphi \right]+\mathrm{\Psi }_s,$$
(7)
where $`a`$ is a normalization length of the order of the size of a molecule, $`f(\varphi )`$ is the Helmholtz free energy density of the mixture, while $`\mu `$ is the chemical potential thermodynamically conjugate to the order parameter $`\varphi `$, and $`\mathrm{\Psi }_s`$ is the surface energy.

The explicit form of the Helmholtz free energy $`f(\varphi )`$ varies depending on the type of mixture. However, the simple observation that the two phases must coexist implies that there are two minima in the free energy at the respective values of the order parameter. We here adopt the mean-field model for a regular (symmetric) mixture reichl.le:1998.a ; rowlinson.js:1969.a

$`f(\varphi )`$ $`=`$ $`{\displaystyle \frac{\chi }{4}}(1\varphi ^2)+`$ (8)
$`k_BT\left[{\displaystyle \frac{1+\varphi }{2}}\mathrm{ln}{\displaystyle \frac{1+\varphi }{2}}+{\displaystyle \frac{1\varphi }{2}}\mathrm{ln}{\displaystyle \frac{1\varphi }{2}}\right],`$
where the first term on the right-hand side corresponds to the excess energy of mixing. Note that this is one of the simplest models to describe unmixing; a more realistic description would need a more sophisticated function, which also takes into account a dependence on the overall density, which can vary along the profile (see Appendix A).

The term $`\left(\varphi \right)^2`$ is needed to provide spatial structure to the theory: at phase coexistence, there are two bulk equilibrium order parameter values $`\varphi _+`$ and $`\varphi _{}`$ with the same free energy density, $`f(\varphi _+)=f(\varphi _{})`$. Without the gradient term, a structure with very many interfaces between the $`\varphi _+`$ and $`\varphi _{}`$ phase would be entropically favored. The term $`(k/2)\left(\varphi \right)^2`$ is the simplest one which penalizes interfaces. While this is justified near the critical point, where interfaces are very wide and the order parameter varies smoothly, a more realistic description at strong segregation (where the interface becomes rather sharp) would require higher–order gradients, too.

### II.2 Surface Free Energy

To describe the interaction with the walls we used the quadratic approximation for the surface energy bonn.d:2001.a ; flebbe.t:1996.a

$$\mathrm{\Psi }_s=\frac{1}{a^2}\left[h\varphi _s\frac{1}{2}\gamma \varphi _s^2\right]𝑑S,$$
(9)
where $`\varphi _s`$ is the surface value of the order parameter and the parameters $`h`$ and $`\gamma `$ are referred to as the short-range surface field and the surface enhancement, respectively. The short-range surface field, $`h`$, is a measure of the attractiveness (or repulsiveness, if negative) of the surface to the component $`1`$. In real systems it can be of either sign and of any magnitude. The surface coupling enhancement, $`\gamma `$, represents the effect that a molecule close to the substrate has fewer neighbors than a molecule in the bulk; $`\gamma `$ is estimated to be small and negative. We have implicitly assumed that all surface effects are of short range, thus $`f_s`$ depends on the local concentration at the walls only. Because of long range van der Waals forces this is not fully realistic. However, for large enough wall separations the differences to the short range case are rather minor.

We introduce dimensionless units by setting $`a=1`$, $`\chi =1`$, and $`k_B=1`$. Hence, energies are given in units of $`\chi `$, temperatures in units of $`\chi /k_B`$, and lengths in units of $`a`$. Furthermore, a value of $`k=1`$ has been used throughout the calculations of this paper.

### II.3 Euler-Lagrange Equations

In thermal equilibrium the grand potential (7) must be minimal. Variation of this functional yields an Euler-Lagrange equation together with two boundary conditions. Due to translational symmetry in $`x`$ and $`y`$ direction, the problem is one-dimensional. In contrast to the situation considered in Introduction, we now focus on the case of two identical walls separated by a distance $`L`$, which is chosen large enough such that the two wetting layers do not overlap, and bulk behavior is established in the center of the slab. It is then convenient to choose the origin of the coordinate system at the center of the film. In this coordinate system the Euler-Lagrange equation and boundary conditions read

$`k{\displaystyle \frac{^2\varphi }{z^2}}+{\displaystyle \frac{1}{2}}\varphi {\displaystyle \frac{1}{2}}T\mathrm{ln}{\displaystyle \frac{1+\varphi }{1\varphi }}+\mu =0,`$ (10)
$`\pm k{\displaystyle \frac{\varphi }{z}}+h+\gamma \varphi |_{z=L/2,L/2}=0.`$ (11)
This boundary-value problem was solved using the relaxation method press.wh:1992.a . In general, it can have up to three different solutions: one stable, one metastable, and one unstable. The relaxation method yields only the metastable and the stable solution. The unstable solution with the highest free energy has negative response function $`(\varphi /\mu )_T`$ and is eliminated by the relaxation method automatically.

To select stable solutions, we calculated the grand potential of the mixture, $`\mathrm{\Omega }`$, for both stable and metastable solutions and chose the solution with the lowest grand potential. This allows the accurate determination of the phase diagram.

### II.4 Velocity Profiles and Slip Length

To calculate the velocity profiles, we assumed that the viscosity of the system is simply a linear combination of that of the individual components

$$\eta (z)=\eta _s\frac{1+\varphi (z)}{2}+\eta _b\frac{1\varphi (z)}{2},$$
(12)
where the viscosity contrast between the two components has been chosen as $`\eta _s:\eta _b=1:3`$. Of course, a more realistic description would have to introduce a more complicated dependence on $`\varphi `$, and also take into account the dependence on the overall density. The stationary velocity profile is the solution of Eq. (4)

$$v(z)=v_\mathrm{w}c^1_{L/2}^z\frac{dz}{\eta (z)},$$
(13)
where

$$c=_{L/2}^{L/2}\frac{dz}{\eta (z)},$$
(14)
and we assumed stick boundary conditions at the walls, $`v(L/2)=0`$ and $`v(+L/2)=v_\mathrm{w}`$. The value of the slip length was obtained by fitting the bulk region of the velocity profile (13) with a linear regression.

## III Results and Discussion

### III.1 Order Parameter Profiles and Prewetting Phase Diagram

Typical order parameter profiles at the same value of the chemical potential $`\mu =0.002`$ and a set of temperatures are shown in Fig. 3. Only the profiles corresponding to the stable thermodynamic state are shown. The profiles have a well-developed flat region in the center of the film, which indicates that there should not be finite size effects for the chosen thickness of the slab ($`L=200`$). This flat region is important for defining the slip length, see Eq. (1). This requires a well-defined linear velocity profile in the bulk.

In order to understand these profiles, one should note that the very small value of $`\mu `$, combined with the low temperature, implies that the system is rather close to bulk coexistence. Bulk coexistence, however, is characterized by two equilibrium values of the order parameter with large absolute value and opposite sign. The small negative value of $`\mu `$ singles out the negative order parameter value, while the absolute value is only slightly changed. Introducing the surface with preference of the other phase, one obtains a profile which is slightly bent up near the surface. The prewetting transition (see Refs. cahn.jw:1977.a ; schick.m:1990.a ) occurs upon increasing the temperature, and results in a sudden increase of the surface excess coverage, i. e. it is a first-order transition between a thin-film and a thick-film state. This jump is directly observed at a temperature near $`T=0.35`$ in the profiles of Fig. 3. Taking into account the remarks made in the introduction and Eq. (6), we can already anticipate small slip lengths for thin films and large slip lengths for thick films (above the prewetting transition temperature). Since the thin-thick film transition is a first-order transition, one can also expect a jump in the slip length at the prewetting transition temperature.

To calculate the thin-thick film transition temperature, we followed the metastable solution up to its stability limit (spinodal) calculating the grand potential for both stable and metastable solutions. The transition line was determined from the intersection of grand potentials of stable and metastable solutions (thick and thin films). Both spinodals as well as the transition line are shown in a prewetting transition phase diagram, Fig. 4. The solid curve is the prewetting curve ending at the prewetting critical point and the dashed curves are spinodals or metastability limits of the metastable states.

### III.2 Slip Length

The temperature dependence of the slip length for several values of the chemical potential (or average volume fraction) is presented in Fig. 5. As we anticipated, $`b(T)`$ has a jump at the prewetting transition temperature, when the thick film is formed (path $`\mathrm{I}`$ in Fig. 4). Before the jump the slip length is small and practically does not depend on temperature. After the transition the slip length decreases with the temperature increasing, even though the prewetting film gets thicker. This is because the bulk volume fraction (order parameter $`\varphi `$) increases, giving rise to a decrease in the bulk viscosity (or, equivalently, viscosity contrast) (see Fig. 3).

If we are above the prewetting critical point (path $`\mathrm{II}`$ in Fig. 4) there is no jump-like transition, but rather a smooth increase of the slip length with the temperature.

The dependence of the slip length on the chemical potential is shown in Fig. 6. If the temperature is below $`T_\mathrm{c}^\mathrm{w}`$ then the thin film is always stable and the slip length is small (path A in Fig. 4). Once we intersect the prewetting transition line, the thick film forms with the increase of the average fraction of the more wettable phase (values of the chemical potential close to zero) and the slip length jumps abruptly to higher values (path B in Fig. 4). Above the prewetting critical point $`T_\mathrm{c}^{\mathrm{pw}}`$ the slip length increases monotonically with the increase of the chemical potential (see path $`\mathrm{C}`$ in Fig. 4).

Finally, the dependence of the slip length on the surface field $`h`$ is shown in Fig. 7. It is clear that there is a threshold value of the surface field (contact angle) when the thick film is formed. At this value of the surface field the slip length increases abruptly. Below the threshold the slip length is small and practically does not depend on the surface field.

## IV Conclusions

To summarize, when the prewetting transition occurs in a flow experiment, it may indeed generate a strong slippage. Prewetting provides a mechanism of generating a macroscopically thick film at the wall. If this film has a lower viscosity than the bulk value, a strong slippage can be observed above the prewetting transition temperature. The value of the slip length has a jump-like dependence on temperature, concentration of the phase-separated liquid, and surface field (contact angle).

The mean-field model of wetting is rather general and can be applied to liquid-gas systems, binary mixtures, as well as incompressible polymer mixtures in the long wavelength approximation flebbe.t:1996.a . This implies that the large boundary slip due to prewetting can be observed in all these systems.

Another goal of this paper is to stimulate accurate quantitative measurements of the slip length combined with measurements of the fluid wetting properties. This will allow to study the underlying microscopic mechanisms of slippage, and choose an adequate model for every experimental situation. In particular, it would be very interesting to measure the slip length as a function of temperature while the system undergoes a prewetting transition, with the adsorbed species having a lower viscosity. Our results indicate that in such a case the slip length should undergo an abrupt change, and we believe that this will be true beyond the limitations of our idealized model.

## Appendix A Hydrodynamics of Slippage in Binary Fluids

We consider a binary fluid (species $`1`$ and $`2`$) confined between two infinite parallel planes to a thin slab of thickness $`L`$. The plane normal is in the $`z`$ direction. We denote the mass densities of the two species by $`\rho _1`$ and $`\rho _2`$, respectively, and introduce the linear combinations $`\rho =\rho _1+\rho _2`$ (total mass density) and $`\mathrm{\Delta }\rho =\rho _1\rho _2`$.

Concerning the thermodynamics of the system, we notice that the pressure tensor $`p_{\alpha \beta }`$ (Greek letters denote Cartesian indices) is anisotropic as a result of the finite size effect and the interface contribution to the free energy allen.mp:2000.a . Since the system is fluid, there is no elastic response to shear, and hence $`p_{\alpha \beta }=0`$ for $`\alpha \beta `$. Furthermore, $`p_{xx}=p_{yy}`$ for symmetry reasons. Apart from this anisotropy, we assume that there is no further inherent anisotropy. In particular, we assume that all transport coefficients (interdiffusion coefficient, thermodiffusion coefficient, viscosities, etc.) retain their simple scalar nature as in the isotropic macroscopic bulk fluid. The viscous stress tensor is hence written as

$$\sigma _{\alpha \beta }=\eta \left(_\alpha v_\beta +_\beta v_\alpha \frac{2}{3}\delta _{\alpha \beta }_\gamma v_\gamma \right)+\zeta \delta _{\alpha \beta }_\gamma v_\gamma ,$$
(15)
where the Einstein summation convention is implied, $`𝒗`$ denotes the fluid flow velocity, while $`\eta `$ and $`\zeta `$ are the shear and the bulk viscosity, respectively. We assume that these parameters depend on density and composition, i. e. $`\eta =\eta (\rho ,\mathrm{\Delta }\rho )`$ and $`\zeta =\zeta (\rho ,\mathrm{\Delta }\rho )`$.

We consider the situation that the fluid is weakly driven such that a shear flow in $`x`$ direction develops, with shear gradient in $`z`$ direction. In this limit of weak driving, it is reasonable to assume that no symmetry is broken except translational invariance in $`z`$ direction. The system remains translationally invariant in $`x`$ and $`y`$ direction, and gradients (of any quantity) occur only in $`z`$ direction. We also assume that the system is kept at constant temperature throughout, i. e. that the heat production is negligibly small. This is reasonable for small shear rates (note that the heat production is proportional to the square of the shear rate). In addition, the heat conductivity is quite large for many real fluids.

Under these assumptions, we seek a stationary solution of the hydrodynamic equations of motion landau.ld:1995.a  for the outlined geometry and symmetry. Since the velocity flow field $`𝒗=(v_x,0,0)`$ is defined as the center–of–mass velocity of a volume element, the dynamics of $`\rho `$ is governed by pure convection:

$$_t\rho +_\alpha \left(\rho v_\alpha \right)=0.$$
(16)
A stationary solution implies $`_t\rho =0`$, while our geometry leads to $`_\alpha \left(\rho v_\alpha \right)=_z\left(\rho (z)v_z\right)=0`$. The continuity equation is therefore identically fulfilled, and a non–constant mass density profile $`\rho (z)`$ is permitted.

For the density difference there is also interdiffusion:

$$_t\mathrm{\Delta }\rho +_\alpha \left(\mathrm{\Delta }\rho v_\alpha \right)=_\alpha j_\alpha ,$$
(17)
where $`𝒋`$ is the interdiffusion current. Again, the left hand side vanishes identically for our flow. Furthermore, the interdiffusion current vanishes in the stationary state:

$$𝒋=0.$$
(18)
The (full nonlinear) Navier–Stokes equation is written as

$$_t\left(\rho v_\alpha \right)+_\beta \left(\rho v_\alpha v_\beta \right)=_\beta p_{\alpha \beta }+_\beta \sigma _{\alpha \beta }.$$
(19)
For our flow, the left hand side vanishes identically, and hence

$`_zp_{zz}`$ $`=`$ $`_z\sigma _{zz}=0,`$ (20)
$`0=_zp_{xz}(z)`$ $`=`$ $`_z\sigma _{xz}=_z\left(\eta _zv_x\right),`$ (21)
$`0=_zp_{yz}`$ $`=`$ $`_z\sigma _{yz}=0.`$ (22)
From Eq. (20) one sees that the pressure profile $`p_{zz}`$ must be constant, while Eq. (21) allows to obtain the velocity profile via integration, as soon as the viscosity profile $`\eta (z)`$ is known from the profiles $`\rho (z)`$ and $`\mathrm{\Delta }\rho (z)`$.

We now turn to the constitutive equation for the interdiffusion current $`𝒋`$. In non–equilibrium thermodynamics, the dissipative currents are assumed to be linear in the gradients of the intensive thermodynamic variables. A binary system has three independent thermodynamic variables, for which we can take any appropriate set. For our purposes it is particularly useful to choose the pressure $`p`$, the temperature $`T`$, and the chemical potential $`\mu `$, which is defined as the variable which is thermodynamically conjugate to the order parameter $`\varphi =(n_1n_2)/(n_1+n_2)`$, where the $`n_i`$ denote the particle number densities (see main text). Macroscopically, the constitutive equation would thus read

$$j_\alpha =\mathrm{\Gamma }_1_\alpha \mu \mathrm{\Gamma }_2_\alpha T\mathrm{\Gamma }_3_\alpha p,$$
(23)
where $`\mathrm{\Gamma }_i`$ are suitable scalar Onsager coefficients. Note that the pressure gradient may appear since $`\mu `$ is not the variable conjugate to $`\mathrm{\Delta }\rho /\rho `$ (which is often used in the literature landau.ld:1995.a ), but rather to $`\varphi `$. We now generalize this equation to the case of an anisotropic pressure tensor, but retain the scalar nature of the Onsager coefficients in accordance with our assumptions (similar to what is done in model H kendon.vm:2001.a ; zhang.zh:2001.a ). We thus find

$$j_\alpha =\mathrm{\Gamma }_1_\alpha \mu \mathrm{\Gamma }_2_\alpha T\mathrm{\Gamma }_3_\beta p_{\alpha \beta },$$
(24)
or, taking into account that there are only gradients in $`z`$ direction, and that $`𝒋`$ vanishes,

$$0=\mathrm{\Gamma }_1_z\mu \mathrm{\Gamma }_2_zT\mathrm{\Gamma }_3_zp_{zz}.$$
(25)
Now, $`_zT`$ vanishes due to our assumption of an isothermal system, while $`_zp_{zz}`$ vanishes as a consequence of the Navier–Stokes equation (see Eq. 20). For this reason, the profile of the chemical potential, $`\mu (z)`$, must be a constant, too.

In summary, we find that under the given assumptions (isothermal system, isotropic Onsager coefficients, weak driving) the conditions for the stationary state are identical to those in thermal equilibrium (all intensive variables must have constant profiles). This permits to first calculate the density profiles just as equilibrium profiles, completely disregarding the flow, and then, in a second step, to calculate the velocity profile by solving Eq. (21). This has been done in the main text for a simple model for the unmixing thermodynamics, and a linear dependence of the viscosity $`\eta `$ on the order parameter. Note also that the restriction to pure equilibrium thermodynamics allows us to introduce a change of ensemble: instead of considering the composition as fixed, we rather view $`\mu `$ as fixed (semi–grand ensemble), which is conceptually and computationally easier. Of course, a quantitative comparison with experiments is not possible for such simple models; one would have to use much more sophisticated free energy functionals, and a better model for the concentration dependence of $`\eta `$ as well. Furthermore, one should expect that at moderate shear rates only hydrodynamic instabilities (e. g. bubble formation near the surfaces) should occur which invalidate the assumption of translational invariance in $`x`$ and $`y`$ direction.

###### Acknowledgements.

We are grateful to R. Evans, F. Feuillebois, K. Kremer, M. Müller, and J. Vollmer for useful discussions. DA acknowledges the support of the Alexander von Humboldt Foundation.