The distribution of a macromolecular solute within an evaporating drop: an exact analytical solution

Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

The distribution of a macromolecular solute within an evaporating drop: an exact analytical solution Th. Klupsch a,∗ , P. Mühlig a , R. Hilgenfeld a,b a

b

Department of Structural Biology and Crystallography, Institute of Molecular Biotechnology, Beutenbergstrasse 11, P.O. Box 100813, D-07708 Jena, Germany Institute for Biochemistry, University of Lübeck, Ratzeburger Allee 160, D-23538 Lübeck, Germany Received 29 January 2003; accepted 13 June 2003

Abstract A space- and time-dependent concentration profile of a low-diffusive solute within an evaporating droplet will be formed due to the forced inward motion of the solute at the confining, inward-moving free surface (solute entrainment), if the solvent evaporates and any convection within the droplet can be neglected. To study this effect theoretically, a generalized diffusion equation including the solute entrainment is derived. An exact analytical solution for spherical and semi-spherical droplet geometries is presented assuming a constant solute diffusion coefficient (linear problem) and an ideal evaporation kinetics of the Langmuir type. Analytical approximations are outlined, and conditions are given to apply the results in practice where the ideal evaporation kinetics is disturbed. Starting from a constant solute concentration, a new spherical–symmetrical concentration profile with maximum at the droplet boundary is approached during evaporation, which is scaled by the time-dependent droplet radius only, and which is independent of the starting droplet radius for small deviations from the Langmuir evaporation kinetics. A concentration enhancement at the droplet boundary by a factor of about 2 compared with the mean concentration level is possible for macromolecular solutes under conditions as used, e.g., in evaporating droplet arrangements for protein crystal growth. In addition, the solute entrainment can act to increase a starting concentration variation parallel to the droplet surface if the droplet kinetics exceeds a certain threshold value. Since any convective flux within the droplet has to be excluded in order to apply our results, the convection-free state has to be proved by observations taking into account that convective instabilities in evaporating droplets are hardly to predict. © 2003 Elsevier B.V. All rights reserved. Keywords: Macromolecular solution; Solute entrainment; Solute concentration profile; Evaporating droplet; Diffusion with moving boundary

1. Introduction In this contribution, modern aspects of droplet physics, namely hydrodynamics, diffusion and evapo∗ Corresponding author. Tel.: +49-3641-206105; fax: +49-3641-206199. E-mail address: [email protected] (Th. Klupsch).

ration kinetics will meet actual problems of biophysics and biochemistry. In particular, we will present an analytical solution for the concentration profile of a low-diffusion coefficient solute within an evaporating droplet. It is a prerequisite that the solvent evaporates only, and that any convective flow within the dropled can be neglected.

0927-7757/$ – see front matter © 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.colsurfa.2003.06.002

86

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

In the following, a solute will be diffusing weakly if the diffusion is too slow to balance the forced inward motion of the solute at the confining, inward-moving free surface (solute entrainment). For this case an equilibration of the concentration across the bulk of the droplet is not reached during evaporation. As we will show in more detail, such a behaviour is typical for macromolecular solutes. Thus, assuming a constant concentration at the beginning of the evaporation kinetics, this profile will undergo a successive deformation with a maximum at the moving droplet boundary, followed by a monotonic decrease until the minimum at the midpoint of the droplet is reached. The maximum may be significantly higher than the mean concentration of the droplet, while the minimum will scarcely exceed the starting concentration level, depending on the special choice of parameters governing the problem. Droplet techniques are widely used in biophysics and biochemistry. For instance, the static and dynamic shape analysis of droplets hanging from a capillary has been developed into a sensitive method to determine the concentration-dependent surface tension of biomolecules at the air–water interface (see, e.g. [1]). This is fundamental to study of the structure and interaction of two-dimensional arrays of macromolecules. Also, the so-called evaporating droplet techniques are established methods for growing protein single crystals from aqueous solutions (see, e.g. [2]) for the purpose of structural analysis by X-ray diffraction. Recently, measurements of thermodynamic bulk properties of protein solutions from evaporating droplets have been reported [3]. It is a requirement for all these methods, in principle, to have a control of the solute concentration. In situ monitoring and control of the spatially averaged lysozyme concentration by means of fibre-optic Raman spectroscopy during crystallization in a hanging drop has recently been achieved [4]. Experimental and theoretical studies of the space-dependent concentration profile of weakly diffusing solutes in evaporating droplets have not been reported so far, although the technique to create concentration gradients of a macromolecular solute by solvent evaporation and to study the phase behaviour by interferometric methods is established for larger-sized cells, where the recording sensitivity was recently improved [5].

Let us discuss the problem of the concentration profile of weakly diffusing solutes in more detail for crystal growth by means of evaporating droplet techniques [2,4,6–9] (vapour diffusion methods). Their sole aim of these experiments is to obtain sufficiently large and defect-free single crystals. To achieve this goal it is suggested, from general principles of crystal growth (see, e.g. [10]), to choose the “reaction path” such that the nucleation as well as the growth of crystallites during the ripening stage takes place under conditions sufficiently close to the thermodynamic equilibrium between crystal and solution. Here, the term “reaction path” may be taken to summarize the time-dependent course of the concentration of all solute components, namely that of the protein con(e) centration c(t) and of the concentrations ci (t), i = 1, 2, . . . of the ith electrolytic component (measured in mole per volume) [2,4]. The solvent (water) will evaporate as long as its equilibrium vapour pressure (which equals its vapour pressure at the droplet surface) exceeds the ambient and adjusted solvent vapour pressure within the crystallization cell, far from the droplet. Then the evaporation kinetics act to increase the concentrations of the solute components and to decrease the equilibrium solvent vapour pressure until the ambient value is approached and the evaporation kinetics stop [2,4,6–9]. The part of the reaction path from start up to the appearance of the first crystallites is determined by the evaporation kinetics only, because no solute molecules are taken from the solution due to crystal growth [4,9]. If there is no further knowledge of what happens within the droplet, all solute concentrations should be assumed, as a first approach to this stage, to be uniquely scaled with the time-dependent droplet volume V(t) according to c(t) ≈ c0

V(0) , V(t)

(e) V(0)

(e)

ci (t) ≈ ci0 (e)

(e)

V(t)

(1)

where c0 = c(0), ci0 = ci (0) are the starting concentrations of the solute components at t = 0. A well-chosen reaction path should start in the undersaturated region and arrive, after crossing the solubility curve, at a moderately supersaturated state before crystallization sets in. We now take into account the fact that the macromolecular solute component may be weakly diffusing. Obviously, (1) does no longer hold because c becomes

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

space- and time-dependent, c = c(r, t). Then, starting from a homogeneous concentration in the undersaturated regime, first a region at the droplet boundary becomes supersaturated. This boundary region now governs the success of any crystallization attempt. The nucleation sets in first near the droplet boundary, as it can be concluded from the rapid decrease of the induction period with increasing supersaturation (see, e.g. [4,11,12]). The crystallites from these first nucleation events will encounter during ripening better conditions for a rapid growth than the remaining crystallites because of their larger size and the enhanced solute concentration in their surroundings (see, e.g. [13]). Therefore, the knowledge of c(r, t) is the base for estimating or modelling optimum crystallization conditions. On the other hand, the thermodynamic data required for a model to include, for instance, the complicated interaction between the macromolecular and electrolytic solute are missing for most of proteins to be crystallized. Therefore, the philosophy of the widely established screening methods is to find appropriate starting conditions of the reaction path by trial and error (see, e.g. [14]). Noting now that the weakly diffusing properties of the macromolecular solute may strongly affect the crystallization conditions, a strict control of all parameters which may influence the evaporation kinetics, and which are commonly treated to be outside the scope of screening is recommended in order to warrant the reproducibility and to cause crystallization attempt to be successful. To get a more qualified understanding of various questions mentioned, we will define and solve the underlying boundary value problem of the diffusion equation for c(r, t) with solute entrainment at moving boundaries. We will analyze a class of solutions derived for two idealized arrangements: (a) the free spherical droplet and (b) the hemispherical droplet fixed on an infinite plane, respectively, both coming into contact, for t ≥ 0, with an infinite reservoir with a given solvent vapour pressure. When the evaporation kinetics start at t = 0, a constant solute concentration c(r, 0) = c0 = const is assumed. The concentration-dependence of the solute diffusion coefficient will be neglected. The related boundary value problem can be solved exactly treating a limiting case of the evaporation kinetics which was first considered by Langmuir and others [15–17], and which will be called the “Langmuir kinetics” in the following. It

87

would occur if (i) the dependence of the solvent vapour pressure on the solute concentration; (ii) the decreasing droplet temperature due to the heat of evaporation; and (iii) the influence of the surface tension on the evaporation kinetics could be neglected. Such conditions are approached for sufficiently small solute concentrations and sufficiently slow evaporation kinetics, if the droplet radius is not too small (minimum about 100 ␮m). Nevertheless, the exact solution given here can be extended, in practice, for conditions which significantly deviate from that of the Langmuir kinetics. However, it is assumed that there is no convective flux in the droplet which would strongly disturb the solute concentration profile. Below we will discuss how to exclude this worst case, where also a condition is given to otherwise justify the application of the formulae. The basic problem described so far is strictly spherical–symmetric; that means the droplet preserves its spherical shape during evaporation, and c(r, t) = c(r, t), with r as the radial coordinate measured from the midpoint of the droplet. A generalization may be of interest assuming that for a spherical droplet with a Langmuir kinetics, the starting solute concentration is slightly space-dependent, c(r, 0) = c0 (r, 0) = c0 (r, ϑ, ϕ) (r, ϑ, ϕ denote spherical coordinates with the origin at the midpoint of the droplet). This may be the case, e.g., if macromolecular fringes (“schlieren”) are initially created in the drop. Then a space and time-dependent concentration is expected according to c(r, t) = c(r, ϑ, ϕ, t). While the solute entrainment always acts to enhance the radial-symmetric part of the concentration profile, there is the question whether a starting tangential concentration-dependence (in the direction parallel to the droplet surface) diminishes due to self-diffusion. We will show that such a tangential balancing only occurs for sufficiently small evaporation rates, meanwhile there exists a critical threshold above which a lateral solute concentration fluctuation increases. This finding is inconsistent with the silent assumption of a constant surface tension and a constant solvent vapour pressure at the surface, as required to preserve the spherical droplet shape during evaporation. Therefore, this threshold may be identified with the onset of, to the best of our knowledge, a new type of instability of the droplet evaporation kinetics which acts to deform the spherical droplet shape. The driving force is the tangential concentration variation at the droplet surface, which

88

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

causes an “anisotropic” evaporation kinetics due to the space-dependence of the outwardly directed diffusive solvent flux at the droplet surface. Of course, this effect is supported by the surface tension which becomes space-dependent, too. Therefore we suppose that this instability may establish a non-spherical, steady-state droplet shape formed through a complicated interplay of non-constant contributions of the surface curvature, surface tension and solute concentration along the droplet surface. We will restrict ourselves to the calculation of the instability threshold without going into details of the droplet kinetics beyond this threshold. This paper is organized as follows. In Section 2, more details of the Langmuir kinetics are given as well as basic parameters entering the calculations are introduced. The exact solution of the basic boundary value problem as well as the instability threshold due to enhanced lateral solute concentration fluctuations is outlined in Section 3. The discussion given in Section 4 focuses on disturbances of the concentration profile created by the solute entrainment, which includes a condition for applying the derived formulae for a weakly disturbed Langmuir kinetics. A mathematical appendix is added where the rigorous formulation of the boundary value problem, its exact solution for the Langmuir kinetics, and its extension to non-Langmuir droplet kinetics applying methods of the singular perturbation theory is presented.

2. Idealized droplet kinetics and basic parameters For the Langmuir kinetics, the diffusion current of the solvent vapour of both the arrangements (a) and (b) introduced above is strictly radial. Therefore, the time-dependent evaporation is governed, for both cases, by the same law, and the time-dependent droplet radius R0 (t) obeys the Langmuir differential equation [15–17] (g)

dR0 (t) 1 (g) ρ = −(1 − f)DA A(l) dt ρ R0 (t)

(2)

A

(g)

(g)

(l)

f, DA , ρA , and ρA denote the relative humidity in the cell (normalized to unity), the solvent–vapour diffusion coefficient, the solvent–vapour density in

thermal equilibrium, and the liquid–solvent density, respectively (Langmuir [15] treated only the case f = 0, but his result is easily generalized to arbitrary f (see, e.g. [16,17])). Eq. (2) is derived assuming that a quasi-static solvent vapour diffusion field accompanies the droplet, which is justified because (g) of |dR0 (t)/dt|  DA /R0 (t), which holds true in practice. Eq. (2) is solved by
t R0 (t) = R0 (0) 1 − 2 (3a) τK where R0 (0) is the starting droplet radius, and the characteristic kinetic time constant τ K is (l)

τK = R0 (0)2

(g)

(ρA /ρA )

(g)

(1 − f)DA

(3b)

Obviously, τ K /2 is the life-time of the evapo(g) rating droplet. With the numerical data DA = (g) (l) 2.49 × 10−5 m2 /s, ρA = 5.70 × 1023 m−3 , ρA = 28 −3 3.34 × 10 m for water at 293 K, and R0 (0) = 7 × 10−4 m, f = 0.85, we find τK = 7.66 × 103 s. It should be noted that, in practice, the decrease of the droplet temperature during evaporation can diminish the evaporation kinetics so that τ K -values may be somewhat higher than calculated. We now introduce the diffusion time τD =

R0 (0)2 (l)

DP

(4)

as the second characteristic parameter. It denotes the order of the time required to achieve a homogeneous concentration profile of the macromolecular solute by diffusion within a droplet with the constant radius (l) R0 (0). DP is the diffusion coefficient of the macro(l) molecular solute. For a rough estimation of DP , we may assume spherically shaped globular protein molecules with unique mass densities so that the hydrodynamic radius is nearly proportional to M1/3 (M: molecular weight). Then from the well-known (l) Stokes–Einstein relation (see, e.g. [18]), DP becomes −(1/3) approximately scaled with M . Inserting known data [19] we arrive, for water as solvent, at an approximate scaling relation valid at room temperature (l) DP ≈ αM −(1/3) , α = 3.24 × 10−9 m2 /s. For M = (l) 1 × 105 Da, we find DP ≈ 3.2 × 10−11 m2 /s, which 4 yields τD ≈ 1.53 × 10 s for R0 (0) = 7 × 10−4 m.

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

Finally we introduce the dimensionless parameter β¯ defined by τD (5) β¯ = − τK where β¯ < 0. As will be shown below, β¯ is the only parameter which governs the long-time behaviour of the concentration profile. From the numerical example given above, we find that β¯ may vary over two orders of magnitude, from about −0.1 to −10, depending on the conditions of the evaporation kinetics (relative humidity) and the molecular weight of the macromolecular solute. The effects discussed here become important if β¯ < −0.5.

3. Exact solution and instability threshold

In this section we will present first the exact solution for the initial condition of a constant starting concentration, that means c(r, t) = c0 = const

for t = 0 and 0 ≤ r ≤ R0 (0) (6)

For this case the solution becomes spherically symmetrical. We introduce the time-dependent scaling factor a(t) of the droplet radius and the stretched radial coordinate x according to
t (7) a(t) = 1 − 2 τK x=

r r = R0 (t) a(t)R0 (0)



= c0 1 + β¯

 ∞   ρn n=0

σn

 (1 − a(t)

σn /β¯

¯ )Φ0 (x, σn , β) (9b)

(10a)

¯ ≥ 0 if β¯ ≤ 0, σ0 (β) ¯ < 0, if β¯ ≤ 0 σn (β)

and

n≥1

¯ > σn+1 (β), ¯ σn (β) ¯ → −∞, if n → ∞ and β¯ < 0 σn (β)

(10b)

(10c)

¯ are given by The functions Φν (x, σ, β)  ¯ 2 σ ν 3 βx ν ¯ Φν (x, σ, β) = x 1 F 1 − + , + ν, − 2 2 2 2β¯ (11) where 1 F 1 is the confluent hypergeometric function (see, e.g. [20]). For ν = 0, β¯ ≤ 0, the basis functions are approached by

(8)

so that the bulk of the droplet is given by 0 ≤ x ≤ 1. The exact solution reads   r ,t (9a) c(r, t) = c¯ R0 (t) c¯ (x, t)

¯ σn = σn (β) ¯ are certain coefficient where ρn = ρn (β), functions defined and outlined, in detail, in Appendix A and B (ρn should be distinguished from the density (g) ρA , etc.). Eq. (9b) involves a series expansion with respect to certain basis functions (“normal modes”) ¯ β), ¯ n = 0, 1, 2, . . . of the stretched space Φ0 (x, σn (β), coordinate x, and time-dependent coefficients where time is involved only implicitly via the scaling function a(t). As a peculiarity, the stretched space coordinate x is time-dependent via a(t), too. The unusual power law of the time-dependent coefficients in (9b) is typical of boundary conditions at moving boundaries; it changes to a common exponential form by carrying out the transform to the stretched coordinates in space and time (see Appendix A). In particular, we have ¯ = −3β¯ σ0 (β)

3.1. General solution for a constant starting concentration

89

¯ β) ¯ ∼ Φ0 (x, σ0 (β), =

 1/2 2 ¯ 2 /4) ¯ −(1/2) e−(βx (−β) 3 ¯ sinh Q0 (x, β) (12a) × x

 ¯    βx 6 3 2 ¯ Q0 (x, β) = − x − + ¯ 4 2 β    β¯ × Ar sinh x −  6

(12b)

90

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

¯ β) ¯ Φ0 (x, σn (β), −(1/2)  ¯ sin Qn (x, β) 3β¯ ¯ 2 ∼ ¯ e−(βx /4) , = −σn (β) − 2 x if n ≥ 1 (13a)  1/2 3β¯ ∼ ¯ ¯ Qn (x, β) = −σn (β) − 2   3 2 x β¯ × x+ , if 0 ≥ β¯ ≥ −15 ¯ + (3β/2) ¯ 24 σn (β) (13b) ¯ β) ¯ are normalized The basis functions Φ0 (x, σn (β), according to ¯ β) ¯ → 1, Φ0 (x, σn (β),

if x → 0

(14)

¯ β) ¯ The first member of the basis functions, Φ0 (x, σn (β), with n = 0, is exponentially increasing in x, while for n ≥ 1, the basis functions are oscillating. The course of the first few basis functions is shown in Fig. 1. 3.2. Long-time behaviour The asymptotic behaviour of c(r, t) for large t is determined solely by the term n = 0 in the sum of (9b), as follows from (10a) and (10b). It is given by

¯ plotted vs. x = r/R0 (t) for various ¯ 0 (x, β) Fig. 2. The function Φ ¯ describes the local concentration enhancement related ¯ Φ ¯ 0 (x, β) β. to the mean concentration in the large-time limit.

  ¯ Φ0 (x, σ0 (β), ¯ β) ¯ ρ0 (β) c¯ (x, t) = c0 + O(1) 3 a(t)3   c0 ¯ ¯ 0 (x, β), Φ if t → ∞ = a(t)3

(15)

where we have introduced ¯ = 1 ρ0 (β)Φ ¯ 0 (x, σ0 (β), ¯ β) ¯ ¯ 0 (x, β) Φ 3 1 ¯ 0 (x, −3β, ¯ β) ¯ = ρ0 (β)Φ 3

¯ β), ¯ Fig. 1. The first members of the basis functions Φ0 (x, σn (β), n = 0, 1, 2, taken for β¯ = −2, and plotted vs. x = r/R0 (t).

(16)

According to (15), the solute concentration profile becomes scaled, for sufficiently large t, with the scaling factor a(t)−3 . Because c0 /a(t)3 = c0 (V(0)/V(t)) is the mean solute concentration within the evaporating ¯ describes the asymptotic behaviour ¯ 0 (x, β) droplet, Φ of the space-dependent, relative solute concentration enhancement with respect to the spatially averaged profile. It is governed by the one parameter β¯ only. ¯ is plotted, for different β, ¯ in Fig. 2. If β¯ < 0, ¯ 0 (x, β) Φ ¯ ¯ we have Φ0 (1, β) > 1 at the droplet boundary, x = ¯ < 1 at the droplet midpoint, ¯ 0 (0, β) 1, and 0 < Φ x = 0. This means that the solute concentration at the droplet boundary increases above its average value, as expected, but it remains below its average at the ¯ and Φ ¯ is monotoni¯ 0 (1, β) ¯ 0 (0, β) droplet midpoint. Φ ¯ cally increasing and decreasing with |β| = τK /τD , re-

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

91

we get m = 4πR0 (0)3 a(t)3

1

dx x2 c¯ (x, t)



0



 1 4π ¯ 0 (x, σ0 (β), ¯ β) ¯ = R0 (0)3 c0 dx x2 ρ0 (β)Φ 3 0  (19) + O(a(t)3 ) = const Because (19) is valid for arbitrary large t where the remainder can be neglected, we find 1 ¯ ¯ β) ¯ ρ0 (β) dx x2 Φ0 (x, σ0 (β), 0

¯ and Φ ¯ Φ ¯ is the ¯ 0 (1, β) ¯ 0 (0, β): ¯ 0 (1, β) Fig. 3. The functions Φ concentration enhancement related to the mean concentration at ¯ is the concentration enhancement ¯ 0 (0, β) the droplet boundary; Φ related to the mean concentration at the droplet midpoint.

spectively, as shown in Fig. 3. Both the functions can be calculated using the approximants

=3

(17a)

¯ ∼ ¯ 0 (0, β) Φ = 0.9808 + 0.298β¯ + 0.0342β + 0.00141β¯ 3 , if − 10 ≤ β¯ ≤ −0.25

(17b)

¯2

In practice, (17a) together with (15) can be applied to find the solute concentration at the droplet boundary if a(t) = R0 (t)/R0 (0) ≤ 0.7. Quite similarly, the solute concentration at the droplet midpoint is found from ¯ ¯ 0 (0, β). (17b) and (15) if a(t)3 < Φ 3.3. Mass conservation and transient behaviour The total solute mass m within the evaporating droplet is conserved, that means   4π R0 (0)3 c0 m= 3 R0 (t) = 4π dr r2 c(r, t) = const (18) 0

Changing to the stretched radial coordinate x using (7)–(9a), and inserting theasymptotic relation (15),

¯ =1 ¯ 0 (x, β) dx x2 Φ

(20)

0

Inserting now (9b) into the first part of (19) and taking into account (20), we get

 1 ∞  ¯ ρ ( β) n 2 ¯ β) ¯ dx x 1 + β¯ Φ (x, σn (β), ¯ 0 σn (β) 0 − β¯

¯ ∼ ¯ 0 (1, β) Φ = 1.016 − 0.134β¯ + 0.00056β¯ 2 − 0.00009β¯ 3 , if − 10 ≤ β¯ ≤ −0.25

1

n=0

1 ¯ ρn (β) ¯ ¯ ¯ β)=0 ¯ a(t)σn (β)/β dx x2 Φ0 (x, σn (β), ¯ σ (β) 0 n=1 n

∞ 

(21) Again taking into account (10a) and (20) and observing that (21) must be fulfilled for arbitrary t, we arrive at a second relation, in addition to (20): 1 ¯ β) ¯ = 0, if n ≥ 1 dx x2 Φ0 (x, σn (β), (22) 0

Inserting now (9a) and (9b) into (18) we find   R0 (t) c0 ¯ m = 4π ρ0 (β) dr r2 a(t)3 0   r ¯ ¯ × Φ0 , σ0 (β), β R0 (t) 0

R0 (t)

 dr r Φ0 2

 r ¯ ¯ , σn (β), β = 0, R0 (t)

(23a)

if n ≥ 1 (23b)

The meaning of (23a) and (23b) is obvious. Although the solute concentration profile c(r, t) appears as superposition of distinct space- and time-dependent

92

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

δc(r, ϑ, ϕ, 0) may be sufficiently small so that, at least during any initial period, a Langmuir kinetics may be assumed. Provided there exists a threshold ¯ > |β¯ c |, the fluctuavalue β¯ = β¯ c so that for |β| tion δc(r, ϑ, ϕ, t) increases for the asymptotic large ¯ < |β¯ c |. Because t-limit, while it vanishes for |β| δc(r, ϑ, ϕ, t) may be represented as a series expansion in spherical harmonics, we can restrict our treatment to the simplified ansatz, which in analogy to (6), (9a) is (µ)

δc0 = δc(r, ϑ, ϕ, 0) = δc0ν (r)Pνµ (cos ϑ)eiµϕ , if 0 ≤ r ≤ R0 (0) δc(r, ϑ, ϕ, t) = δ¯cν(µ)



(24)

 r , t Pνµ (cos ϑ)eiµϕ , R0 (t) if 0 ≤ r ≤ R0 (t) (25)

µ

Fig. 4. The concentration enhancement related to the mean concentration, c(r, t)/a(t)3 , at the transitional stage as it develops from a homogeneous starting concentration, plotted vs. x = r/R0 (t). (l) Parameters: R0 (0) = 8 × 10−4 m, DP = 3 × 10−11 m2 /s, 4 4 τD = 2.13 × 10 s, τK = 1.06 × 10 s corresponding to β¯ = −2. Curve 1: t = 8.0 × 102 s corresponding to a(t) = 0.92. Curve 2: t = 1.6 × 103 s corresponding to a(t) = 0.83. Curve 3: t = 3.2 × 103 s corresponding to a(t) = 0.63. Curve 3 is very close to the asymptotic profile for β¯ = −2 shown in Fig. 2.

contributions according to (9b), the total solute mass m involves only the normal mode with n = 0, as Eq. (23a) illustrates. All other modes in the sum in (9b) can be understood to describe “solute mass redistribution effects”, which do not contribute to the total solute mass, according to (23b). But the latter govern, of course, the transient behaviour from the constant concentration at the start towards the asymptotic, exponential concentration profile. In Fig. 4, an example of the transient concentration profile is shown. 3.4. An instability threshold of the Langmuir kinetics Generalizing the problem outlined so far, we are now interested in the behaviour of any superposed angular-dependent solute concentration fluctuation δc = δc(r, ϑ, ϕ, t). At the starting time, δc0 =

where the Pν (cos ϑ) are the associated Legendre functions (ν, µ integer, ν > 0, |µ| ≤ ν). Then, the asymptotic behaviour for large t is found to be ¯

δc(r, ϑ, ϕ, t) ∝ a(t)σν0 /β   r ¯ × Φν , σν0 , β Pνµ (cos ϑ)eiµϕ , R0 (t)

if t → ∞ (26)

¯ is given by (11), and σν0 = σν0 (β) ¯ where Φν (x, σ, β) are certain coefficient functions obeying ¯ > σ10 (β) ¯ > σ20 (β) ¯ > ··· σ0 (β)

(27)

¯ is monotonically increasing In particular, σν0 (β) ¯ starting from σν0 (0) < 0 for ν > 0. The with |β|, ¯ is changed at β¯ = β¯ cν , where |β¯ c1 | < sign of σν0 (β) ¯ and σ20 (β) ¯ are |β¯ c2 | < · · · . The functions σ10 (β) shown in Fig. 5. According to (26), the change of the ¯ indicates the transition of δc(r, ϑ, ϕ, t) sign of σν0 (β) from the decreasing to the increasing behaviour for large t. Obviously, the threshold β¯ c is given by β¯ c = β¯ c1 , which means that the part of starting perturbation δc(r, ϑ, ϕ, 0) with the lowest symmetry (ν = 1) governs the large-time behaviour of the perturbation. From numerical calculations we find β¯ c1 = −1.303

(28)

We note that this is the necessary condition only for the Langmuir kinetics to become destabilized against

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

¯ for ν = 1, 2, 3, plotted vs. Fig. 5. The coefficient functions σν0 (β) ¯ The amplitude of any concentration fluctuation ∝ Pν(µ) (cos ϑ) β. ¯ > 0. increases, due to the solute drag, if σν0 (β)

non-spherical deformations of the droplet shape according to the mechanism proposed in Section 1. The question whether this instability is strong enough to be observed in practice can be answered only by carrying out a complete stability analysis taking into account the coupling with the evaporation kinetics, the solvent convection, and the surface tension, which is outside the scope of this paper.

4. Discussion 4.1. Strong convective disturbances The physics of evaporating droplets is, in reality, much more complicated than prescribed by the model of the idealized Langmuir kinetics. First, we have to discuss the occurrence of effects which strongly disturb or destroy a solute concentration profile created by diffusion. But if this can be excluded, there remain several “weak” perturbations which lead to more or less strong deviations from the Langmuir kinetics. Therefore, subsequent conditions are required to guarantee that the formulae outlined above can be applied, in practice.

93

Any convection within the evaporating droplet will destroy the concentration profile due to diffusion if the convective velocity uc reaches or exceeds the order dR0 /dt, which is typical 10−7 m/s. Two effects are known in the literature, which can create such a convection. The first one is the Marangoni–Bénard instability [21,22] where the driving force is the gradient of the surface tension created by thermal or compositional inhomogeneities. The second effect is responsible for the formation of the so-called “coffee ring pattern” which appears for sessile droplets with a pinned contact line [23,24]. The appearance of both the effects leading to strong disturbances is difficult to predict. For the second one, the convection results from a geometrical constraint only. uc is of the order 10−7 m/s if the initial contact angle starts at π/2 and significantly decreases during evaporation. The contact line pinning may be caused by the stochastic roughness of the substrate surface [23,24]. With respect to the Marangoni instability, some investigations and estimations for droplets have been carried out [21,22] although a more detailed mathematical analysis for the spherical geometry including calculations of critical Marangoni numbers are to the best of our knowledge missing so far. Characteristic sizes of the Bénard cells and convection velocities of about 10−8 m and 10−1 m/s, respectively, have been reported [22] for investigations of the surfactant-enhanced Marangoni instability on sessile droplets of aqueous solutions, with droplet volumes of about 7 × 10−9 m3 , and total evaporation times of the order of 1000 s corresponding to a relative humidity of about 50%. The Marangoni convection can superpose the convective flow due to contact line pinning [23,24]. The Marangoni instability in droplets is expected to depend very sensitively upon geometrical peculiarities and initial conditions. For instance, the temperature gradient between the evaporating droplet boundary and the substrate in sessile droplet arrangements is strongly affected by the heat transport through the substrate [24]. Otherwise, for free droplets under steady-state conditions, the heat of evaporation would be completely balanced through the heat flux from the exterior assuming, for the moment, the droplet radius to be constant, which means that the temperature gradient within the droplet is only related to the time-dependent change of the droplet radius

94

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

(assuming the droplet radius to be constant, the heat flux from outside would, after a short transitional period, completely balances the loss due to the heat of evaporation, see also [25]). Similar conditions are met, e.g., for droplets hanging on a fine wire loop [21] or on a fine capillary [25]. The critical Marangoni number [21] (see also [22] and references sited therein) can be estimated only on the basis of numerical values of governing material parameters (temperature-derivative and the compositional derivative of the surface tension) which are not known for most of the protein solutions of interest. The segregation of macromolecules and fine particles at the droplet surface may significantly influence the convection [24]. The kinetics of the formation of 2D ordered structures of macromolecules at the surface (see, e.g. [1]) may act in the same way. The enhancement of the Marangoni–Bénard convection due to surfactants as well as its suppression due to the adsorption of trace surface-active impurities is reported for aqueous systems [22]. 4.2. Weak perturbations Assuming now that there is no convection within the droplet, there remain in practice a lot of “weak” perturbations which affect the idealized Langmuir evaporation kinetics and the solute concentration profile depending on them. We may roughly distinguish between effects due to geometrical restrictions and thermodynamic perturbations. The former are related to deformations of the idealized spherical–symmetrical diffusion field of the solvent vapour around the droplet (see, e.g. [6,8]). For droplets positioned, for example, in a closed, narrow container, they become significant, if the distance between the droplet midpoint and the container wall, and the droplet midpoint and the surface of the reservoir that is introduced to stabilize the relative humidity, respectively, is less than about ten droplet diameters. The thermodynamic disturbations are connected with the change of the solvent vapour pressure at the droplet surface during evaporation, as it results from the decrease of the droplet temperature and the enhanced solute concentration (see, e.g. [2,6,26]), respectively. In particular, we observed that the strong temperature decrease that occurs in evaporating droplets, which consist of aqueous solutions, and which hang from a fine capillary, can cause

a buoyancy-driven convection of the surrounding air. This affects the droplet kinetics [25]. In practice, direct observation using optical microscopy seems the only realistic possibility to prove the absence of convection within the evaporating droplet. It is recommended that such observations be combined with the measurement of the evaporation kinetics in order to quantify the deviations from Langmuir kinetics. Some arrangements for carrying out such combined investigations have been applied to sessile droplets, for example, using video technique [23,24] and microscopy together with weighing [22]. Droplets hanging from a capillary have also been studied in this way [1,21,25]. The usual technique used to identify convection cells consists of the observation of patterns formed by dispersed, fine particles on the droplet surface [21–24] or, after evaporation, on the substrate [22]. This technique is demonstrated, for evaporating droplets hanging from a capillary, in Fig. 6a and b. Details will be given elsewhere. The formulae derived for the Langmuir kinetics, outlined above, can be extended to describe the concentration profile in a spherical or semi-spherical droplet in the presence of weak perturbations. Instead of complicated calculations of the droplet kinetics, the experimental data of the time-dependent droplet radius R0 (t) for a certain observation period 0 ≤ t ≤ tobs will be taken as input, because not all effects entering the droplet kinetics may be precisely known. In particular, R0 (t) may be given as a sufficiently smooth function where rapid fluctuations due to stochastic errors of measurement are averaged out. This works provided we introduce the new time-dependent functions  2 ˜ = 1 d R0 (t) β(t) (29a) 2 dt R0 (0) ˜ ¯ = τD β(t) β(t)

(29b)

¯ ˜ = which reduce to β(t) = −1/τK = const, β(t) −τD /τK = const for the Langmuir kinetics. As shown in Appendix A, the formulae derived for the Langmuir kinetics remain valid for weak disturbations if the scaling function a(t) given by (7) is substituted by  ˜ (30) a(t) = 1 − 2β(t)t and the parameter β¯ is everywhere substituted by the ¯ given by (29b). In addition, the following function β(t)

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

condition must be satisfied:      (dβ(t)/dt)  ¯ ˜ 2 3|β(t)| σ1 β(t)  ˜  , , ,    min ˜  β(t)  t a(t)2 τD a(t)2 0 < t ≤ tobs

95

(31)

where min denotes the minimum of the quantities within brackets for a given t. According to (31), the ˜ ˜ “drift” |(dβ(t)/dt)/ β(t)| of the evaporation conditions is the governing parameter which measures the influence of the weakly disturbed droplet kinetics on the concentration profile. For sessile or hanging droplet arrangements, as used for protein crystal growth, a sharp transition from the evaporation stage to the equilibrium ˜ ˜ stage is found (see, e.g. [2]) where |(dβ(t)/dt)/ β(t)| strongly increases. Then according to (31), the maximum observation period tobs where the derived formulae may be applied is limited from the start of the droplet evaporation to onset of this transitional stage. 4.3. Numerical estimations

Fig. 6. Microscopic images of evaporating droplets of aqueous protein solution hanging from a fine capillary (visible on top). White spots in the middle are reflections caused by illumination. (a) With a pattern from B´enard convection cells. Droplet of aqueous Des-ThrB30 human insulin solution under preparation conditions described in [19]. Temperature: 296.8 K, relative humidity 68%. Starting insulin concentration 4.0 mg/ml, starting droplet radius 0.040 cm. Recorded 3360 s after start of the evaporation, actual droplet radius 0.031 cm. Visualizing of the convection pattern by means of insulin precipitates in the supersaturated solution [19]. (b) Without B´enard convection cells. Droplet preparation: human serum albumine (HSA) dissolved in aqueous NaCl solution, starting concentrations 100 mg/ml HSA, 0.9% NaCl. Latex particles (0.08%) with 1.06 ␮m diameter (Molecular Probes Inc.) added to visualize convection patterns. Temperature 296.7 K, relative humidity 77%. Starting droplet radius 0.042 cm. Recorded 4620 s after start of the evaporation, actual droplet radius 0.031 cm. Homogenous intensity of scattered light across the droplet due to a homogeneous distribution of the latex particles in absence of convection cells.

The formulae (17a) and (31) combined with (3b), (4), (5) and (29b), respectively, represent main results of our paper. First we see from (3b) and (4) that, if the droplet kinetics approach the Langmuir kinetics, the characteristic parameter β¯ is given by

(g)  (g)  DA ρA ∼ ¯ |β| = (1 − f) (32) (l) (l) DP ρA For this limiting case, β¯ becomes independent of the starting droplet radius R0 (0) and can be adjusted only via the relative humidity f within the chamber. Us¯ ∼ ing the numerical values of Section 2, we find |β| = 1 1.41 × 10 (1 − f). Then the only possibility to avoid a significant solute concentration enhancement at the droplet boundary consists of choosing a large relative ¯ ≤ humidity f ≥ 0.93, which is required to get |β| 1 (see (17b) and Fig. 5, respectively). Such a scenario is typical in the case of protein crystal growth from the solution for systems with low concentrations of the precipitant electrolyte and a low protein solubility. Low solute concentrations are necessary, not only, to achieve a value of the solvent vapour pressure required for evaporation, but also to be sufficient to approach the Langmuir kinetics of the evaporating droplet.

96

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

There are some data in the literature allowing ¯ estimations of τ K and β(t) for stronger deviations from the Langmuir kinetics, as they occur, in practice, for the transitional stage in evaporating droplet arrangements before the approach to equilibration. We will consider a typical semi-spherical droplet with the starting volume V(0) ≈ 2.5 × 10−8 m3 , and τK ≈ 7.2 × 104 s (20 h). The Marangoni convection is assumed to disappear. We note that the evaporation kinetics are much slower by at least one order of magnitude compared with the above-mentioned arrangements where the Marangoni convection was observed. ¯ ≈ 2.4 with D(l) ≈ 3.2 × 10−11 m2 /s We find |β| P from above. Although the reported droplet volume change remains restricted to typical V(t)/V(0) = a(t)3 ≈ 1/2 so that the asymptotic concentration profile is not completely approached, we see that during evaporation, the concentration enhancement at the droplet boundary for the solute in question and, of course, for any lower diffusive solute is not negligible.

5. Conclusion The entrainment of a non-evaporating solute within an evaporating droplet is a simple effect. This effect exists even if there is no other superimposed molecular interaction at the droplet boundary, and which has to be taken into account, as the simplest approach, to analyze the local solute concentration, if nothing else is known about what happens at the moving droplet boundary. It has been shown here that a solute with a sufficiently low-diffusion coefficient will always exhibit a space and time-dependent concentration profile due to the solute entrainment, provided that there is no convective flux within the droplet. This concentration profile is determined by the interplay of the evaporation kinetics and the solute diffusion, where the former creates the solute concentration gradient, while the latter acts to equilibrate the enhanced concentration at the droplet boundary across the bulk of the droplet. It is shown here that the mathematical description is based on a generalized diffusion equation with boundary conditions at a moving boundary, which involves the effect of solute entrainment as an additional source term. This equation can be exactly solved assuming a

droplet evaporation kinetics of the Langmuir type. The singular perturbation theory is appropriate to generalize this exact solution for weak perturbations of the droplet kinetics as occurs in practice. It is useful to introduce two time constants to characterize the evaporation kinetics and the diffusion, namely the kinetic time constant τ K and the diffusion time τ D . A starting homogeneous solute concentration always develops into a monotonically increasing concentration profile with a maximum at the droplet boundary. In particular, this deformation is significant if τD /τK approaches or exceeds unity. If the droplet radius diminishes to about 0.7 R0 (0), the asymptotic profile is approached where the local concentration as related to the mean concentration becomes a function of the stretched radial coordinate x = r/R0 (t) only, and is governed solely by the parameter β¯ = −τD /τK . This asymptotic state becomes independent of the starting droplet radius for the limiting case of the Langmuir evaporation kinetics. It is shown further that the solute entrainment can act to enhance a starting tangential concentration variation if the evaporation kinetics are sufficiently fast. The effect of the concentration enhancement at the moving droplet boundary should not be underestimated. For a macromolecular solute under conditions occurring in practice (e.g., τD ≈ 105 s, τK ≈ 104 s), an enhancement by a factor of about 2 compared with the mean concentration level is possible, while at the midpoint of the droplet, the concentration scarcely exceeds its starting value. The results presented are of particular interest in various evaporating droplet techniques as used for the crystallization of biological macromolecules, because the nucleation and ripening will start near the droplet boundary. It should be noted that a similar concentration enhancement can also be met for arbitrarily small droplets if the evaporation kinetics are governed by the solvent vapour diffusion, because for this case, τ K is scaled with R0 (0)2 , in the same way as τ D is. Finally, the results presented here may be useful if a simple method is wanted to create a specific solute concentration gradient. The model assumptions as well as the derived formulae are valid only if any convective flux within the droplet can be excluded. Because, according to the present state of the art, the appearance of convections during droplet evaporation is difficult to predict, the convection-free status has to be checked by direct observation.

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

Acknowledgements This work was carried out within the frame of the JBCC project (Jena Biocrystallogenesis Centre) and supported by the Government of Thuringia through grant B378-000024 and by the Deutsche Forschungsgemeinschaft (SFB 604, TP C2). R.H. acknowledges support by The Fonds der Chemischen Industrie. The authors thank A. Walter for valuable discussions.

Appendix A. Mathematical appendix A.1. The diffusion equation with solute entrainment We consider the spherical shell R1 ≤ r ≤ R0 (t). df and df denote the surface element and the oriented surface element in the outer normal direction, respectively. In the following, D denotes the solute diffusion coefficient. The solute current according to the first Fick’s law is j (r, t) = −D grad c(r, t)

(A.1)

Taking into account the solute drag, the mass balance of the spherical shell reads dfj (r, t) + dfj (r, t) r=R0 (t)

d + dt



r=R1

r=R0 (t)

r=R1

d rc(r, t) = 0

(A.2)

where dfj (r, t) vanishes for r = R0 (t). The volume integral of (A.2) is evaluated according to d r=R0 (t) d rc(r, t) dt r=R1 r=R0 (t) ∂ dr c(r, t) = ∂t r=R1 dR0 (t) + dfc(r, t) dt r=R0 (t) r=R0 (t)  ∂ c(r, t) dr = ∂t r=R1  dR0 (t) + c(r, t)δ(r − R0 (t) + ε1 ) (A.3) dt where δ(r − R0 (t) + ε1 ) is Dirac’s delta function, and ε1 with ε1 → 0 is a small positive quantity. Con-

97

verting the surface integral on the spherical shell in (A.2) into a volume integral using the Gauss theorem, we find

r=R0 (t) r=R1

 d r div j (r, t) +

∂ c(r, t) ∂t

 dR0 (t) + c(r, t)δ(r − R0 (t) + ε1 ) = 0 dt

(A.4)

which must be satisfied for arbitrary R1 and R0 (t), so that the integrand of (A.4) vanishes. Taking into account (A.1), we arrive at the modified diffusion equation with a source term on the RHS:   1 ∂ 1 dR0 (t) ∆− c(r, t) = c(r, t)δ(r − R0 (t)) D ∂t D dt (A.5) which has to be solved with the boundary conditions c(r, t) continuous for r → 0, ∂ c(r, t) = 0 for r = R0 (t) + ε1 ∂r

(A.6)

and the initial condition given below. ∆ is the Laplacean operator. To solve (A.5) and (A.6) in spherical coordinates (r, ϑ, ϕ), we introduce the usual ansatz c(r, t) = c(r, ϑ, ϕ, t) ∞   = cν(µ) (r, t)Pνµ (cos ϑ)eiµϕ

(A.7)

ν=0 −ν≤µ≤ν

We will carry out a coordinate transformation from (r, t) to (r1 , t1 ) to replace the second boundary condition (A.6) by a new one at a fixed boundary: cν(µ) (r, t) = cˆ ν(µ) (r1 , t1 ) r1 = p(r, t) =

t

t1 = q(t) = 0

a¯ (t) =

R0 (t) R0 (0)

r a¯ (t) dt 

1 a¯ (t  )2

(A.8a) (A.8b)

(A.8c)

(A.8d)

98

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

ˆ ν by Defining the differential operator D   2 ˆ ν ≡ d + 2 + βr1 d − ν(ν + 1) D r1 dr1 dr12 r12

Now we have

  2 ∂2 ∂ + + β(t )r 1 1 r1 ∂r1 ∂r12  ν(ν + 1) 1 ∂ − − cˆ (µ) (r1 , t1 ) D ∂t1 ν r12 = R10 β(t1 )ˆcν(µ) (r1 , t1 )δ(r1 − R10 ) cˆ ν(µ) (r1 , t1 ) continuous for r1 → 0, ∂ (µ) cˆ (r1 , t1 ) = 0 for r1 = R10 + ε1 ∂r1 ν

(µ)

the Laplace transform c˜ ν (r1 , s) obeys the differential equation and boundary condition (A.9a)

lim cˆ ν(µ) (r1 , t1 ) = cˆ ν(µ) (r1 , 0) = cν(µ) (r, 0)

(A.9c)

(A.10a)

the function β(t1 ) entering (A.9a) is found from β(t1 ) =

1 ˜ −1 β(q (t1 )) D



 1 cˆ (µ) (r1 , 0) D ν (A.14a)

(A.9b)

˜ via observing that r = r1 at t = 0. Defining β(t) d ˜ a¯ (t)2 = 2β(t) dt

ˆ ν − s¯ )˜cν(µ) (r1 , s) (D = R10 β˜cν(µ) (r1 , s¯ )δ(r1 − R10 ) −

where R10 = R0 (0). Eqs. (A.9a) and (A.9b) have to be solved with the initial condition t1 →0

(A.13)

(A.10b)

where q−1 (t1 ) = t is the inverse function of (A.8c). A.2. The solution for the Langmuir kinetics

c˜ ν(µ) (r1 , s) continuous, for r1 → 0,   d c˜ (µ) (r1 , s) = 0, for r1 = R10 + ε1 dr1 ν (A.14b) s (A.14c) s¯ = D The formal solution of (A.14a) and (A.14b) is R10 dr2 Gν (r1 , r2 , s¯ )f(r2 , s¯ ) c˜ ν(µ) (r1 , s) = (A.15) 0

where f(r1 , s¯ ), 0 ≤ r1 ≤ R10 + ε1 , abbreviates the RHS of (A.14a), and Gν (r1 , r2 , s¯ ) is Green’s function ˆ ν − s¯ which satisfies of D ˆ ν − s¯ )Gν (r1 , r2 , s¯ ) = δ(r1 − r2 ) (D (A.16)

For the special case of the Langmuir kinetics we have

and the boundary condition (A.14b). Inserting the RHS (µ) of (A.14a) into (A.15), c˜ ν (R10 , s) has to be determined first solving a linear algebraic equation. The final result is

˜ = − 1 = const, β(t) τK 1 β(t1 ) = − = const DτK  ˜ a¯ (t) = 1 + 2βt

c˜ ν(µ) (r1 , s)  R10 βR10 Gν (r1 , R10 , s¯ )Gν (R10 , r2 , s¯ ) =− dr2 1 − βR10 Gν (R10 , R10 , s¯ ) 0  (µ) cˆ ν (r2 , 0) + Gν (r1 , r2 , s¯ ) (A.17) D

t1 =

1 ˜ ln(1 + 2βt) 2β˜

(A.11a) (A.11b) (A.11c)

We will change to the Laplace transform (see, e.g. [27]) ∞ (µ) c˜ ν (r1 , s) = dt1 e−st1 cˆ ν(µ) (r1 , t1 ) (A.12) 0

ˆ ν can be rewritten to become self-adjoint Because D (see, e.g. [28,29]), Gν (r1 , r2 , s¯ ) is represented as (see, e.g. [28])

 ∞  βr22 ψνi (r1 )ψνi (r2 ) 2 Gν (r1 , r2 , s¯ ) = − r2 exp , λνi + s¯ 2 i=0

if 0 ≤ r1 ≤ R10 + ε1 , 0 ≤ r2 ≤ R10 + ε1

(A.18)

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

ψνi (r1 ), i = 0, 1, 2, . . . denote the eigenfunctions of ˆ ν , which satisfy D ˆ ν ψνi (r1 ) = −λνi ψνi (r1 ) D

(A.19a)

and the boundary condition (A.14b). The eigenvalues λνn , n = 0, 1, 2, . . . are real and discrete, where λν0 < λν1 < λν2 < · · · → ∞, and λ00 = 0, otherwise λν0 > 0 if ν > 0 and β ≤ 0. The eigenfunctions are orthonormalized according to (see, e.g. [28])

 R10 βr12 2 dr1 r1 exp ψνi (r1 )ψνj (r1 ) = δij (A.19b) 2 0 with δij being the Kronecker symbol. Instead of (A.18), we also use the simple expression (see, e.g. [28]) valid for r2 = R10 Gν (r1 , R10 , s¯ ) = −

Ψν (r1 , s¯ , β) , Ψν (R10 , s¯ , β)

if r1 < R10 (A.20a)

Gν (r1 , R10 , s¯ ) = 0,

if r1 > R10

(A.20b)

ˆν− where Ψν (r1 , s¯ , β) is that particular solution of (D s¯ )Ψν (r1 , s¯ , β) = 0 which is continuous for r1 → 0 (see, e.g. [20]):   r1 ¯ , σ, β (A.21a) Ψν (r1 , s¯ , β) = Φν R10 d 1 = Ψν (r1 , s¯ , β) = dr1 R10    σ − νβ¯ 1+ν σ β¯ ν 5 × x 1F 1 1 − + , + ν, − x2 3 + 2ν 2 2 2 2β¯   ¯ σ ν 3 β + νxν 1 F 1 − + , + ν, − x2 (A.21b) ¯ 2 2 2 2β

Ψν (r1 , s¯ , β)

1F 1

is the confluent hypergeometric function. We introduced the dimensionless Laplace variable σ = s¯ R210 = s

R210 D

(A.21c)

as well as the dimensionless radial coordinate x = r1 /R10 and the dimensionless parameter β¯ = βR210 , in accordance with (5) and (8). The eigenfunctions ψνi (r1 ) and theeigenvalues λνi are

99

given by ψνi (r1 ) = Cνi Ψν (r1 , −λνi , β), Ψν (R10 , −λνi , β) = 0

(A.22)

where Cνi has to be determined by (A.19b). The relations (12a) to (13b) are found if Ψν (r1 , s¯ , β) is determined in the WKB approximation (see, e.g. [29]). For ν = 0, (A.17) becomes simplified if, according to (6) and (A.9c), a constant starting con(0) (0) centration c0 (r, 0) = cˆ 0 (r1 , 0) = c0 = const is ˆ 0 − s¯ )ψ = assumed. Observing that the equation (D 1 together with (A.14b) has the trivial solution R ψ = −1/¯s = 0 01 dr2 G0 (r1 , r2 , s¯ ), we find, applying the second part of this relation only and using (A.20a) to (A.21b): c0 (0) c˜ 0 (r1 , s) = D¯s   Ψ0 (r1 , s¯ , β) × 1 − βR10  Ψ0 (R10 , s¯ , β) + βR10 Ψ0 (R10 , s¯ , β)

 ¯ R210 c0 Φ (x, σ, β) 0 = 1 − β¯  ¯ + βΦ ¯ 0 (1, σ, β) ¯ D σ Φ0 (1, σ, β) (A.23a) ¯ = d Φν (x, σ, β) ¯ Φν (x, σ, β) dx

(A.23b)

For further evaluation we note that Gν (r1 , r2 , s¯ ) taken as analytic function of s¯ is represented, according to (A.18), as a partial fraction series with first-order poles only at s¯ = −λνn , n = 0, 1, 2, . . . . The residues of Gν (R10 , R10 , s¯ ) are strictly negative. (µ) Therefore the only singularities of c˜ ν (r1 , s) given by (A.17) are first-order poles from Gν (r1 , r2 , s¯ ) at s¯ = −λνn , n = 0, 1, 2, . . . , and first-order poles from the zeros of the denominator at s¯ = s¯νn , where s¯νn > −λνn > s¯ν,n+1 , n = 0, 1, 2, . . . . (0) For c˜ 0 (r1 , s) given by (A.22a), this reduces to first-order poles at s¯ = −λ00 = 0 and s¯ = s¯0n . Taking into account (A.17), (A.20a), (A.23a), (A.23b), s¯νn = s¯νn (β) is found from s¯νn (β) = 2 , where the σ (β) ¯ σνn (β)/R νn ¯ are the roots of the 10 equation ¯ + βΦ ¯ ν (1, σ, β) ¯ =0 Φν (1, σ, β)

(A.24)

100

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

Thus we arrive at ¯ Φ0 (x, σ, β)  ¯ + βΦ ¯ 0 (1, σ, β) ¯ Φ0 (1, σ, β) ∞  ρn (β) ¯ ¯ β) ¯ = Φ (x, σ0n (β), ¯ 0 σ − σ0n (β) n=0

(A.25a)



 ∞   1 1 1 + β¯ − ¯ σ σ σ0n (β) n=0  ¯ ρn (β) ¯ β) ¯ × Φ (x, σ0n (β), (A.25b) ¯ 0 σ0n (β)

(0) c˜ 0 (r1 , s)

R2 = c0 10 D

¯ denotes the residuum of 1/(Φ (1, σ, β)+ ¯ where ρn (β) 0 ¯ 0 (1, σ, β)) ¯ at σ = σ0n (β). ¯ σνn (β) ¯ and ρn (β) ¯ have βΦ to be calculated numerically (see Appendix B) except of (10a) which follows from the mass conservation. Transforming from (A.25b) to the original function (see, e.g. [27]) of original coordinates, we get (9b), ¯ is omitted, for simplicwhere the first index of σ0n (β) ity.

where τ is treated as parameter in each perturba˜ denotes the known (and tion order. If h = h(t, β) exact) result for β˜ = const, the zeroth-order according to (A.26b) is defined straightforwardly via ˜ β˜ = β(τ). We will restrict ourselves to deriving sufficient conditions to neglect the first perturbation order compared with the zeroth one. Starting with (A.10a) and (A.8c) (converted into a differential equation), we find that (A.11b) and (A.11c) can be used inserting ˜ if β˜ = β(t)   ˜  t  dβ(t) ˜ (A.27) for 0 < t < tobs    |β(t)|, 2  dt  To solve (A.9a) to (A.9c), we introduce the slow transformed time variable τ1 = εt1 writing β = β(τ1 ), so ˆ ν given by (A.13) now becomes D ˆ ν (τ1 ). We that D start with the ansatzes (µ) cˆ ν(µ) (r1 , t1 ) = c¯ˆ ν (r1 , t1 , τ1 ) (A.28a) ∞ (µ) (µ) dt1 e−st1 c¯ˆ ν (r1 , t1 , τ1 ) (A.28b) c¯˜ ν (r1 , s, τ1 ) = 0

(µ)

A.3. Deviations from the Langmuir kinetics In the presence of weak disturbances of the Langmuir kinetics, β˜ given by (A.10a) becomes a ˜ time-dependent function β˜ = β(t). It is obvious that ˜ the formulae derived for β = const remain approx˜ imately valid if β(t) depends upon t “being sufficiently slow”. Such considerations suggest to study ˜ the general case β˜ = β(t) applying the methods of the singular perturbation theory (see, e.g. [30]). We introduce the “slow” time variable τ = εt, where ε = 1 is a formal perturbation parameter and τ is ˜ ˜ defined by writing now β(τ) instead of β(t). Any other time-dependent function h = h(t) that depends upon ¯ τ) so that β˜ will be rewritten to become h = h(t, ¯ τ). dh(t)/dt has to replaced by ((∂/∂t) + ε(∂/∂τ))h(t, Starting with a differential equation to determine h = ¯ τ) can be successively calculated from h(t), thus h(t, the perturbation ansatz ¯ τ) = h¯ 0 (t, τ) + εh¯ 1 (t, τ) + ε2 h¯ 2 (t, τ) + · · · h(t, (A.26a) ˜ h¯ 0 (t, τ) = h(t, β(τ))

(A.26b)

The Laplace transform c¯˜ ν (r1 , s, τ1 ) obeys, in generalizing (A.14a):   ε ∂ (µ) ˆ − s¯ c¯˜ ν (r1 , s¯ , τ1 ) D ν (τ1 ) − D ∂τ1 (µ)

= R10 β(τ1 )c¯˜ ν (r1 , s¯ , τ1 )δ(r1 − R10 ) 1 (µ) − cˆ¯ ν (r1 , 0, 0) D

(A.29)

ˆ ν (τ1 ) is given where (A.14b) and (A.14c) holds, and D by (A.13) with β = β(τ1 ). If the RHS of (A.29) is abbreviated, for a moment, by f(r1 , s¯ , τ1 ), the pertur(µ) bation series of c¯˜ ν (r1 , s, τ1 ) according to the scheme of the singular perturbation theory reads (µ) c¯˜ ν (r1 , s, τ1 ) R10 ε dr2 Gν (r1 , r2 , s¯ , τ1 )f(r2 , s¯ , τ1 ) + = D 0 R10 × dr2 Gν (r1 , r2 , s¯ , τ1 ) 0

∂ × ∂τ1



R10 0

dr3 Gν (r2 , r3 , s¯ , τ1 )f(r3 , s¯ , τ1 ) + · · · (A.30)

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

Gν (r1 , r2 , s¯ , τ1 ) is given by (A.18), (A.20a) and (A.20b), respectively, where β = β(τ1 ). Insert(µ) ing the RHS of (A.29) into (A.30), c¯˜ ν (r1 , s, τ1 ) is found as before. Restricting to ν = 0 and (0) c¯ˆ 0 (r2 , 0, 0) = c0 = const, and taking into account  R01 0 dr2 G0 (r1 , r2 , s¯ , τ1 ) = −1/¯s, we arrive at the first-order result (0) c˜¯ 0 (r1 , s, τ1 )   c0 β(τ1 )R10 L0 (r1 , R10 , s¯ , τ1 ) = 1+ s¯ 1 − β(τ1 )R10 L0 (R10 , R10 , s¯ , τ1 ) (A.31a)

L0 (r1 , r2 , s¯ , τ1 )

ε = G0 (r1 , r2 , s¯ , τ1 ) + D R01 × dr3 G0 (r1 , r3 , s¯ , τ1 ) 0   dβ(τ1 )/dτ1 ∂ × + G0 (r3 , r2 , s¯ , τ1 ) β(τ1 ) ∂τ1 (A.31b)

∂ 1 dβ ∂ G0 (r1 , r2 , s¯ , τ1 ) = r1 G0 (r1 , r2 , s¯ , τ1 ) ∂τ1 s¯ dτ1 ∂r1 (A.31c) where (A.31c) is found from differentiation of ˆ 0 (τ1 ) − s¯ )G0 (r1 , r2 , s¯ , τ1 ) − δ(r1 − r2 ) ≡ 0 with (D respect to τ 1 . (0) c¯˜ 0 (r1 , s, τ1 ) given by (A.31a) is the generaliza(0) tion of c˜ 0 (r1 , s) given by (A.23a), and exhibits the same analytical properties, taken as function of s¯ . In particular, there are first-order poles only at s¯ = 0, ¯ 1 )) = and s¯ = s¯0n (β(τ1 )) + δ¯s0n , with σ0n = σ0n (β(τ s¯0n (β(τ1 ))R210 as root of (A.24). The shift δ¯s0n of the pole positions and the corrected residues of (A.31a) compared with (A.23a) indicate the change of the (µ) time and space-dependence of c¯ˆ ν (r1 , t1 , τ1 ), respectively, as they appear in the first perturbation order. We will take δ¯s0n , n = 0, 1 as measure of the strength of all first-order corrections because (A.31a) is governed by the only function L0 (r1 , R10 , s¯ , τ1 ). First (µ) analyzing the asymptotic behaviour of c¯ˆ ν (r1 , t1 , τ1 ) for large t1 , the condition to neglect the first

101

perturbation order compared with the zeroth one reads |δ¯s00 |  |¯s00 (β(τ1 ))| = 3|β(τ1 )|

(A.32a)

(µ) Quite similar, the transitional stage of c¯ˆ ν (r1 , t1 , τ1 ) is not affected by first-order perturbation corrections if

|δ¯s01 |  |¯s01 (β(τ1 ))|

(A.32b)

δs0n , n = 0, 1 is found from ∂ G0 (R10 , R10 , s¯0n , τ1 ) ∂¯s ε R01 ∼ dr3 G0 (R10 , r3 , s¯0n , τ1 ) = − D 0   dβ(τ1 )/dτ1 ∂ G0 (r3 , R10 , s¯0n , τ1 ) × + β(τ1 ) ∂τ1 (A.33)

δ¯s0n

The numerical evaluation of (A.33) using (A.18) and (A.31c) shows that the second term on the RHS of (A.33) is small compared with the first one, so that (A.33) strongly simplifies to 1 dβ/dτ1 δ¯s00 ∼ = δ¯s01 ∼ = D β

(A.34)

Inserting (A.34) into (A.32a) and (A.32b) and carrying out the back transformation to t, we get (31), where (A.27) is taken into account, too.

Appendix B. Analytic approximants ¯ = σn (β). ¯ Approximations Abbreviation: σ0n (β) valid for −10 ≤ β¯ ≤ −0.25. ¯ = − 20.1996 − 2.5156β¯ − 0.1004β¯ 2 − 0.0010β¯ 3 σ1 (β) ¯ = − 59.6782 − 2.4976β¯ − 0.0857β¯ 2 + 0.0002β¯ 3 σ2 (β) ¯ = − 118.898 − 2.4973β¯ − 0.0839β¯ 2 + 0.0002β¯ 3 σ3 (β) ¯ = − 197.857 − 2.4981β¯ − 0.0836β¯ 2 + 0.0001β¯ 3 σ4 (β) ¯ = ρ0 (β)

8.7702 152.42 413.87 345.14 + + + 2 3 ¯ ¯β − 3 ¯ (β − 3) (β − 3) (β¯ − 3)4

¯ =− ρ1 (β)

27.096 552.20 1828.2 1997.9 − − − β¯ − 3 (β¯ − 3)2 (β¯ − 3)3 (β¯ − 3)4

102

¯ ρ2 (β)=

Th. Klupsch et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 231 (2003) 85–102

44.246 1003.4 3605.2 4270.1 + + + 2 3 ¯β − 3 ¯ ¯ (β − 3) (β − 3) (β¯ − 3)4

¯ ρ3 (β)= − ¯ ρ4 (β)=

59.537 1407.8 5133.1 6155.8 − − − 2 3 β¯ − 3 (β¯ − 3) (β¯ − 3) (β¯ − 3)4

75.050 1808.3 6629.9 7986.5 + + + β¯ − 3 (β¯ − 3)2 (β¯ − 3)3 (β¯ − 3)4

References [1] P. Chen, O.I. del Rio, A.W. Neumann, in: A. Baszkin, W. Norde (Eds.), Physical Chemistry of Biological Interfaces, Marcel Dekker, New York, 2000, p. 523. [2] A. Ducruix, R. Giege, in: A. Ducruix, R. Giege (Eds.), Crystallization of Nucleic Acids and Proteins, A Practical Approach, second ed., Oxford University Press, Oxford, 1999, p. 130. [3] P. Mühlig, Th. Klupsch, A. Walter, R. Hilgenfeld, in: Proceedings of the Contribution P-B.9 to the ICCBM9 Conference, Jena, 23–28 March 2002 (book of abstracts), and to be published. [4] A.M. Schwartz, K.A. Berglund, J. Cryst. Growth 210 (2000) 753. [5] P. Kékicheff, R.G. Laughlin, R.L. Munyon, Langmuir 17 (2001) 4696. [6] W.W. Fowlis, L.J. DeLucas, P.J. Twigg, S.B. Howard, E.J. Meehan, J.K. Baird, J. Cryst. Growth 90 (1988) 117. [7] G.T. DeTitta, J.R. Luft, Acta Cryst. D51 (1995) 786. [8] J.R. Luft, D.T. Albright, J.K. Baird, G.T. DeTitta, Acta Cryst. D52 (1996) 1098. [9] J.K. Baird, J. Cryst. Growth 204 (1999) 553. [10] A.A. Chernov, in: B.K. Vainshtein, A.A. Chernov, L.A. Shuvalov (Eds.), Modern Crystallography. III. Crystal Growth, Springer-Verlag, Berlin, 1984 (Chapters 4 and 6). [11] J. Drenth, C. Haas, Acta Cryst. D54 (1998) 867.

[12] O. Galkin, P.G. Vekilov, J. Phys. Chem. B 103 (1999) 10965. [13] L.D. Landau, E.M. Lifshitz, Lehrbuch der Theoretischen Physik, Band X: Physikalische Kinetik, Akademie Verlag, Berlin, 1983, Paragraphs 99 and 100, p. 457. [14] S. Cudney, K. Patel, K. Weisgraber, Y. Newhouse, A. McPherson, Acta Cryst. D50 (1994) 414. [15] I. Langmuir, Phys. Rev. 12 (1918) 368. [16] N.A. Fuchs, Evaporation and Droplet Growth in Gaseous Media, Pergamon Press, Oxford, 1959. [17] H.Y. Erbil, M. Dogan, Langmuir 16 (2000) 9267. [18] K.S. Schmitz, An Introduction to Dynamic Light Scattering by Macromolecules, Academic Press, Boston, 1990 (Chapters 2 and 3). [19] P. Mühlig, Th. Klupsch, U. Schell, R. Hilgenfeld, J. Cryst. Growth 232 (2001) 93. [20] H. Bateman, A. Erdelyi, Higher Transcendental Functions, vol. I, McGraw-Hill, New York, 1953 (Chapter 6). [21] A.T. Chai, N. Rashidnia, V.S. Arpaci, in: Proceedings of the VIIIth European Symposium on Materials and Fluid Sciences in Microgravity, SP 333, ESA, Brussels, Belgium, August 1992, p. 187. [22] V.X. Nguyen, K.J. Stebe, Phys. Rev. Lett. 88 (2002) 164501. [23] R.D. Deegan, O. Bakajin, T.F. Dupont, G. Huber, S.R. Nagel, T.A. Witten, Nature 389 (1997) 827. [24] R.D. Deegan, O. Bakajin, T.F. Dupont, G. Huber, S.R. Nagel, T.A. Witten, Phys. Rev. E 62 (2000) 756. [25] P. Mühlig, Th. Klupsch, R. Hilgenfeld, in: Proceedings of the Contribution P-B.10 to the ICCBM9 Conference, Jena, 23–28 March 2002 (book of abstracts), and to be published. [26] J.R. Luft, G.T. DeTitta, Acta Cryst. D51 (1995) 780. [27] G. Doetsch, Anleitung zum praktischen Gebrauch der Laplace-Transformation, Oldenbourg, München, 1961. [28] W.I. Smirnov, Lehrgang der höheren Mathematik, Teil IV, VEB Deutscher Verlag der Wissenschaften, Berlin, 1958 (Chapter IV). [29] M.V. Fedoryuk, Asymptotic Analysis, Springer-Verlag, Berlin, 1993. [30] A.H. Nayfeh, Introduction to Perturbation Techniques, Wiley, New York, 1981 (Chapter 10).