Article pubs.acs.org/ac

Frequency Shifts of a Quartz Crystal Microbalance Calculated with the Frequency-Domain Lattice−Boltzmann Method: Application to Coupled Liquid Mass Diethelm Johannsmann*,† and Gunther Brenner‡ †

Institute of Physical Chemistry, Clausthal University of Technology, 38678 Clausthal-Zellerfeld, Germany Institute of Applied Mechanics, Clausthal University of Technology, 38678 Clausthal-Zellerfeld, Germany

‡

S Supporting Information *

ABSTRACT: In recent years the quartz crystal microbalance (QCM) has seen an impressive evolution from a film-thickness monitor to a surface-analytical instrument with capabilities much beyond gravimetry. In particular, the instrument has often been applied to adsorbates from a liquid phase and, also, to samples with structure in the surface plane. In order to quantitatively predict frequency shifts induced by such samples from a model, one needs to compute the in-phase component of the area-averaged periodic tangential stress at the resonator surface. A method is described which performs this task, making use of a variant of the Lattice−Boltzmann (LB) method. The algorithm differs from the conventional LB method in that it deals with oscillatory flows and only covers linear hydrodynamics. The adsorption of small particles (mimicking proteins) was chosen as an example to test the performance. These samples are acoustically thin, which simplifies the calculations. Also, the material’s finite compliance can be neglected in this limit. The simulations predict the amount of solvent trapped between neighboring particles, which contributes to the adsorbate’s apparent mass. The unknown amount of hydrodynamically coupled liquid is a serious problem in the interpretation of QCM experiments. On an experimental level, the amount of trapped solvent can be estimated from the comparison of the optical layer thickness (determined with ellipsometry) and the acoustic layer thickness (determined with a QCM). Since the amount of trapped liquid decreases when neighboring particles aggregate into clusters, this analysis can lead to a statement on the degree of clustering. The LB-based simulations show, though, that the relation between the cluster geometry and the amount of trapped solvent is highly nontrivial. The details of the geometry do matter. The LB-based algorithm can calculate the amount of trapped solvent for userspecified particle shapes, orientations, interparticle distances, and also distributions thereof. It is an essential step in the quantitative interpretation of QCM results obtained on thin samples with in-plane structure.

T

he quartz crystal microbalance (QCM) is an established tool to measure film thickness in the nanometer range.1 For sufficiently rigid layers, the frequency shift is negative and proportional to the mass per unit area,2 which motivates the term “microbalance”. Importantly, the QCM can reveal information about the layer of interest much beyond the layer’s mass. When a QCM is brought into contact with a soft sample, there is a shift in the resonance bandwidth, 2Γ, in addition to a shift in resonance frequency, f. Also, one may determine Δf and ΔΓ on a number of different overtones. While this set of experimentally accessible parameters opens up numerous novel ways of studying soft interfaces,3,4 it is not always clear how to physically interpret the values of Δf and ΔΓ. Explicit equations exist which predict Δf and ΔΓ for planar layer systems,5 but in order to analytically describe samples with in-plane structure, simplifying assumptions are needed, which only are valid in special situations.6,7 There are numerous interesting samples with in-plane structure, sketched in Figure 1. These include sparse protein layers,8 adsorbed nanoparticles,9,10 bacteria,11 cells,12 vesicles,13 and nanobubbles.14 All of these are highly relevant to the field of soft matter research, and all of these require numerical modeling. © XXXX American Chemical Society

Figure 1. Different types of structured samples encountered in QCM experiments. Sketched are a protein, a micelle, a liposome, a nanoparticle, a biological cell, and a bacterium (not to scale).

Received: May 22, 2015 Accepted: June 29, 2015

A

DOI: 10.1021/acs.analchem.5b01912 Anal. Chem. XXXX, XXX, XXX−XXX

Article

Analytical Chemistry In principle, the prediction of Δf and ΔΓ from a sample’s geometry and materials parameters is straightforward. Δf and ΔΓ are proportional to the in-phase and the out-of-phase components of the area-averaged tangential stress at the surface of the resonator plate. This is the content of the small load approximation.15,16 More quantitatively, one has the relation Δf + iΔΓ i i ⟨σS⟩area = ZL = f0 πZq πZq uS

Lattice−Boltzmann Method” (FD-LBM). The algebra goes back to Shi and Sader.23 We go through a brief explanation below. The modification comprises two steps. First, nonlinear terms are omitted throughout. “Nonlinear terms” here are all terms which are of second or higher order in velocity. With nonlinear terms omitted, the calculation evidently is limited to flows at low Reynolds numbers. This is not a problem because QCM experiments indeed occur in this regime. There is a second consequence of linearization, which is that Euler coordinates and Lagrange coordinates are equivalent. The issue of coordinates moving with the medium or resting in space is irrelevant in linear acoustics (that is, in the smallamplitude limit). One may therefore leave all boundaries fixed in space. Even though the boundaries themselves undergo an oscillatory movement, the field of oscillation amplitudes does not. FD-LBM entails a further modification of standard LBM, in that it computes complex amplitudes of a periodic flow field, rather than the flow field itself. One writes the vector field of velocities, u(r,t), as

(1)

Δf + iΔΓ is a complex frequency shift, which is the prime motivation to quantify dissipative processes with the parameter Γ. Γ is the half-bandwidth at half-maximum. It is related to the inverse Q-factor by the relation Γ = Q−1f/2. Q−1 is sometimes also called the “dissipation factor”, D.17 f 0 is the resonance frequency of the fundamental (5 MHz in the simulations reported below), and Zq = 8.8 × 106 kg/(m2 s) is the shearwave impedance of the crystal. σS and uS are the tangential stress and the tangential velocity at the resonator surface, respectively. Angle brackets (⟨...⟩) denote an area-average. The ratio of stress and velocity is the load impedance, ZL. The smallload-approximation entails a few difficulties,18 but it is still widely used in the analysis of QCM-experiments. This leaves the problem of calculating the area-averaged periodic stress for structured adsorbates. In the past, the Finite Element Method (FEM) has been used to this end.7 The formalism (employing COMSOL’s Multiphysics module) is successful, in principle, but only works in two dimensions. FEM software solving similar problems in 3D exists,19 but all attempts to apply such codes to QCM-related problems have failed so far. Calculations of Δf and ΔΓ carried out on 2D geometries are not necessarily meaningless. There are a number of general problems in the physics of the QCM, which can be addressed in the frame of 2D simulations. An example is the coupled resonance.20 But, of course, such calculations are not quantitatively predictive for the geometries encountered in realworld experiments. In the work described below, we propose to model QCM experiments with a modification of the Lattice−Boltzmann method (LBM).21 The Lattice−Boltzmann method implements a computationally convenient discretization of the Boltzmann equation, using a lattice in space and time and assigning distribution functions to the velocities at each site. The number of vectors sampling velocity space varies, but it is always much smaller than the number of lattice sites sampling real space. The number is 19 in the present work, which employs the D3Q19 scheme (following the nomenclature proposed in ref 22, “D3” denotes a calculation in three dimensions). The program basically consists of alternating steps of “streaming” and “relaxation”. In the streaming step, each subpopulation with a certain velocity is shifted to the neighboring lattice site corresponding to the respective velocity and the time step. In the relaxation step, the populations on the individual sites are redistributed, such that they relax toward local equilibrium. After incorporation of the boundary conditions, the algorithm produces solutions of the Navier−Stokes equation, valid in the limit of small Mach numbers. (The Mach number is the ratio of velocity and speed of sound.) Simple LBM code can run on personal computers. LBM code can be parallelized rather easily because the computationally most expensive step (the relaxation) occurs locally on each lattice site. We make use of a modified version of the Lattice− Boltzmann method, which we call the “Frequency-Domain

u(r , t ) = Re(uE(r , t ) exp(iωt ))

(2)

ω is the frequency of excitation. The source of excitation is turned on at t = 0 and is strictly time-harmonic from thereon. In the problems treated below, excitation occurs from the bottom surface, which is a wall undergoing tangential oscillation. uE(r,t) is an Envelope. The relaxation step in FDLBM is formulated such that the code computes uE rather than u. In general, uE depends on space and time, but the time dependence disappears once the simulation reaches the stationary state. Since eq 2 is complex, uE must also be complex. All steps of computation are carried out on both the in-phase and the out-of-phase components of the flow. The added computational effort roughly compensates what was gained by omission of the nonlinear terms in the relaxation step. In terms of computational efficiency, FD-LBM and standard LBM are comparable. Turning the velocity, u(r,t), into a product as in eq 2 is advantageous because the function uE(r,t) becomes constant in time after ringing-in is completed; the solution evolves toward constant uE. uE then turns into a vector field of complex oscillation amplitudes, uE,∞(r). This stationary state (with regard to the envelope) is characterized by an entropy production rate attaining its minimum value. That such a state is always reached is proven in nonequilibrium thermodynamics.24 The time scale of ringing-in is given by hC2/ν, with hC being the thickness of the simulation volume and ν the kinematic viscosity. In FD-LBM, streaming and relaxation are iterated until the steady state is reached. (The bounce-back step accounting for the boundary conditions may in this context be viewed as part of the streaming step.) All further analysis is based on the properties of the steady state, alone; it does not rely on the time-dependent processes leading to it. One may therefore start the simulation from unphysical initial conditions. For instance, the decaying shear wave (which solves the problem in the absence of particles or adsorbates) usually is closer to the shape of the final flow pattern than the quiescent fluid. Starting the simulation from the decaying shear wave shortens the time needed to reach stationarity. At this point, we slightly digress and mention a further key advantage of the FD-LBM scheme, which is played out when the liquid contains particles with dynamical behavior of their own. Such particles might be floating in the bulk of the liquid (in which case they undergo translation and rotation) or be B

DOI: 10.1021/acs.analchem.5b01912 Anal. Chem. XXXX, XXX, XXX−XXX

Article

Analytical Chemistry

Sauerbrey,2 thin, rigid layers decrease the resonance frequency, following the equation

adsorbed to the surface (in which case they may be deformed by the shear stress acting onto that surface). FD-LBM simulations covering particle dynamics will be described in a later paper. Generally speaking, the simulation of immersed particles with the LB method is plagued by numerical instabilities.25 These are immediately evident in FD-LBM, as well. If one attempts to regularly update a particle’s amplitude of motion such that mechanical equilibrium holds at the particle surface, the simulation quickly diverges. (“Mechanical equilibrium” here includes the inertial forces linked to the term exp(iωt).) Importantly, all such instabilities can be damped out by updating the particle velocity gradually. If the algorithm finds that a particle’s oscillation amplitude is different from what would be required by mechanical equilibrium, the algorithm shifts the amplitude toward equilibrium in small steps, rather than immediately imposing equilibrium. It thereby slows down the approach to mechanical equilibrium and avoids the instability. Since the amplitude is stationary in the final state, the delayed updating does not affect the properties of the final state. Even though the FD-LBM simulation no longer faithfully reproduces the ringing-in process, it still faithfully reproduces the stationary state (which is all that matters). In the work described below, we apply the Shi−Sader formalism to the geometry of the QCM, but we avoid the consequences of nontrivial particle dynamics by only considering small, rigid particles tightly adsorbed to the resonator surface. For such samples, the velocity at the particle surface is strictly equal to the velocity of the substrate; there is no need to separately calculate it. The geometry is depicted in Figure 2. Excitation occurs from the lower boundary of the cell,

Δf =

if0 πZq

iωΔm = −

2nf02 Zq

Δm = − f

Δm mq

(3)

n is the overtone order and Δm is the mass per unit area. The Sauerbrey equation can be quickly deduced from the small load approximation eq 1 by inserting the term iωΔm uS for the stress, ⟨σS⟩. The modal mass per unit area, mq, was introduced in the last step on the right-hand side. For the ideal thicknessshear resonator, the modal mass is mq = Zq/(2f 0). (The modal mass slightly differs from mq for real resonators,26 but this problem is outside the scope here.) Importantly, the Sauerbrey relation per se does not assume planar layer systems; Δm is an area-averaged mass per area. The sample must be rigid and much thinner than the wavelength of sound. Otherwise, there are no constraints on geometry. A well-known problem with the Sauerbrey equation when applied to adsorbates in the liquid phase is the unknown amount of trapped liquid. If the adsorbate has in-plane structure, the liquid contained in the gaps moves together with the adsorbate, thereby contributing to the adsorbate’s effective mass (also called “Sauerbrey mass”). The fact that the Sauerbrey mass is larger than the true mass makes quantitative interpretation of adsorption experiments difficult. Clearly, knowledge of the effective mass can amount to a valuable piece of information on the sample, if the true adsorbed mass is known independently. An experimental assessment of the true mass can be based on optical techniques such as ellipsometry or surface plasmon resonance (SPR) spectroscopy. Combined optical and acoustic measurements on thin adsorbates have a long history.27−29 It turns out that the adsorbed mass as determined optically is closer to the true mass than the adsorbed mass determined acoustically. The reasons are discussed in refs 26 and 30. If one accepts the “optical mass” as representative of the true mass, one can infer the amount of trapped water from the comparison of optical and acoustic mass. In the introduction, it was emphasized that the QCM carries information beyond gravimetry because it reports the bandwidth in addition to the resonance frequency and because it does so on many overtones. With thin rigid layers, this added information comes into play, indirectly: Whether or not an experiment may be analyzed using the Sauerbrey equation can be checked. If the Sauerbrey equation holds, the change in bandwidth is much smaller than the change in frequency and, also, the normalized frequency shifts, Δf/n, are the same on all overtones. This check is frequently undertaken in experiment and will be undertaken on the simulation results in the Test of the Sauerbrey Model section. We will find these two predictions to be largely confirmed. This is not a trivial result because the liquid contained in the pockets between particles might affect the bandwidth and the overtone dependence of Δf/n. (It does, according to both the simulation and to analytical theory, but these effects are small.) The choice of the geometry to be modeled was motivated by the research reported in ref 31. These authors combined a quartz crystal microbalance with spectroscopic ellipsometry with the goal to, first, determine the amount of coupled water and to, second, use this information to discriminate between clustered aggregates and adsorbates obeying random sequential adsorption.32 The protein Annexin A5, adsorbed to a supported

Figure 2. Sketch of the two geometries under consideration. The resonator surface is located at the bottom. The resonator surface undergoes a tangential time-harmonic oscillation. Periodic boundary conditions are applied on the four sides. The upper boundary of the cell is a wall (transparent yellow), the velocity of which is updated at each simulation step such that the tangential stress equals the product of the tangential velocity and the liquid’s shear-wave impedance. Two configurations of adsorbed cylinders are compared, which are individual cylinders randomly placed on the surface (left) and cylinders adsorbed in clusters (right). Clustering changes the flow pattern (double arrows).

which is a wall undergoing transverse oscillation. In a semiinfinite Newtonian liquid, the resonator surface launches a plane shear wave, which decays within at distance roughly equal to its wavelength. This case is well understood. If the surface is covered with adsorbate particles, these distort the shear wave and induce a nontrivial stress distribution at the resonator surface. It is this stress distribution which needs to be calculated in order to link the properties of the sample to the shifts of frequency and bandwidth as experimentally determined with a QCM. The geometry chosen here is simple in many regards, but still of considerable practical relevance. As first formulated by C

DOI: 10.1021/acs.analchem.5b01912 Anal. Chem. XXXX, XXX, XXX−XXX

Article

Analytical Chemistry lipid bilayer (SLB) was selected as the model system. The SLB was formed on the resonator surface following standard protocols33 prior to the adsorption experiment. Annexin A5 has cylindrical shape with a diameter of about 5.4 nm and a height of 2.7 nm (values taken from Table 2 in ref 31). The adsorption of Annexin A5 to SLBs has been studied in detail before.34 AFM images show Annexin A5 to adsorb face-on. For the work reported in ref 31, Annexin A5 was employed in native form and in a genetically modified form. Under the conditions chosen, the native form mostly adsorbs in the form of trimers. The genetically modified form lacks certain interaction sites, and the molecules therefore adsorb as single units. As expected, the adsorbates of the native form (trimers) were found to contain less coupled water than adsorbates of the mutant (adsorbed as individual particles). While this seems plausible, the Lattice−Boltzmann simulation reveals that trimers should contain more coupled water than dimers and, also, that the experimentally determined differences in coupled water between single particles and trimers match the differences found by simulation between single particles and dimers (not trimers). Clearly, the amount of trapped water depends on the details of the geometry.

hn =

w0 =

(6)

12 2 1 , wn = for n = 1...6, wn = for n = 7...18 36 36 36

Somewhat loosely, we call hn “populations” in the following. The more correct term would be “normalized deviations from the populations at global equilibrium”. Formulating the model in terms of hn rather than f n does not entail any assumptions. It only simplifies the algebra. The amplitude of the local velocity, uE, is related to the populations as

uE , i =

∑ wnhncn,i

(8)

n

The index i labels the direction in space. The deviation of the density from its equilibrium value is given by

δρ =

∑ wnhn

(9)

n

Apart from the global equilibrium at vanishing velocity, there also are local equilibria, accounting for finite velocity and a nonzero δρ. The local equilibrium values, hn,eq, are given as uE , i δρ hn , eq = + ∑ 3cn , i ρ0 cu (10) i

(4)

ρ0 is the density at equilibrium. Note again that all populations are linear functions of uE and δρ. With boundary conditions included, the simulation iterates three steps, which are “streaming”, “bounce-back”, and “relaxation”. In the streaming step, all values of hn are shifted to the neighboring lattice site corresponding to the respective velocity. The bounce-back step accounts for the boundary conditions and occurs after the streaming step (see eq 14 below). Relaxation occurs according to the equation

The combination of parameters ν/cs2 is a relaxation time. The additive term of 1/2 is a consequence of the discrete lattice. The kinematic viscosity was chosen as ν = 10−6 m2/s (resulting from a dynamic viscosity of η = 1 mPa s and a density of ρ = 103 kg/m3). The D3Q19 scheme was employed. The 19 velocities in units of cu are given as c0 = (0, 0, 0)

1 (hn(r , t ) − hn , eq(r , t )) λ 1 + iωΔthn(r , t ) − (ωΔt )2 hn , eq(r , t ) 2

hn(r , t + Δt ) = hn(r , t ) −

c1 = (1, 0, 0), c 2 = ( −1, 0, 0), c3 = (0, 1, 0), c4 = (0, −1, 0), c5 = (0, 0, 1), c6 = (0, 0, −1),

(11)

c 7 = (1, 1, 0), c8 = (1, −1, 0),

(The left-hand side of eq 11 looks different from the corresponding equation in ref 23 because eq 11 does not include the streaming step, while the corresponding equation in ref 23 includes it.) The first two terms on the right-hand side in eq 11 are the same as in conventional LBM. The two last terms account for the oscillatory nature of the flow field. Technically speaking, these two terms are not part of the relaxation process. They are only incorporated into this step for the sake of computational convenience. For the derivation of eq 11, see ref 23. Critical to the success of the simulation is an accurate determination of the tangential stress at planar interfaces. The stress at a planar boundary is equal to the density of momentum flow in the plane adjacent to the surface, which is35

c 9 = (1, 0, 1), c10 = (1, 0, −1), c11 = ( −1, 1, 0), c12 = ( −1, −1, 0), c13 = ( −1, 0, 1), c14 = ( −1, 0, −1), c15 = (0, 1, 1), c16 = (0, 1, −1), c17 = (0, −1, 1), c18 = (0, −1, −1)

−1

(7)

DETAILS OF THE FD-LBM ALGORITHM The time step (Δt) and the unit length of the grid (Δr) are linked to each other by the unit velocity, cu. As in conventional Lattice−Boltzmann methodology, cu is chosen as 31/2cs, where cs is the speed of sound (the speed of compressional waves, cs = 1500 m/s in water). The kinematic viscosity, ν, enters the calculation in the form of a normalized relaxation time, λ, given as 1 1 ν + 2 Δt cs2

fn , globeq

The index n = 0...18 labels the different velocities. The values of hn are complex. The denominators are the population densities in global equilibrium, which are equal to the weight functions, wn, from standard LBM. In the D3Q19 scheme, the weight functions are

■

λ=

fn

(5)

Differing from conventional LBM, the populations corresponding to the various velocities are not quantified by probabilities, f n, but rather by deviations of these probabilities from their equilibrium values at vanishing velocity, f n,globeq. The deviations are called hn and are defined as D

DOI: 10.1021/acs.analchem.5b01912 Anal. Chem. XXXX, XXX, XXX−XXX

Article

Analytical Chemistry σij =

cu2ρ 1−

1 2λ

Equation 14 is formulated such that there are flows on both sides of the boundary. Of course, the flow inside the particle has no physical relevance. Keeping this flow in the formulation has the advantage that the momentum transfer at the interface (equal to the force density) can be computed from the populations on both sides of the interface. For the details see ref 24.The index K labels the boundary nodes. The total force onto a particle then follows from summation over the boundary nodes. Unfortunately, there are a few intricacies connected to the determination of the force onto particles with the bounceback scheme. The main problem is that the result is not accurate enough to be used in the calculation of Δf. (The problem is less severe in the context of particles dynamics.) The bounce-back scheme lacks what is called “second-order accuracy”. This contrasts to eq 13, which does have secondorder accuracy. When simulating plane waves (section I in the Supporting Information), the stress as determined with eq 13 is found to agree with the analytical result within about 10−3. Converted to frequency shift, the accuracy is around 0.5 Hz, which happens to be the accuracy actually needed for the comparison with experiment. Unfortunately, an accuracy at this level is not reached when measuring forces using the bounceback scheme. When particles are adsorbed to the resonator surface, the area-averaged stress cannot be computed with eq 13 because parts of the resonator surface are not actually exposed to the fluid phase. In principle, one might add the transverse force exerted onto the surface-adsorbed particles into the calculation of the area-averaged stress. However, the lack of accuracy in the determination of the force from the bounceback scheme spoils this method of determining ⟨σxz⟩area. The calculation of Δf relies on the small dif ference in stress between the two states with and without adsorbed particles. A 10%-error in the calculation of the force onto the particles translates into a much larger error in the calculation of this small difference. A precision of the order of 10−3 is needed for the calculation of the stress at the resonator surface. It turns out that equivalent information can also be obtained from the stress at the upper cell boundary, if the sample is acoustically thin, if the sample is structured on a scale much below the wavelength of sound, and if the sample has a density close to the density of the liquid. Of course, the upper boundary must then be void of particles, which is no problem. We describe this method of calculating Δf in the next section. There are representations of moving interfaces, which are more accurate than the bounce-back scheme.36,37 One of these was incorporated into the FD-LBM scheme in ref 38. Application of these concepts to the simulation of QCM experiments is left for future work. The validation steps reported in the Supporting Information prove that bounce-back combined with determination of Δf from the stress at the opposing cell wall is a reasonable option for a start. It should be emphasized that the code underlying the current simulation is extraordinarily simple. The core of the program (written in Delphi) contains less than 200 lines, strictly following eq 5 to eq 14 above. Applied to thin layers, the agreement with the analytical results and the expectations following from the Sauerbrey equation is impressive.

∑ cn,icn,jwnhn (12)

n

The term in the denominator accounts for incomplete equilibration in the relaxation step. Applied to the shear stress at the resonators surface, eq 12 reads as ⟨σxz⟩area =

cu2ρ 1−

1 2λ

1 ⟨h9 + h14 − h10 − h13⟩area 36

(13)

The area-averaged stress must be computed because this is the essential input to the small-load approximation (eq 18) on the right-hand side. Periodic boundary conditions were applied at the four sides of the simulation volume. All other boundaries were accounted for by the bounce-back scheme, following ref 24. Bounce-back occurs at the boundary nodes (circles in Figure 3). For an

Figure 3. Geometry underlying the bounce-back scheme. The particle was depicted as a truncated sphere (rather than a cylinder as in Figure 2). Squares denote the lattice sites; small circles denote boundary nodes. To the simulation, the particle is represented by the boundary nodes. Changing the resolution of the grid can therefore change the effective particle size and shape (see also Figure S3A in the Supporting Information).

interface at rest, the populations corresponding to velocities pointing away from the interface (index n) are replaced by the populations traveling in the opposite direction (index n′). If the interface moves, the term 6/cu ∑i(uicn,i) is added, which results in u hn ′(r , t + 1) = hn(r + cnΔt , t ) − 6 ∑ cn , i i c u i hn(r + cnΔt , t + 1) = hn ′(r , t ) + 6 ∑ cn , i i

ui cu

(14)

ui is the i-th component of the velocity of the boundary node. The bounce-back step is carried out after the streaming step. The population at r + cnΔt (first term on the right-hand side in line 1) has moved to this site from the site at r in the step preceding bounce-back. For that reason, hn′ only depends on the flow outside the boundary. Likewise, the artificial flow inside the boundary (second line in eq 14) is decoupled from the outside flow. To see that the factor of 6 must be applied on the right-hand side, consider a homogeneous flow with an immersed particle moving at the same velocity as the ambient fluid. For such flows, all populations are in local equilibrium. Following eq 10 we have hn = 3/cu ∑i(uicn,i). hn′ is equal and opposite to hn. If hn shall be equal to 3/cu ∑i(uicn,i) after bounce-back, the term 6/cu ∑i(uicn,i) must be added to hn′. This term in effect reverses the sign, as it must be.

■

IMPEDANCE BOUNDARY CONDITION AT THE OUTER BOUNDARY, DETERMINATION OF Δf AND ΔΓ FROM THE STRESS AT THE OUTER BOUNDARY We will later be interested in adsorbed objects with a size of a few nanometers. A mesh size of Δr < 1 nm is therefore needed. E

DOI: 10.1021/acs.analchem.5b01912 Anal. Chem. XXXX, XXX, XXX−XXX

Article

Analytical Chemistry If the simulation cell is made to cover a z-range corresponding to a multiple of the penetration depth (the latter being 250 nm for 5 MHz resonators in water), a large part of the simulation effort is invested into the homogeneous liquid. One can make more efficient use of the computational resources by placing the outer boundary closer to the sample and oscillating this boundary tangentially, such that it emulates a Newtonian liquid.39 One updates the velocity of this wall in each simulation step, such that the ratio of stress and velocity equals the shearwave impedance of the liquid Zliq = (iωρη)1/2. A complication arises from the fact that the stress-velocity ratio must take the correct value af ter bounce-back. Going through the calculation, one finds Zliquupb = σ13 =

must be interpreted with care. The most critical parameter is the grid resolution. While grid resolution is critical, it cannot be easily enhanced to better than Δr = 0.25 nm because both memory space and simulation time scale as Δr−3. Test of the Sauerbrey Model. The Sauerbrey equation predicts ΔΓ to be zero and Δf/n to be the same on all overtones. This prediction has been tested in Figure 4. Panel A

cu2ρ 1 1 − 1/(2λ) 36

× (2h9 − 6uupb − (2h13 + 6uupb)) =

h9 − h13 cu2ρ cu2ρ 1 − uupb 1 − 1/(2λ) 18 1 − 1/(2λ) 3

(15)

From this relation, the tangential velocity of upper boundary, uupb, follows as ⎛ 1 cu2ρ 1⎞ uupb = ⎜⎜ + ⎟⎟(h9 − h13) 6⎠ ⎝ Zliq 18(1 − 1/(2λ))

(16)

Equation 16 invokes an approximation insofar as the flow is assumed to be in the xy-plane at the boundary. If the outer boundary is close to the sample, the deviations from in-plane flow will be noticeable and the simulation results will depend on the height of the simulation volume (see Figure S3 in the Supporting Information). More elaborate boundary conditions can be used to account for such effects,40 but for the simulations reported here, we rather removed the outer boundary far enough from the sample to let nonplanar flows at the boundary be negligible. Impedance boundary conditions at z = zmax carry a second benefit unrelated to computational efficiency, which is the possibility to infer Δf and ΔΓ from the stress at this boundary. Why this is necessary was explained at the end of the Details of the FD-LBM Algorithm section. The equation to be applied (proven in the Supporting Information) is ⎛ uupb ⎞ Δf + iΔΓ i ≈ − 1⎟⎟ Zliq⎜⎜ πZq ⎝ uupb , ref f0 ⎠

Figure 4. Test of the predictions of the Sauerbrey model. Panel A shows the normalized frequency shifts Δf/n for harmonics 1, 3, 5, 7, and 9. The legend gives the corresponding resonance frequencies. The differences between the overtone orders are hardly discerned. Panel B shows the shifts in bandwidth. These are negative, but very small. This result should be viewed as an artifact. In panel C, the results from A have been normalized to the result at n = 1 and plotted versus n. There indeed are small differences between the overtones, but their significance is uncertain.

shows the normalized frequency shifts at harmonics n = 1 to n = 9 versus coverage. Since the differences between overtones cannot be discerned, they are again displayed in panel C in normalized form. All values have been divided by the value at n = 1. The normalized values are not strictly the same; there is variability at the level of below 1%, corresponding to frequency shifts in the range of 0.1 Hz. Likewise, the values of ΔΓ are not strictly zero. The values are negative and of the order of magnitude of the order of 0.1 Hz. Neither the variability in Δf/n nor the finite values of ΔΓ are significant because errors of the same magnitude are found for the semi-infinite liquid (section SI in the Supporting Information). It is still interesting to compare the simulation results to analytical predictions of these parameters from ref 6. Ref 6 is concerned with random rough surfaces. (The geometry studied here is not strictly a random rough surface in the sense of ref 6. However, periodicity is invoked at a scale larger than the single-particle level. There is a limited degree of randomness.) For the bandwidth, ref 6 predicts

(17)

uupb,ref is the amplitude of motion in the absence of the sample. Again, eq 17 holds for thin samples with in-plane structure on a scale much below the wavelength of sound and with a density close to the density of the liquid.

■

RESULTS AND DISCUSSION We have gone through a few steps of validation. The details are provided in the Supporting Information (section SI). There are two classical test cases for simulations of QCM experiments, which are the plane shear wave in a semi-infinite Newtonian liquid and the rigid thin film. We go through these tests in sections I.1 and I.2 of the Supporting Information. Section I.3 is concerned with the limitations resulting from finite grid resolution, finite cell height, and finite ringing-in time. The accuracy achieved is about 0.5 Hz. The simulations do show systematic trends at levels below 0.5 Hz. Clearly, these

ΔΓ 1 ≈ πZq f0

2 ρωη hr 2 2 2 δ

(18)

hr is the vertical scale of roughness, assumed to be much smaller than the penetration depth, δ. Choosing hr as 2 nm for the sake F

DOI: 10.1021/acs.analchem.5b01912 Anal. Chem. XXXX, XXX, XXX−XXX

Article

Analytical Chemistry

Figure 5. Frequency shift and fractional trapped solvent versus coverage for single particles and for clusters containing 2, 3, or 4 particles. Error bars are errors of the mean from 10 simulation runs with randomly positioned particles. “Clustered” particles have an interparticle distance of exactly twice the cylinder radius. Trimers (center) form equilateral triangles; tetramers (right) form two such triangles connected to each other at one side. The fractional trapped solvent is obtained by converting Δf to a volume using the Sauerbrey equation (assuming ρparticle = ρliq) and dividing this volume by the volume of the particles themselves. The straight lines in panels B and D are the results for singlets (upper, black line) and trimers (lower, red line) from Figure 2 in ref 31, where full coverage was identified with an optical mass of 360 ng/cm2.

of comparison, one finds ΔΓ ∼ +0.1 Hz n1/2. Clearly, the sign is opposite to what was seen in the simulation (Figure 4B). Why this discrepancy arises is unclear at this point. Again, the absolute values are close to the limit of what experiment can resolve. Reference 6 also makes predictions on the overtones dependence of Δf/n, but the details of the argument are so complicated that the comparison is not worthwhile. Details aside, Figure 4 shows that the layer under study follows the Sauerbrey model within 0.2 Hz. Influence of Clustering on the Amount of Coupled Water. We now come to the example chosen to demonstrate the usefulness of the technique, which is the amount of coupled water and its dependence on the lateral organization of adsorbed particles.31 Richter and co-workers divide the coupled mass by the true mass of the adsorbate and arrive at a dimensionless parameter, which they call “fractional trapped solvent”. Typical values are between 0.5 and 0.9. The fractional trapped solvent is a parameter specific to the structure of the adsorbate. Its interpretation is subject to debate. In an earlier paper, Richter and co-workers calculated the expected amount of fractional trapped solvent from a simple geometric model. They equate the effective volume of coupled water to the volume inside a cone around the particle in question.41 For the sake of this simple model, all particles are represented as cylinders. The height of the cylinder and the angle of the cone are free parameters of the model. They are determined by fitting the dependence of optical and acoustic mass on coverage with the results from the model. Once the apparent height of the cylinders and the cone angle are fixed, the model can be expanded to predict the consequences of clustering. The line of reasoning is that clustering will reduce the amount of coupled water because neighboring cones of coupled water overlap. While this model has many free parameters, it makes one prediction which is not subject to fitting: Comparing single particles, dimers, trimers, and tetramers, the amount of coupled water decreases with cluster size. This hypothesis was tested with the calculations shown in Figure 5. As in the previous calculations, Δr and hC were 0.25 and 10 nm, respectively. The ringing-in time was 0.22 ns. The width of the cell, wC, was

adjusted such that wC2 was equal to the area covered by particles divided by the desired coverage. The principal source of uncertainty is the limited grid resolution. All calculations were carried out 10 times, randomly varying particle positions. Note: What was randomly varied were the positions of 2, 3, or 4 particles. This procedure is to be distinguished from random sequential adsorption because of the periodic boundary conditions at the sides of the cell. Random sequential adsorption uses much larger particle ensembles.32 The results of these simulations are shown in Figure 5. Only Δf as calculated on the fundamental at 5 MHz is shown. Changes in bandwidth and the dependence of Δf/n on overtone order were negligible. The frequency shift in all three cases (dimers, trimers, tetramers) decreases when comparing clustered particles to randomly positioned particles. Interestingly, the decrease is strongest for the dimers. It is hardly noticeable for dimers and tetramers. These results impressively demonstrate that hydrodynamic effects are not easily captured by simple models. The “intersection of cones” invoked in ref 41 does not capture the entire picture. When comparing dimers to trimers the volume affected by the hydrodynamic interaction between the liquid and the neighboring particles decreases, but the strength of this coupling is strong in the space between the members of a trimer because there is no escape route to the side. Figure 6 sketches the essence of this argument. The water contained in the interstitial area between three particles is

Figure 6. Sketch illustrating why the efficiency of hydrodynamic coupling of liquid to neighboring particles is more efficient for trimers than for dimers. The water contained in the interstitial space between the members of a trimer can only escape to the top (color online). G

DOI: 10.1021/acs.analchem.5b01912 Anal. Chem. XXXX, XXX, XXX−XXX

Analytical Chemistry

Article

■

CONCLUSIONS A method was described which simulates QCM experiments, based on a frequency domain version of the Lattice Boltzmann method. The code is simple; it can run on office PCs. Also, it can be parallelized in order to increase the size of the simulation volume or to improve the resolution of the grid. For the semi-infinite liquid, the simulation reproduces the analytical predictions within an accuracy of about 0.1%. This accuracy suffices for the comparison with experiment. The current implementation of the method represents embedded particles with the bounce-back scheme. Since the bounce-back scheme lacks second order accuracy, the stress at the resonator surface (needed in the calculation of Δf) can no longer be evaluated with sufficient accuracy at the resonator surface. For acoustically thin layers with a density equal to the density of the liquid, equivalent information can be obtained from the stress at the outer boundary of the cell. Applied to the problem of the apparent mass of structured adsorbates in the liquid phase, the simulations show that details of the geometry do matter. Clustering of adsorbed cylinders reduces the amount of coupled liquid, but contrary to intuition, the effect is substantially stronger for dimers than for trimers. This can be tentatively explained with the fact that liquid contained in the interstitial space between trimers can only escape to the top, while the liquid around dimers can also escape to the side. This caging effect increases the efficiency of hydrodynamic coupling.

tightly coupled to its neighbors. On a more fundamental level, the cylinders-and-cone model attempts to capture the situation by volume-considerations only. It neglects the efficiency of coupling, which depends on the details of the molecular arrangement. The results shown in Figure 5 are not only at variance with the model from ref 41 but also at variance with experiment, which leaves the question what features of the experiment are not captured by the FD-LBM simulation. A first guess is an unreasonably simple representation of clusters as cylinders just barely touching each other. Representing the protein as a cylinder is a simplification in the first place. When proteins aggregate, there are parts at the rim, which fit well to each other, and there are interactions, which will deform the protein. In all likelihood, the size of a trimer is less than what results from arranging undeformed cylinders in triangular form. Figure 7 shows results from calculations looking into this hypothesis. Single particles were compared to trimers, but

■

ASSOCIATED CONTENT

* Supporting Information S

Certain steps of validation of the code as well as the proof of eq 17. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.analchem.5b01912.

Figure 7. Frequency shift and fractional trapped solvent for assemblies of three particles, which were not clustered at all, were clustered as trimers, or were clustered as trimers with overlap between neighboring cylinders, such that the distance between the particle centers was lowered by 8% or 20%. With a contraction of 20%, the simulation comes close to the experimental result (straight lines in Figure 5D).

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: + 49-532372 3768. Fax + 49-5323-72 2863. Notes

The authors declare no competing financial interest.

■ ■

differing from Figure 5, the interparticle distance was not equal to twice the cylinder radius, but was rather lowered relative to this value by an amount indicated in the legend. A contraction by about 20% is needed in order make the simulation results match the experiment (straight lines in Figure 5D). One should interpret this number with some care, keeping the finite resolution of the grid in mind. Also, adsorption in experiment occurred to a slightly flexible substrate (a supported lipid bilayer, SLB), which also may be of influence. A final word of caution: The FD-LBM algorithm is strong when the experimental configuration offers a clear-cut hypothesis on the geometry or offers clear-cut alternatives between different geometrical configurations. These can then be translated to predictions on Δf and ΔΓ, thereby allowing for comparison with experiment. If there are no such clear-cut alternatives, it will always be difficult to run the formalism backward. Many different geometries can lead to the same values of Δf and ΔΓ; inferring geometry from QCM-data alone without further constraints is impossible. This limitation must be kept in mind for all models of the kind described above.

ACKNOWLEDGMENTS We thank R. P. Richer (CIC Biomagune, San Sebastian) for stimulating discussions. REFERENCES

(1) Lu, C.; Czanderna, A. W. Applications of Piezoelectric Quartz Crystal Microbalances; Elsevier: Amsterdam, 1984. (2) Sauerbrey, G. Eur. Phys. J. A 1959, 155 (2), 206−222. (3) Steinem, C.; Janshoff, A. Piezoeletric Sensors; Springer: Heidelberg, 2007. (4) Lucklum, R.; Hauptmann, P. Anal. Bioanal. Chem. 2006, 384 (3), 667−682. (5) Johannsmann, D. Phys. Chem. Chem. Phys. 2008, 10 (31), 4516− 4534. (6) Daikhin, L.; Gileadi, E.; Katz, G.; Tsionsky, V.; Urbakh, M.; Zagidulin, D. Anal. Chem. 2002, 74 (3), 554−561. (7) McHale, G.; Newton, M. I. J. Appl. Phys. 2004, 95 (1), 373−380. (8) Johannsmann, D.; Reviakine, I.; Rojas, E.; Gallego, M. Anal. Chem. 2008, 80 (23), 8891−8899. (9) Olsson, A. L. J.; Quevedo, I. R.; He, D.; Basnet, M.; Tufenkji, N. ACS Nano 2013, 7 (9), 7833−7843. H

DOI: 10.1021/acs.analchem.5b01912 Anal. Chem. XXXX, XXX, XXX−XXX

Article

Analytical Chemistry (10) Berglin, M.; Olsson, A.; Elwing, H. Macromol. Biosci. 2008, 8 (5), 410−416. (11) Olsen, E. V.; Pathirana, S. T.; Samoylov, A. M.; Barbaree, J. M.; Chin, B. A.; Neely, W. C.; Vodyanoy, V. J. Microbiol. Methods 2003, 53 (2), 273−285. (12) Wegener, J.; Janshoff, A.; Galla, H. J. Eur. Biophys. J. 1998, 28 (1), 26−37. (13) Keller, C. A.; Kasemo, B. Biophys. J. 1998, 75 (3), 1397−1402. (14) Zhang, X. H. Phys. Chem. Chem. Phys. 2008, 10 (45), 6842− 6848. (15) Johannsmann, D.; Mathauer, K.; Wegner, G.; Knoll, W. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 46 (12), 7808−7815. (16) Reviakine, I.; Johannsmann, D.; Richter, R. P. Anal. Chem. 2011, 83 (23), 8838−8848. (17) Rodahl, M.; Hook, F.; Krozer, A.; Brzezinski, P.; Kasemo, B. Rev. Sci. Instrum. 1995, 66 (7), 3924−3930. (18) Johannsmann, D. Macromol. Chem. Phys. 1999, 200 (3), 501− 516. (19) Part of the Multiphysics Package by COMSOL, http://www. comsol.com/. (20) Pomorska, A.; Shchukin, D.; Hammond, R.; Cooper, M. A.; Grundmeier, G.; Johannsmann, D. Anal. Chem. 2010, 82 (6), 2237− 2242. (21) Succi, S. The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond; Oxford Science Publications: 2013. (22) Qian, Y. H.; Dhumieres, D.; Lallemand, P. Europhys. Lett. 1992, 17 (6BIS), 479−484. (23) Shi, Y.; Sader, J. E. Phys. Rev. E 2010, 81, 036706. (24) Groot, S. R. D.; Mazur, P. Non-Equilibrium Thermodynamics; Dover Publications: 1985. (25) Ladd, A. J. C. J. Fluid Mech. 1994, 271, 285−309. (26) Johannsmann, D. The Quartz Crystal Microbalance in Soft Matter Research: Fundamentals and Modeling; Springer: 2014. (27) Hook, F.; Kasemo, B.; Nylander, T.; Fant, C.; Sott, K.; Elwing, H. Anal. Chem. 2001, 73 (24), 5796−5804. (28) Kosslinger, C.; Uttenthaler, E.; Drost, S.; Aberl, F.; Wolf, H.; Brink, G.; Stanglmaier, A.; Sackmann, E. Sens. Actuators, B 1995, 24 (1−3), 107−112. (29) Domack, A.; Prucker, O.; Ruhe, J.; Johannsmann, D. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1997, 56 (1), 680−689. (30) Plunkett, M. A.; Wang, Z. H.; Rutland, M. W.; Johannsmann, D. Langmuir 2003, 19 (17), 6837−6844. (31) Carton, I.; Brisson, A. R.; Richter, R. P. Anal. Chem. 2010, 82 (22), 9275−9281. (32) Adamczyk, Z.; Siwek, B.; Zembala, M.; Belouschek, P. Adv. Colloid Interface Sci. 1994, 48, 151−280. (33) Richter, R. P.; Berat, R.; Brisson, A. R. Langmuir 2006, 22 (8), 3497−3505. (34) Richter, R. P.; Him, J. L. K.; Tessier, B.; Tessier, C.; Brisson, A. R. Biophys. J. 2005, 89 (5), 3372−3385. (35) Krueger, T.; Varnik, F.; Raabe, D. Phys. Rev. E 2009, 79, 046704. (36) Latt, J.; Chopard, B.; Malaspinas, O.; Deville, M.; Michler, A. Phys. Rev. E 2008, 77, 056703. (37) Niu, X. D.; Shu, C.; Chew, Y. T.; Peng, Y. Phys. Lett. A 2006, 354 (3), 173−182. (38) Shi, Y.; Yap, Y. W.; Sader, J. E. Phys. Rev. E 2014, 89, 3. (39) Sun, C.; Pérot, F.; Zhang, R.; Freed, D. M.; Chen, H. Commun. Comput. Phys. 2013, 13 (3), 757−768. (40) Mecke, K. R. J. Phys.: Condens. Matter 2001, 13 (21), 4615− 4636. (41) Bingen, P.; Wang, G.; Steinmetz, N. F.; Rodahl, M.; Richter, R. P. Anal. Chem. 2008, 80, 8880−8888.

I

DOI: 10.1021/acs.analchem.5b01912 Anal. Chem. XXXX, XXX, XXX−XXX