Near-wall dynamics of concentrated hard-sphere suspensions: comparison of evanescent wave DLS experiments, virial approximation and simulations

Abstract

In this article we report on a study of the near-wall dynamics of suspended colloidal hard spheres over a broad range of volume fractions. We present a thorough comparison of experimental data with predictions based on a virial approximation and simulation results. We find that the virial approach describes the experimental data reasonably well up to a volume fraction of ϕ ≈ 0.25 which provides us with a fast and non-costly tool for the analysis and prediction of evanescent wave DLS data. Based on this we propose a new method to assess the near-wall self-diffusion at elevated density. Here, we qualitatively confirm earlier results [Michailidou, et al., Phys. Rev. Lett., 2009, 102, 068302], which indicate that many-particle hydrodynamic interactions are diminished by the presence of the wall at increasing volume fractions as compared to bulk dynamics. Beyond this finding we show that this diminishment is different for the particle motion normal and parallel to the wall.

---

1 Introduction

Soft matter at interfaces is an essential component of many biological, chemical, and industrial processes. The effect of interactions with system boundaries is even more pronounced as modern technology zooms into smaller length scale, where confinement geometry is comparable to particle size. Some interesting examples are antifreeze proteins at water–ice interface,¹ ‘contact killing’ of bacteria by copper surfaces,² design of food with novel texture by structuring water–water interfaces,³etc. The particular case of colloidal particles near a flat solid wall is prominent in the reduced-scale world of micro- and nanofluidics, for example in lab-on-chip applications, where colloidal particles may be used to manipulate fluid flow. The phase behaviour and structures formed in colloidal suspensions near walls have been investigated thoroughly using X-ray and neutron scattering techniques with grazing incidence,⁴ while static interactions of colloidal particles with solid surfaces were investigated using e.g. total internal reflection microscopy.⁵ On the other hand, experimental investigations of near wall colloidal dynamics have been lagging behind theoretical developments for a long time. The first theoretical considerations of the problem of slow viscous motion of a sphere close to a wall date back to the early twentieth century by Lorentz⁶ and Faxén,⁷ while first experiments were achieved only in the 1980's.⁸ Only during the last twenty years, dynamics at interfaces has developed into a major research branch.

The motion of colloidal particles is known to be hindered⁹ in the vicinity of a wall due to hydrodynamic interactions (HI). Their effect may be probed by a number of experimental techniques, with the method of choice depending on the system, its size and optical properties. For a review of these methods, see ref. 10. Evanescent wave dynamic light scattering (EWDLS) is an important tool to study near-interface dynamics of colloids, and it is the only technique which is available for the investigation of colloids with a size in the 100 nm range. In a typical EWDLS experiment, a laser beam is totally reflected off a glass–solution interface, and an evanescent wave is then created as illumination source. The penetration depth of the evanescent wave can be tuned by varying the incident angle. Particles located in the volume illuminated by the evanescent wave scatter light which is collected by a detecting unit and passed down to a correlator to generate the intensity time autocorrelation function (IACF). Since the method has been devised by Lan et al.,¹¹ it has witnessed rapid development. In early attempts, EWDLS has been applied to study translational diffusion of spherical colloids in dilute solutions.^12–17 With a set-up which allows independent variation of the components of the scattering vector Q_∥, Q_⊥, parallel and perpendicular to the surface, respectively, it is possible to distinguish between the anisotropic diffusivity of colloidal spheres in these directions experimentally.^18,19 EWDLS has also been employed to study the dynamics of stiff polymers adsorbed to the interface,²⁰ as well as the collective motion of end-grafted polymer brushes;^21,22 near-wall diffusion of a spherical particle in a suspension of rod-like depletants²³ and colloidal dumbbells;²⁴ and rotational diffusion of optically anisotropic spheres.^25,26 Notably, evanescent waves have also been used for near-wall nano-velocimetry,^27,28 and to probe dynamics at liquid–liquid interfaces.²⁹

Recent years have brought increasing interest into the effects of confinement on collective dynamics of colloids. To this end, EWDLS experiments have been performed on hard-spheres suspensions with volume fractions up to 42 percent by Michailidou et al.,^30,31 along with theoretical developments.³² On the basis of a heuristic approximation for the near wall self-diffusivity, these works suggest that for a concentrated suspension, many-particle hydrodynamic interactions are diminished at high volume fractions due to the presence of the wall, which is there referred to as ‘screening out’. In this paper we qualitatively confirm this observation. However, we provide a more quantitative method to determine the near-wall self-diffusion coefficients and we are able to show the diminishment of HI affects the self-diffusivity normal and parallel to the interface to a different extent. This becomes possible using the virial approximation for the initial decay rate of the scattered electric field autocorrelation function (EACF), which we described in our earlier contribution.³³ There we presented a detailed discussion of the derivation of exact expressions for the first cumulant (i.e. initial decay rate) of the EACF in a concentrated suspension of hard spheres. After constructing an appropriate theoretical framework based on the Smoluchowski equation, we have presented two methods for practical calculations of the first cumulant: the virial expansion, and precise multipole simulations. While the latter may be used for high accuracy calculations at any volume fraction of the suspension, the virial expansion is expected to correctly reproduce the experimentally measured cumulants up to moderate concentrations.

The aim of this paper is to present our results on near-wall dynamics in a model hard-sphere system, viewed in EWDLS experiments. For the first time we provide a thorough analysis of the scattering vector dependence of the first cumulant, which allows us to assess the range of volume fractions where the virial approximation can be used to describe the experimental data. By tuning the suspension volume fraction and the penetration depth, we are able to investigate the effect of HI-diminishment for high concentrations and its anisotropy in a more convincing way than proposed earlier.^30,31 Comparing to virial expansion results and simulations, we are able to assess the validity of the former approximate scheme at higher volume fractions. We also discuss in detail the colloidal near-wall self-diffusivity which may be determined from our results in a similar way as proposed for the corresponding bulk property by Pusey,³⁴ Segré³⁵et al. and Banchio et al.³⁶ and analysed theoretically by Abade et al.³⁷

The paper is organised as follows. After a short review of the theoretical foundations (Section 2) and the details of the numerical simulations (Section 4), we describe the details of sample preparation, the evanescent wave light scattering setup and the data analysis in Section 5. The experimental EWDLS data are compared to the theoretical predictions in the result and discussion Section 6 where we also confront our predictions to data published earlier and discuss the progress we make here beyond the state of earlier contributions.^30,31 Finally we summarize our results in the conclusion Section 7.

2 Theoretical description

We consider an ensemble of N identical spherical particles of radius a immersed in a Newtonian solvent of viscosity η. The fluid is bounded by a planar no-slip wall at z = 0.

In EWDLS experiments, the scattered light intensity time autocorrelation function g₂(t) is measured, from which the normalised scattered electric field correlation function ĝ₁(t) (EACF) is calculated. Since the scattered electric field E_s depends on the configuration of the system, i.e. the positions of the particles, its fluctuations can be related to the diffusive dynamics of near-wall particles. The initial decay of the EACF is exponential in time

with Γ being the first cumulant, similarly to bulk DLS.³⁸ However, there are two important differences to the bulk case. Firstly, the sample is illuminated by a non-uniform evanescent wave. Its intensity decays exponentially with the distance z from the wall as exp(−κz), thus restricting the scattering volume to a wall-bounded region with a thickness of order κ⁻¹. The particles staying closer to the boundary receive more intensity and yield the strongest signal. The instantaneous scattered electric field is then given for an ensemble of N particles as³³where r_j is the position of the centre of sphere j, Q is the scattering vector and z_j = r_j·ê_z, with ê_z being a unit vector normal to the wall.

Secondly, the mobility of the particles is strongly hindered by the presence of the wall. The boundary reflects the flow incident upon it, leading to an increase of friction, and thus a slow-down of colloidal dynamics. The effect is more pronounced for particles staying close to the surface where their mobility becomes distance-dependent. This information is encoded in the hydrodynamic mobility tensor μ^w_ij which describes the velocity U_i the particle i acquires due to the force F_j applied to the particle j

U_i = μ^w_ij·F_j. (3)

For non-interacting spheres in a wall-bounded fluid, the tensors μ^w_ij become diagonal in particle indices, but retain the anisotropic structure which follows from the invariant properties of the system,

μ^w_ij = δ_ij[μ^w_∥(1 − ê_zê_z) + μ^w_⊥ê_zê_z]. (4)

where 1 is the unit tensor, and μ^w_∥,⊥ are scalar mobilities for motion parallel and perpendicular to the boundary. In the absence of the wall, the mobility tensor becomes isotropic, with μ_∥ = μ_⊥ = μ₀ = 1/6πηa being the Stokes mobility of a spherical particle. It follows from the fluctuation-dissipation theorem that the Stokes–Einstein diffusion coefficient D₀ is given by k_BTμ₀, where k_B is the Boltzmann constant, and T denotes the temperature. The same relation holds between the many-particle diffusion matrix D and the mobility tensor μ^w.

Using the Smoluchowski equation formalism, Cichocki et al.³³ derived an analytical expression for the first cumulant measured in an EWDLS experiment for a suspension of spherical particles,

where the hydrodynamic function readsand the wall-structure factor is given by Here, nA/κ is the number of particles within the illuminated scattering volume, with n being the bulk particle number density, and A is the illuminated wall area. The brackets 〈⋯〉 denote ensemble averaging. Eqn (5) is a generalisation of the bulk result for concentrated suspensions³⁹which corresponds to the limit of infinite penetration depth or κ → 0.

Decomposing the scattering vector into components parallel and perpendicular to the wall,

Q = Q_∥ + Q_⊥ = Q_∥ê_∥ + Q_⊥ê_z, (9)

where ê_∥ is a unit vector in the direction of Q_∥ and using the invariant properties of the system, we arrive at the following structure of the first cumulantwhere

H_⊥ = ê_z·H_w(κ,Q)·ê_z, (11)

H_∥ = ê_∥·H_w(κ,Q)·ê_∥, (12)

H_I = ê_z·2Im[H_w(κ,Q)]·ê_∥, (13)

H_R = ê_∥·2Re[H_w(κ,Q)]·ê_z, (14)

with Im and Re standing for the imaginary and real part, respectively. The coefficients H as well as S_w may be either evaluated numerically using the virial expansion approach, or by extracted from numerical simulations. Both techniques are briefly described in the course of this work. The expressions given above are valid for an arbitrary wall–particle interaction potential. Further on, we restrict to hard-core interactions.

In the dilute regime, the hydrodynamic function and structure factor have only single-particle contributions, from which it follows that H_I = H_R = 0. The surviving parts D₀H_∥/S_w and D₀H_⊥/S_w in eqn (10) simplify then to the single-particle average diffusion coefficients 〈D_∥〉_κ and 〈D_⊥〉_κ, respectively, in agreement with the notation proposed in earlier works.^17–19 In the case of hard-core sphere-wall interactions, the penetration-depth average (in the dilute limit) reads

We may now explicitly write the first cumulant in this case as^18,19 The averaged diffusivities 〈D_∥,⊥〉_κ have been calculated as functions of κa in ref. 17.

Importantly, this is also the case in the limit of Q_∥ → ∞ or Q_⊥ → ∞, where only the self-parts of H_w and S_w survive. The cumulant may then be expressed using the self-diffusion tensor D^s which is defined as the initial slope of the mean square displacement tensor of a tracer particle located at a height z at t = 0, viz.

where Δr(t) is the displacement vector of the tracer particle during the time t. The tensor D^s may be expressed in terms of the mobility matrix μ^w, as we have shown in ref. 33. Thus, the cumulant may be approximated for sufficiently large Q_∥ or Q_⊥ bywhere D^s_∥,⊥ are the components of the self-diffusion tensor D^s(z), with the average given byand g(z) being the single-particle distribution function. Its definition readswhere P^w_eq(R) is the equilibrium probability density function (in the presence of a wall) for the system to be at a configuration R = {r₁,…,r_N}. The quantities in eqn (6), (7), and (20), are taken in the thermodynamic limit, which has been discussed for a wall-bounded system in ref. 33. In a dilute system with hard sphere-wall interactions, and when interactions between the particles may be neglected, the average (19) reduces to the formula (15).

3 Virial expansion

For moderately concentrated systems, calculations of the wall structure factor S_w and the components of the wall-hydrodynamic tensor H_w may be performed by expanding them in terms of powers of bulk particle concentration n far from the wall. The procedure has already been elaborated in great detail in ref. 33. Thus, we refrain here from the technical aspects, focusing on the resulting expressions.

The small dimensionless parameter in the density expansion is the bulk volume fraction,

instead of the concentration n. The virial expansion of the wall-structure factor (7) reads The coefficient S⁽¹⁾ and the self-part of S⁽²⁾ may be found analytically as

S⁽¹⁾ = e^−κa, (23)

with the particle diameter d = 2a. The distinct part of S⁽²⁾ has to be evaluated numerically. The analogous virial expansion of the wall hydrodynamic tensor requires a cluster decomposition of the mobility matrix,⁴⁷ and has a similar form In this case in order to obtain the terms H⁽¹⁾ and H⁽²⁾ we need the one- and two-particle cluster components of the mobility matrix. Explicit expressions for S^(1),(2) and H^(1),(2) are rather complex, and have been given explicitly in ref. 33. To calculate them, the hydromultipole code, implemented according to ref. 48, has been used.

Inserting the expansions (22) and (25) into eqn (5), we find the following virial expansion for the first cumulant

where the factoris the infinite dilution prediction, given explicitly by eqn (16), while the second term reads These virial expansion results, together with simulations that are also valid at high concentrations, will be compared to experiments in Section 6. The relations above may be transformed using eqn (11)–(14) and expressed in terms of the tensorial components of the hydrodynamic function H_w. The subsequent section contains the details of simulations.

4 Numerical simulations

To determine the equilibrium wall-structure factor (7), the hydrodynamic functions (11)–(14), and the first cumulant (5), we have carried out a series of numerical simulations for a wall-bounded hard-sphere system with particle volume fractions in the range 0 < ϕ ≤ 0.3. Key elements of our numerical techniques are summarized below; a more detailed description is provided in our previous paper.³³

Since hydrodynamic-interaction algorithms are unavailable for a single-wall system with periodic-boundary conditions, the calculations were performed for a suspension confined between two well separated parallel walls. The equilibrium particle distributions were determined using a standard Monte-Carlo (MC) algorithm, and the multiparticle mobility matrix μ^w_ij was evaluated using the periodic version⁴⁰ of the Cartesian-representation algorithm^41–43 for a suspension of spheres in a parallel-wall channel.

Most of our calculations were carried out for a wall separation h = 13d (where d is the sphere diameter). By comparing results for different values of h, we have established that the above wall separation is sufficient to obtain accurate one-wall results, provided that the particle volume fraction is adjusted for the excess particle density in the near-wall regions.

To evaluate the required volume-fraction correction, we constructed the equilibrium ensemble for a reference system with a large wall separation h = h₀ and the assumed particle number density n in the middle of the channel. The excess particle number per unit area, n_ex, was determined using the formula

N = Ahn + 2An_ex, (29)

where N is the number of particles in the periodic cell, A is the wall area, and h = h₀ is the wall separation in the reference system. The particle number N = N(h) for channels with different widths h is obtained from expression (29), with known reference values of n and n_ex.

Since the evanescent wave scattering occurs only near the illuminated surface, and the hydrodynamic field associated with the periodic forcing ∼exp(iQ·r) decays on the length scale l ∼ Q_∥⁻¹ with the distance from the wall, the effect of the second wall of the channel on the multiparticle mobility is small.⁴² We find that for the evanescent wave parameters corresponding to our experiments, the effect of the second wall on the hydrodynamic functions (11)–(14) is smaller than the statistical simulation inaccuracies.

The hydrodynamic tensor H_w was determined as an average over M independent MC trials. To obtain statistical accuracy of the order of 2%, we have used M in the range from M = 30 for large systems with N ≈ 10³ particles to M = 400 for N ≈ 200 particles.

5 Experimental details

5.1 Hard-sphere sample and preparation

As model systems for the EWDLS experiments, we used two batches of poly(methyl methacrylate) (PMMA) particles, named ASM470 and ASM540 in the following, which were purchased from Andrew Schofield, University of Edinburgh. The spherical particles are covered with a thin poly-12-hydrostearic acid layer to stabilize them against aggregation in organic solvents. To allow scattering experiments at high volume fractions the particles were transferred from a cis-decaline suspension (as received) to a refractive index matching cis-decaline/tetraline mixture by spinning and re-dispersing them. The solvent used had a cis-decaline mass fraction of w = 0.2, a refractive index of n₂ = 1.498 and a viscosity of η = 2.658 mPa s at temperature of T = 298 K as measured using an Abbemat RXA156 and an Automated Microviscometer AMVη from Anton Paar, Graz, Austria. To determine the particle radius, we employed standard Dynamic Light Scattering (DLS) measurements. The recorded time autocorrelation functions of the scattered intensity g₂(t) (IACF) were analysed by three different methods, namely cumulant analysis, stretched exponential fitting and inverse Laplace transformation. The three methods yield hydrodynamic radii of R_H = 98 nm (ASM470) and R_H = 144 nm (ASM540) with a variation of less than 1 nm in both cases. These values are assumed to be identical with the hard sphere particle radius a in the following. Further, the size distributions obtained from inverse Laplace transformation showed a full width at half maximum of less than five percent. The negligible size polydispersity is confirmed by the observation that the suspensions crystallize at sufficiently large particle volume fractions.

Prior to the scattering experiments, the suspensions were filtered through PTFE syringe filters with a nominal pore size of 1 μm directly into the measurement cells to minimize parasitic scattering from dust particles. To reduce the number of necessary alignment processes of the EWDLS measuring cell, this was filled with the hard sphere suspension of highest volume fraction, and further dilution was achieved by removing a part of the sample and replacing it by pure solvent. The exact volume fraction was determined a posteriori by drying a 250 μl aliquot and determining the mass of the remaining particles. Further, the EWDLS sample cell was equipped with a small magnetic stirrer bar with which the samples were homogenised before each angular scan to minimize the influence of particle sedimentation.

5.2 EWDLS set-up

EWDLS experiments were performed with a home-built instrument, based on a triple axis diffractometer by Huber Diffraktionstechnik, Rimsting, Germany, which has been described in detail elsewhere.¹⁹ The setup is equipped with a frequency doubled Nd/Yag Laser (Excelsior; Spectra Physics) with a vacuum wavelength of λ₀ = 532 nm and a nominal power output of 300 mW as a light source. Scattered light is collected with an optical enhancer system by ALV Lasververtriebsgesellschaft, Langen, Germany, which is connected to two avalanche photo diodes by Perkin Elmer via an ALV fiber splitter. The TTL signals of the diode were cross-correlated using an ALV-6000 multiple tau correlator. The scattering geometry and the definition of the scattering vector and its component parallel and normal to the interface are sketched in Fig. 1. The sample cell (custom-made by Hellma GmbH, Muellheim, Germany) consists of a hemispherical lens as the bottom part, made of SF10 glass, with an index of refraction n₁ = 1.736 at λ₀ = 532 nm. The hard sphere suspension is contained in a hemispherical dome sitting on top of the lens. The primary beam is reflected off the interface between the glass and the solution, by that creating an evanescent wave in the solution which is used as the illumination for the scattering experiment. The evanescent wave has a wave vector k_e pointing in the direction of the reflected beam's projection onto the reflecting interface. The scattering vector is given by Q = k_s − k_e, where the scattered light wave vector, k_s, is defined by the two angles θ and α_r which describe the position of the detecting unit.

Fig. 1 Scattering geometry in EWDLS setup with a spherical geometry. Left: For the definition of angles and wave vectors. Right: For the illustration of the primary beam being back reflected at exit of the hemispherical lens to air, thereby creating a second evanescent wave.

The inverse penetration depth of the evanescent wave depends on the angle of incidence α_i as . The magnitudes of the scattering vector components parallel and normal Q_⊥ = 2πn₂ sin α_r/λ₀ to the interface can be varied by changing θ and α_r. In a typical EWDLS experiment, we record series of correlation functions at fixed Q_⊥ varying Q_∥ (Q_∥-scan) or vice versa (Q_⊥-scan).

5.3 Details of data analysis

The analysis of the scattered intensity time autocorrelation function g₂(t) from EWDLS is much less straightforward than in conventional bulk dynamic light scattering (DLS), mainly for two reasons. The first major complication occurs from the fact that a simple quadratic Siegert relation between g₂(t) and the correlation function of the scattered field g₁(t) which is usually assumed in DLS does not apply in most cases in EWDLS. As described in Section 5.2, the incident laser beam is totally reflected from the glass/solution interface in EWDLS. As this interface is always corrugated, it acts as a static scatterer which in general contributes significantly to the observed signal. Therefore a mixed homodyne/heterodyne detection scheme has to be taken into account, and the generalized Siegert relation⁴⁴for the conversion from g₂(t) to ĝ₁(t) has to be applied. Here, and C₁ = C₂ − C₂², with A being the experimental intercept of g₂(t). Further, in many cases EWDLS intensity-autocorrelation functions exhibit a very slow decay at large times. The physical origin of this slow relaxation is not clear yet. While Garnier et al. conjecture that it is due to a slow reversible adsorption of the particles to the wall due to van der Waals attraction,¹² Steffen⁴⁶ and Lisicki et al.¹⁷ argue that it is also caused by the unavoidable stray-light from surface defects, which is scattered by colloids in the bulk of the suspension into the detector. Since these slow modes are in general well separated from the relaxation rates of interest, we approximate their contribution by an additional baseline B₁ to ĝ₁(t).

Thus, to determine the initial slope Γ of ĝ₁(t), which is related to the dynamic properties of interest, we chose to non-linear least squares fit the experimental correlation functions to eqn (30), where g₁(t) is modelled as a decaying single exponential function in time

ĝ₁(t) = (1 − B₁) exp{−Γt} + B₁. (31)

According to eqn (30) and (31), B₁ is related to B₂, the baseline of g₂(t), by . Consequently there are three fit parameters A, Γ and B₂. Since an erroneous baseline value will lead to a deviation of Γ from its true value, due to a normalization error, we fitted g₂(t) repeatedly, starting with a number of data points, N_p. After a single fit had converged, two data points at the long time end of g₂(t) were removed, reducing N_p by two, and the remaining data points were fitted again. This procedure was repeated until N_p < 20. With this technique it was possible to identify a limited range of N_p's where the values of B₂ and Γ are essentially independent of N_p. The Γ values determined in this range are considered to be the initial slope or the first cumulant of ĝ₁(t). Where error bars are presented with values of Γ, they reflect the standard deviation of repeated measurements.

5.4 Effect of back-reflection

An additional difficulty in EWDLS stems from the fact that, differently from bulk DLS, it is not possible to apply a refractive index-matching batch around the sample cell without fundamental changes of our instrument design. Therefore the primary beam will be back-reflected at the exit from the semi-spherical lens with a reflectance R. In the present case the semi-spherical lens has a refractive index of n₁ = 1.736, which leads to a reflectance of R = 0.072 according to Fresnel's equations.⁴⁵ As sketched in Fig. 1, the back-reflected beam will also be reflected off the glass sample interface, thereby causing a second evanescent wave with wave vector k_R = −k_P, where k_P is the wave vector of the evanescent wave caused by the original primary beam. In what follows, the subscript P will refer to the evanescent wave caused by the primary beam, while R will be associated with the evanescent wave due to the back-reflected beam. The latter gives rise to a second scattering process, for which the in-plane scattering angle is θ_R = 180 − θ_P. Consequently the scattering vector components parallel to the interface are given bywhere i ∈ {R,P}. Differently, the component normal to the interface remains unchanged in the two cases The normalized field correlation functions in such a situation should be considered as a weighted sum of two individual correlation functions from two scattering experimentswhere , and P(Q_i) is the particle scattering factor of a sphere.

In evanescent illumination, the scattering factor is affected by the non-uniform character of the electric field and becomes penetration-depth dependent. For an optically uniform particle, the scattering amplitude in the evanescent field reads . Thanks to the high symmetry, for a spherical particle of radius a, B(Q,κ) can be explicitly calculated as

with . The particle scattering factor is then found as P(Q) = |B(Q,κ)|².

In order to illustrate the effect of back-reflection for a dilute suspension, we analyse eqn (34) using the field correlation function given by the first two cumulants:

with the first cumulant Γ_i given in the dilute regime by eqn (16), and the second cumulant can be calculated as

Γ_2,i = γ_i − Γ_i², (37)

where the second moment γ_i is defined as¹⁷ Here . The resulting IACF has to be calculated from ĝ₁(t) using the generalized Siegert relation, eqn (30). The averaged diffusion coefficients, which are required for the calculations of Γ and Γ₂ at a given value of κa were calculated in ref. 17.

The calculated model correlation functions are now evaluated according to the same analysis procedure as the experimental data to obtain initial slopes Γ as a function of Q_∥. In the top part of Fig. 2 we compare initial relaxation rates of model correlation function, which were calculated in this way at infinite dilution. The model calculations coincide perfectly at low scattering vectors. However, those data, which were calculated taking into account the effect of back-reflection, strongly decrease at larger scattering vectors. Here and in the following, we will present the results in dimensionless form, i.e. relaxation rates in units of D₀/d² and scattering vectors in units of 1/d where D₀ is the particles' bulk diffusion coefficient at infinite dilution.

Fig. 2 Top: Result of model calculations for zero particle density without (dashed dotted line) and with (full line) taking into account the effect of back-reflection, showing a considerable difference in the high-Q range. Bottom: Comparison of experimental data obtained at ϕ = 0.175 (full circles) to simulations (line with triangles) and virial approximations at the same concentrations with (full line) and without (dashed dotted line) taking into account the effect of back-reflection. The experimental parameter Q_⊥d = 3.15 and κd = 2.08 are the same for both graphs.

For a concentrated suspension, both Γ_P and Γ_R may be calculated from the virial expansion and from simulations. The first cumulant of the EACF including the back-reflection effect may thus be written from eqn (34) as

again without any free parameter. In the bottom part of Fig. 2 the same experimental data are compared to simulation results and to virial calculations for ϕ = 0.175. In both, the virial calculations and the simulations, the effect of back-reflection can be included as described above. It turns out that up to Q_∥d ∼ 7 the first term in eqn (39) dominates, so that Γ ≈ Γ_P, and the back-reflection effect need not be taken into account. However, in the high-Q range, the back-reflection is essential to correctly reproduce the experimental data, as seen in Fig. 2. We are therefore led to conclude that the first cumulants obtained experimentally at high in-plane angles, i.e. θ > π/2 should be considered with extreme care and potentially discarded when comparing experimental data to theoretical predictions and simulations.

6 Results and discussion

To illustrate the influence of the particle volume fraction, we display experimental data of Γ versus the scattering vector from Q_∥-scans with ASM470 suspensions at different volume fractions in Fig. 3. It is obvious that at high Q_∥ the experimental data deviate from the virial approximation displayed as full lines in Fig. 3 for all concentrations, which is fully explained by the effect of back-reflections, discussed in Section 5.4. Apart from this high-Q_∥ deviation, the virial approximation predicts the experimentally observed data correctly, even at a sphere volume fraction of almost 25 percent.

Fig. 3 Relaxation rates versus parallel component of the scattering vector. Symbols represent experimental data obtained from ASM470 (R_H = 98 nm) suspensions at different volume fractions at Q_⊥d = 3.0 and κd = 2.6, lines are prediction by the virial approximation and open symbols refer to data points which are obscured by the back-reflection effect discussed in Section 5.4.

The same degree of agreement between virial approximation and experimental data is observed in Q_⊥-scans, which are shown in Fig. 4. Here we display the experimental data obtained from the ASM540 suspension with ϕ = 0.175 performed at the same penetration depth κd = 2.08 but with extremely different values of parallel scattering vector component, i.e. Q_∥ = 1.83 and Q_∥ = 5.7. Together with these data we also present the results of simulations which were obtained for a set of similar parameters, i.e. Q_∥d = 1.83, κd = 2.08 and ϕ = 0.15. At low Q_∥, the results from all three methods agree very well, and at large Q_∥, where no simulation data are available, the agreement between experiment and virial approximation is also within the experimental error.

Fig. 4 Relaxation rates versus normal component of the scattering vector. Symbols represent experimental data obtained from an ASM540 (R_H = 144 nm) suspension with ϕ = 0.175 at κd = 2.08 at different values of the parallel scattering vector component, i.e. Q_∥d = 5.7 (open circles) and Q_∥d = 2.18 (full squares). Full lines are predictions by the virial approximation for the same experimental parameters and the dashed-dotted line refers to simulation results obtained for Q_∥d = 2.18 and κd = 1.8.

Only at the highest volume fraction (ϕ = 0.3) for which experimental data and simulations are available there is a significantly better agreement between simulation data and experiments than between virial approximation and measured data. This is shown in Fig. 5 where we display data from a Q_∥ scan, obtained from an ASM540 suspension with ϕ = 0.3 at Q_⊥d = 2.36 and κd = 2.08 together with the corresponding predictions. At this high volume fraction the deviation between virial approximation and simulations is comparable or even larger than experimental error bars.

Fig. 5 Relaxation rates versus parallel component of the scattering vector. Symbols represent experimental data obtained at Q_⊥d = 2.36 and κd = 2.08 from an ASM540 (R_H = 144 nm) suspension with ϕ = 0.30. The full lines are prediction by the virial approximation and the line with triangles refers to simulation results. Experimental data points, which are obscured by the back-reflection effect discussed in Section 5.4, are omitted in this graph.

It is interesting to investigate the limit of large scattering vectors, where self-diffusion is probed. As discussed in Section 5.4, the relaxation rates determined at the largest scattering vectors (and thus the largest angle θ) are not reliable. However, as first suggested by Pusey,³⁴ self-diffusion in bulk can be probed approximately at a wave vector Q* such that S(Q*) ≈ S(Q → ∞). This observation has been supported theoretically by Abade et al.³⁷ It is expected that at this point the distinct structure factor vanishes, and likewise does the distinct hydrodynamic function, so that only the self-parts contribute to the dynamic properties at this point. In the bulk case, this statement was later corroborated by extensive numerical simulations.^35,36 Segré et al.³⁵ stated that in a bulk suspension of hard spheres, this point is found for Q*a ∼ 4.0, where S(Q*a) = 1, to the right of the main peak of S(Q). As shown by Banchio et al.,⁵² bulk structure factors of hard sphere suspensions with different volume fractions show an isosbestic point at S(Qa = 4.02) = 1 and at the same value of Qa the corresponding hydrodynamic functions attain their high-Q limit. Michailidou et al.³⁰ used the EWDLS experimental data at Qa = 4.58 arguing that this should not be too far from Q*a, thus providing a good estimate of the near-wall self-diffusion coefficient. Here, we propose a more thorough way to determine the particles' near-wall self-diffusion properties which follows the same line of arguments as discussed for bulk systems above. We note here that in EWDLS both the structure factor and the hydrodynamic function become penetration-depth dependent.³³ However, upon re-scaling by their asymptotic values, both S_w(Q) and the components of H_w(Q) exhibit an isosbestic point at which they attain their asymptotic values. We compute them using the virial expansion, and plot the results in Fig. 6. Like for bulk experiments, first cumulants obtained at the Q_∥,⊥a values of the isosbestic point provide a good approximation for the near wall self-diffusion coefficients.

Fig. 6 The wall structure factor S_w(Q) and the hydrodynamic function H_∥(Q) in a Q_∥-scan at fixed Q_⊥d = 2.36 and κd = 2.08 for a selection of volume fractions. Both functions are normalised by their self-values at Q → ∞. At (Q_∥d)² ≈ 53 we find an isosbestic point for both functions (marked by the dotted vertical line), suggesting that the self-diffusion coefficients may be determined from the data collected in the vicinity of this point. The statement also holds sway for other components of H.

However, as the first isosbestic point right of the structure factor main maximum is found approximately at Q_∥*d = 7.3, we could determine experimental data of the first cumulant at this scattering vector mainly from the ASM540 suspensions. For the smaller ASM470 particles the data at Q_∥*d is distorted by the back-reflection effect (except for ϕ = 0.1), and thus it may not be used to experimentally determine the self-diffusion coefficient parallel to the wall 〈D^s_∥〉_κ. Further, the experimentally accessible range of Q_⊥d is in all cases much smaller than Q_⊥*d such that we can not get reliable experimental information on the self-diffusion properties normal to the wall.

In Fig. 7, we present the normalized ratios of 〈D^s_∥〉 over the bulk self-diffusion constant. The latter was calculated according to the semi-empirical formula⁴⁹

which includes two virial coefficients due to Batchelor⁵¹ and Cichocki et al.,⁴⁷ and is expected to be accurate up to ϕ ∼ 0.45. Its validity has been extended by Riest et al.⁵⁰ up to ϕ = 0.5 by modifying the coefficient of the last term to −0.70. We compare experimental data to predictions by virial approximation and simulations. The theoretical values for 〈D^s_∥〉_κ were determined by linearly extrapolating the high-Q range of the Γ vs. Q_∥² dependence, making use of eqn (18). Our experimental data confirm the trend predicted by both methods and show that the near-wall dynamics approach the bulk behaviour at high particle volume fractions. With this observation we qualitatively confirm the earlier results by Michailidou et al.^30,31 Since the virial approach allows quick calculation of Γ vs. Q_∥² data, we can easily predict near wall self-diffusion coefficients for a variety of parameters, by using the slope in the high Q-range. We use this possibility to quantitatively compare self-diffusion properties predicted by the virial approximation to the data by Michailidou. For this purpose we calculate 〈D^s_∥〉_κ and 〈D^s_⊥〉_κ for a series of volume fractions and average them as 〈D^s_w〉_κ = (〈D^s_∥〉_κ + 〈D^s_⊥〉_κ)/2 according to their experimental procedure. Their choice of Qa = 4.58 is determined by the fact that they measured with a geometry which corresponds to θ = 0° and α_r = 90°, thus at a scattering vector which makes an angle of 45° with the interface. In this configuration the parallel contribution and the normal contribution to self-diffusivity are weighted equally in the experiment. The comparison in Fig. 8 shows that the prediction calculated by a 1 : 1 weighing of the normal and the parallel component are deviating systematically from the experimental data in the range of volume fractions, where the virial approach should hold. Only at very high volume fractions, where the virial approximation is certainly not valid the experimental data appear to agree with it. This is probably due to the effect that first cumulants obtained at Qa = 4.58 are not a good approximation for the self-diffusion properties. Actually simulations of bulk properties⁵² show that even at moderate volume fractions, both the structure factor and the hydrodynamic function are significantly different from their value at Q*a = 4.02. For the sake of completeness we also show predictions for the self-diffusion constants in Fig. 8, which are averaged according to 〈D^s_w〉_κ = (2〈D^s_∥〉_κ + 〈D^s_⊥〉_κ)/3. These agree reasonably well with the earlier experimental data, which is probably a coincidence.

Fig. 7 Comparison of experimental data (full circles) obtained from ASM470 (R_H = 98 nm, ϕ = 0.1) and ASM540 (R_H = 144 nm) suspensions, virial calculations (full line, the dashed dotted line represents virial calculations in a range of volume fractions where the approximation is not considered valid) and simulation results (triangles) for the self-diffusion coefficient parallel to the wall. Experimental parameters are at Q_⊥d = 2.36 and κd = 2.08 for all cases.

Fig. 8 Comparison of virial predictions for the self-diffusion coefficient with experimental data by Michailidou et al.³⁰ measured at κd = 0.89. The predicted data for normal and parallel contribution were averaged as indicated in the legend.

Nevertheless, we confirm the earlier conjecture that particle–particle hydrodynamic interactions in the presence of a wall are diminished at high volume fractions as compared to bulk dynamics. However, here we can show that effect influences the diffusion parallel to the wall and normal to the wall differently. As discussed in ref. 33, the anisotropic self-diffusion coefficients have the following virial expansion

The coefficients of this expansion have been presented graphically in Fig. 3 of ref. 33. Here we have tabulated them for a selection of penetration depths in Table 1. The coefficients have a clear interpretation: G⁽¹⁾_∥,⊥ refers to single-particle dynamics at infinite dilution, while G⁽²⁾_∥,⊥ bears information on the effect of the wall on two-particle interactions. All coefficients decrease with increasing κd, but the effect is stronger for the motion perpendicular to the wall. The behaviour of G⁽¹⁾_∥,⊥ follows from the single-particle physical picture,¹⁷ in which motion normal to the interface is suppressed more than in the parallel direction. This is due to the fact that perpendicular motion generates ‘squeezing’ flows which lead to stronger hydrodynamic resistance as compared to ‘shearing’ flows induced by parallel motion.⁵³ The particle–particle HI are affected in the same way, which explains the faster decay of G⁽²⁾_⊥ as compared to G⁽²⁾_∥. Thus, the coefficients corresponding to the normal motion are affected more strongly. However, the near-wall self-diffusivity is frequently written in the formwith 〈D_∥,⊥〉_κ = D₀G⁽¹⁾_∥,⊥(κd). The coefficientbecomes a result of an interplay between the single- and two-particle effects. In Fig. 9 we show normalized ratios of 〈D^s_∥〉_κ and 〈D^s_⊥〉_κ over the bulk self-diffusion as a function of volume fractions for two different penetration depths of the evanescent wave. The curves are calculated using the virial approach up to a volume fraction of 25%. First we observe that the self-diffusion coefficient (averaged over the illumination profile) normal to the wall is smaller than that parallel to the wall and that both components increase with penetration depth κ⁻¹, similarly to the components of the near-wall diffusion coefficients at infinite dilution. The variation of these ratios over the range of volume fractions covered is indicated by the numbers on the far right of Fig. 9, which are the ratios of the values obtained at ϕ = 10⁻³ and ϕ = 0.25. It is important to note, that although 〈D^s_∥〉_κ/D^s_b varies stronger with increasing ϕ as compared to 〈D^s_⊥〉_κ/D^s_b, this does not imply that the wall diminishes the particle–particle HI more in the parallel direction, as we discussed above.

Table 1 The coefficients of the virial expansion of anisotropic self-diffusivity, defined in eqn (41). The decay of the ⊥ elements is faster with increasing penetration depth, indicating that both single- and two-particle mobilities are hindered more for motion in the direction normal to the interface

κd	G ^(1)⊥(κd)	G ^(2)⊥(κd)	G ^(1)∥(κd)	G ^(2)∥(κd)
0	1.0	−1.832	1.0	−1.832
0.2	0.781	−1.371	0.884	−1.535
0.5	0.644	−1.117	0.810	−1.357
1.0	0.516	−0.871	0.736	−1.160
2.0	0.383	−0.588	0.654	−0.903
5.0	0.227	−0.250	0.547	−0.550

---

Fig. 9 Calculated data for the self-diffusion coefficients parallel and normal to the wall for different penetration depths: κd = 2.08 full lines and κd = 1.3 dashed dotted lines. The numbers on the far right represent the ratio of the values at ϕ = 10⁻³ over ϕ = 0.25, which are an indication that the diminishment of hydrodynamic interaction is more pronounced for particle motion parallel to the wall than normal to the wall.

7 Conclusions

In this paper we describe our EWDLS investigations of the near wall dynamics of colloidal hard spheres in suspensions with volume fractions up to ϕ = 0.3. We thoroughly compare experimental data for the dependence of the first cumulant on the scattering vector components parallel and normal to the interface to corresponding predictions based on a second order virial approximation and to simulation results, where the full hydrodynamic interaction is taken into account. Up to volume fractions of about fifteen to twenty percent we find perfect agreement between the three methods. Above this range, the predictions by the virial approach deviate discernibly from the simulation data,³³ however this deviation is still in the range of experimental error bars. Therefore we conclude that the virial approach provides a good approximation for the prediction and analysis of experimental data up to a volume fraction of about 25%, which is much less time consuming and elaborate than full scale simulations. Only at ϕ ≥ 0.3 the virial approximation is clearly not anymore able to capture the details of the dependence of the first cumulant on the scattering vector. Further we introduce a new method to assess the particles' near wall self-diffusivity from experimental data. This method follows the same line of argument, which is used to assess bulk self-diffusivity in cases where the limit of sufficiently high scattering vector cannot be reached experimentally. We identify an isosbestic point of the near-wall structure factors right to the first maximum, where near wall structure factor and hydrodynamic function attain their asymptotic values. Diffusion data measured at the scattering vector of the isosbestic point are a good approximation for the self-properties. Comparison of experimental data with predictions, based on the virial approach and on simulations, show that this method yields better estimates of the self-diffusivity as methods used earlier. Finally we confirm earlier data which show that the diminishment of particle–particle hydrodynamic interactions due to the presence of the wall is less pronounced at high volume fraction compared to bulk dynamics. Beyond that, we show (see Table 1) that the observed effect is weaker for the mobility parallel to the wall as compared to motion in the normal direction. In conclusion, with the virial approximation, we have a method at hand, which qualitatively supports earlier data, but provides significant further insight into the near wall dynamics of colloidal hard spheres. This is especially important since this approach can be easily adopted to systems with long ranging static interaction, providing a quick and non-costly method for the prediction and analysis of EWDLS results obtained from e.g. charged colloids.

Acknowledgements

YL thanks for support from Marie Sklodowska Curie initial Training network SOMATAI under EU Grant Agreement No. 316866. JB was supported by NSF Grant No. CBET-1059745. He would also like to acknowledge the hospitality and financial support from IPPT PAN during his summer visits. ML wishes to acknowledge support from the National Center of Science grant no. 2012/07/N/ST3/03120. EW was supported by the National Center of Science grant no. 2012/05/B/ST8/03010. The work was supported by the Foundation for Polish Science International PhD Projects Programme co-financed by the EU European Regional Development Fund. Part of the experimental data presented was obtained with financial support from the European Commission under the Seventh Framework Program by means of the grant agreement for the Integrated Infrastructure Initiative No. 262348 European Soft Matter Infrastructure (ESMI) which is gratefully acknowledged.

References

1. K. Meister, S. Strazdaite, A. L. De Vries, S. Lotze, L. C. Olijve, I. K. Voets and H. J. Bakker, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 17732.
2. M. Hans, A. Erbe, S. Mathews, Y. Chen, M. Solioz and F. Mücklich, Langmuir, 2013, 29, 16160.
3. R. H. Tromp, M. Vis, B. H. Erné and E. M. Blokhuis, J. Phys.: Condens. Matter, 2014, 26, 464101.
4. P. Lang, J. Phys.: Condens. Matter, 2004, 16, R699 , and references therein.
5. D. Kleshchanok, R. Tuinier and P. R. Lang, J. Phys.: Condens. Matter, 2008, 20, 073101 , and references therein.
6. H. A. Lorentz, Abhandlungüber Theoretische Physics, Leipzig und Berlin, 1907.
7. H. Faxén, Ark. Mat., Astron. Fys., 1923, 17, 1.
8. Z. Adamczyk, M. Adamczyk and T. G. M. van de Ven, J. Colloid Interface Sci., 1983, 96, 204,  DOI:10.1016/0021-9797(83)90022-X.
9. J. Happel and H. Brenner, Low Reynolds Numbers Hydrodynamics, Kluwer, Dordrecht, 1991.
10. R. Sigel, Curr. Opin. Colloid Interface Sci., 2009, 14, 426.
11. K. H. Lan, N. Ostrowsky and D. Sornette, Phys. Rev. Lett., 1986, 57, 17.
12. N. Garnier and N. Ostrowsky, J. Phys. II, 1991, 1, 1221.
13. N. Ostrowsky and N. Garnier, Prog. Colloid Polym. Sci., 1991, 84, 371.
14. N. Ostrowsky and L. Lobry, Nuovo Cimento D, 1994, 16, 1111.
15. M. I. M. Feitosa and O. N. Mesquita, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 44, 6677.
16. H. Matsuoka, H. Morikawa, S. Tanimoto, A. Kubota, Y. Naito and H. Yamaoka, Colloid Polym. Sci., 1998, 276, 349.
17. M. Lisicki, B. Cichocki, J. K. G. Dhont and P. R. Lang, J. Chem. Phys., 2012, 136, 204704,  DOI:10.1063/1.4720069.
18. P. Holmqvist, J. K. G. Dhont and P. R. Lang, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 021402.
19. P. Holmqvist, J. K. G. Dhont and P. R. Lang, J. Chem. Phys., 2007, 126, 044707.
20. B. Loppinet, G. Petekidis and G. Fytas, Langmuir, 1998, 14, 4958.
21. G. Fytas, S. H. Anastasiadis, R. Seghrouchni, D. Vlassopoulos, J. Li, B. Factor, W. Theobald and C. Topracksioglou, Science, 1998, 274, 2041.
22. G. E. Yakubov, B. Loppinet, H. Zhang, J. Rühe, R. Sigel and G. Fytas, Phys. Rev. Lett., 2004, 92, 115501.
23. P. Holmqvist, D. Kleshchanok and P. R. Lang, Langmuir, 2007, 23, 12010.
24. M. Haghighi, M. N. Tahir, W. Tremel, H.-J. Butt and W. Steffen, J. Chem. Phys., 2013, 139, 064710.
25. S. A. Rogers, M. Lisicki, B. Cichocki, J. K. G. Dhont and P. R. Lang, Phys. Rev. Lett., 2012, 109, 098305.
26. M. Lisicki, B. Cichocki, S. Rogers, J. K. G. Dhont and P. R. Lang, Soft Matter, 2014, 10, 4312.
27. B. Loppinet, J. K. G. Dhont and P. R. Lang, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 62.
28. Z. Li, L. D'eramo, C. Lee, F. Monti, M. Yonger, P. Tabeling, B. Chollet, B. Bresson and Y. Tran, J. Fluid Mech., 2015, 766, 147.
29. A. Stocco, T. Mokhtari, G. Haseloff, A. Erbe and R. Sigel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 011601.
30. V. N. Michailidou, G. Petekidis, J. W. Swan and J. F. Brady, Phys. Rev. Lett., 2009, 102, 068302,  DOI:10.1103/PhysRevLett.102.068302.
31. V. N. Michailidou, J. W. Swan, J. F. Brady and G. Petekidis, J. Chem. Phys., 2013, 139, 164905,  DOI:10.1063/1.4825261.
32. J. W. Swan and J. F. Brady, J. Chem. Phys., 2011, 135, 014701.
33. B. Cichocki, E. Wajnryb, J. Bławzdziewicz, J. K. G. Dhont and P. R. Lang, J. Chem. Phys., 2010, 132, 074704,  DOI:10.1063/1.3305328.
34. P. N. Pusey, J. Phys. A: Math. Gen., 1978, 11, 119.
35. P. N. Segré, O. P. Behrend and P. N. Pusey, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 5070,  DOI:10.1103/PhysRevE.52.5070.
36. A. J. Banchio and G. Nägele, J. Chem. Phys., 2008, 128, 104903,  DOI:10.1063/1.2868773.
37. G. C. Abade, B. Cichocki, M. L. Ekiel-Jeżewska, G. Nägele and E. Wajnryb, J. Chem. Phys., 2010, 132, 014503.
38. B. J. Berne and R. Pecora, Dynamic Light Scattering: With Applications to Chemistry, Biology, and Physics, Dover Books on Physics Series, Dover Publications, 2000.
39. A. Z. Akcasu and H. Gurol, J. Polym. Sci., Polym. Phys. Ed., 1976, 14, 1.
40. J. Bławzdziewicz and E. Wajnryb, Phys. Fluids, 2008, 20, 093303.
41. S. Bhattacharya, J. Bławzdziewicz and E. Wajnryb, J. Fluid Mech., 2005, 541, 263.
42. S. Bhattacharya, J. Bławzdziewicz and E. Wajnryb, Physica A, 2005, 356, 294.
43. S. Bhattacharya, J. Bławzdziewicz and E. Wajnryb, J. Comput. Phys., 2006, 212, 718.
44. J. K. G. Dhont, An Introduction to Dynamics of Colloids, in Studies in Interface Science, ed. D. Möbius and R. Miller, Elsevier, 1996.
45. M. Born and E. Wolf, Principles of optics, Cambridge Univ. Press, 6th edn, 1980.
46. W. Steffen, MPI Mainz, Germany, personal communication, 2012.
47. B. Cichocki, M. L. Ekiel-Jeżewska and E. Wajnryb, J. Chem. Phys., 1999, 111, 3265,  DOI:10.1063/1.479605.
48. B. Cichocki, R. B. Jones, R. Kutteh and E. Wajnryb, J. Chem. Phys., 2000, 112, 2548,  DOI:10.1063/1.480894.
49. G. C. Abade, B. Cichocki, M. L. Ekiel-Jeżewska, G. Nägele and E. Wajnryb, J. Chem. Phys., 2011, 134, 244903.
50. J. Riest, T. Eckert, W. Richtering and G. Nägele, Soft Matter, 2015, 11, 2821,  10.1039/c4sm02816c.
51. G. K. Batchelor, J. Fluid Mech., 1976, 74, 1.
52. A. J. Banchio, M. G. McPhie and G. Nägele, J. Phys.: Condens. Matter, 2008, 20, 404213,  DOI:10.1088/0953-8984/20/40/404213.
53. S. Kim and S. J. Karrila, Microhydrodynamics: Principles and Selected Applications, Butterworth-Heinemann, Boston, 1991.

---