Stiffness of Contacts between Adsorbed Particles and the Surface of a QCM‐D Inferred from the Adsorption Kinetics and a Frequency‐Domain Lattice Boltzmann Simulation

A simulation based on the frequency‐domain lattice Boltzmann method (FreqD‐LBM) is employed to predict the shifts of resonance frequency, Δf, and half bandwidth, ΔΓ, of a quartz crystal microbalance with dissipation monitoring (QCM‐D) induced by the adsorption of rigid spheres to the resonator surface. The comparison with the experimental values of Δf and ΔΓ allows to estimate the stiffness of the contacts between the spheres and the resonator surface. The contact stiffness is of interest in contact mechanics, but also in sensing because it depends on the properties of thin films situated between the resonator surface and the sphere. The simulation differs from previous implementations of FreqD‐LBM insofar, as the material inside the particles is not included in the FreqD‐LBM algorithm. Rather, the particle surface is configured to be an oscillating boundary. The amplitude of the particles' motions (displacement and rotation) is governed by the force balance at the surface of the particle. Because the contact stiffness enters this balance, it can be derived from experimental values of Δf and ΔΓ. The simulation reproduces experiments by the Krakow group. For sufficiently small spheres, a contact stiffness can be derived from the comparison of the simulation with the experiment.


Introduction
The QCM-D as a surface-analytical instrument is attractive not only because it determines the thickness of the layer under study with submonolayer precision, but also because it is sensitive to non-gravimetric effects. [2]In this regard, the applications of the QCM-D fall into two classes.On the one hand, researchers detect minute amounts of adsorbed material and immobilize certain receptors to the QCM surface, thereby adding specificity chanical" may be replaced by "acoustic" for planar layer systems.When the sample is a planar layer system, the QCM works similar to optical reflectometry, where the refractive index profile is replaced by the profile of the material's shear wave impedance (equal to ((z)•G(z)) 1/2 with  the density and G the shear modulus). [12]For structured samples, one should stick to "mechanical", but keep the high frequency in mind.Objects, which are soft at low frequency (meaning, transport little stress at low frequency), are less soft at a few MHz because the fast deformation prevents some of the relaxation processes.Phrased differently, the shear modulus depends on frequency in soft matter.
In these cases, the challenge is not so much the sensitivity, but rather the interpretation.There are two problems, which are, first, the targeted configuration being underdetermined and, second, the lack of models, which would predict the QCM response from a given configuration.This work contributes to the second problem.Within the standard formalism, the QCM exerts a small perturbation onto the sample by forcing the material in contact with the resonator to take part in the oscillatory movement of the surface.The stress onto the resonator surface results from the sample's displacement pattern.The small load approximation states that the movement of the resonator surface is almost independent of the stress.In rheological terms, the sample is under strain control.The parameter to be measured is the stress.It follows that the complex frequency shift, Δf + iΔΓ (with f the resonance frequency and Γ the half bandwidth, both determined on a few different overtones), is given as f 0 (often 5 MHz) is the frequency of the fundamental, ΔD is the shift in dissipation factor (related to the shift in half bandwidth as ΔD = 2ΔΓ/f res ), Z q = 8.8 ⋅10 6 kg (m 2 s) −1 is the resonator's shear wave impedance, < S > area is the area-averaged tangential stress at the resonator surface (a complex amplitude), and v x , S = 1 is the tangential velocity.The stress-velocity ratio is also called load impedance, Z L .For proof of Equation ( 1) see Section 8.3 in Ref. [13] The task is to compute the displacement pattern inside the sample with the oscillating substrate as a boundary condition.Once this solution is found, the transverse oscillatory stress at the surface follows from the shear gradient at the surface.The stress is averaged over the area, and this average is inserted into Equation (1).
When the sample constitutes an inertial load, the load impedance simply equals im f with m f the area-averaged mass per unit area.This situation is characteristic of gravimetry, where the Sauerbrey relation connects the frequency shift to the mass. [14]In gravimetry, the shifts of bandwidth are zero and the overtone-normalized frequency shifts, Δf/n, are the same on all overtones.Anything deviating from this situation is a "nongravimetric effect". [15]Non-gravimetric effects are reflected in the QCM data in two separate forms.First, they increase the half bandwidth, ΔΓ.Second, the non-gravimetric effects cause a dependence of the overtone-normalized frequency shift, Δf/n, on overtone order, n (often n = 3,5,7,9).
For planar films, non-gravimetric effects mostly go back to softness.For sufficiently stiff films, the elastic compliance, J′, is proportional to the ratio ΔΓ/(−Δf).Clearly, the bandwidth is of much importance, as also expressed in the name "QCM-D" for "QCM with dissipation monitoring".One may also infer nongravimetric parameters from the dependence of Δf/n on n.The overtone-normalized frequency shifts are no longer the same on the diffenent n's.However, the dependence of Δf/n on overtone order is less easily analyzed than ΔΓ because the viscoelastic parameters themselves depend on frequency.The work reported below mostly builds on ΔΓ, more specifically on the ratio ΔΓ/(−Δf).
[18] The ratio of the two, ΔΓ/(−Δf), therefore is an intrinsic property of the material, unrelated to thickness.It has become customary to discuss and analyze this ratio even for samples, which are not thin in the above sense.Following Ref. [19] ΔΓ/(−Δf) is called "acoustic ratio" in the following.More precisely, Ref. [19] defines the acoustic ratio as ΔD[10 −6 ]/(−Δf/n), which differs from ΔΓ/(−Δf) by a factor of 2.5 (for 5 MHz resonators).Defining the acoustic ratio as ΔΓ/(−Δf) avoids a dependence on the resonator's fundamental frequency and it makes the acoustic ratio dimensionless.
For planar layers with a thickness much below the wavelength of sound (a few nanometers for soft adsorbates), the viscoelastic parameters of the layer (such as the shear modulus G = Gʹ + iGʹʹ or the shear compliance, J = 1/G = Jʹ -iJʹʹ) can be derived from the set of parameters {Δf/n, ΔΓ/n}, even including an estimate of their dependence on frequency. [20]No such analytical models exist for structured samples.These require a numerical simulation, which predicts Δf/n and ΔΓ/n from the geometry and the viscoelastic parameters in the different domains.The simulation solves the continuum equations of hydrodynamics and viscoelasticity.The solution leads to an area-averaged complex amplitude of the oscillatory transverse stress at the surface, < S > area , to be inserted into Equation (1).
As shown in previous work, [21,22] running the simulation in the forward direction leads to rather reliable results for soft matter.The inversion, namely the determination of a sample's configuration from the sets of Δf/n and ΔΓ/n is problematic.Usually, there will be numerous different configurations of a sample (particle size, particle shape, coverage, softness, …) leading to the same values of Δf/n and ΔΓ/n.Such calculations only provide an insight, if the geometry is well-defined to the extent that the answer to the question posed by the experimentalist can be condensed into few numbers.The contact stiffness plays this role, here.
The working hypothesis underlying this study is that QCM-D data obtained while adsorbing particles to the resonator can be exploited to estimate the stiffness of the contact between the particle and the substrate.The simulation below targets spherical particles, [1,23,24] but similar simulations might address non-spherical nanoparticles, [25] proteins, [26] dendrimers, [27] or viruses. [28]The contact stiffness is of interest for contact mechanics, but it may also be an attractive parameter for sensing.This occurs when the contact is made across layers of (bio-) molecules, the properties of which are the target of the study.Such layers will take an influence on the contact stiffness.Whether they increase or decrease it, is not easily predicted.The Gizeli group has emphasized that spheres (vesicles, in this case) can amplify a QCM response. [29]The amplification with spheres reveals properties of DNA strands (part of this experimental setting), which would otherwise be inaccessible to the QCM.
We briefly digress on the stiffness of a contact between a sphere and a plate as inferred from QCM experiments on adsorbed spheres in air.For spheres in air, Δf and ΔΓ can be calculated analytically (Section 6.2 in Ref. [13]).There is an instructive analogy to the thin, soft film.For the thin, soft film in air, one plots Δf/n and ΔΓ/n versus n 2 and infers the real and the imaginary part of the layer's complex viscoelastic compliance from the slopes in these plots.The analytical model of soft adsorbed spheres in air predicts an analogous dependence of Δf/n and ΔΓ/n on n 2 .The difference simply is that Jʹ − iJʹʹ is replaced by the inverse contact stiffness, (1/)ʹ − i(1/)ʹʹ.Generally, nongravimetric effects go back to softness.The larger the softness of the sample, the larger are the deviations from Sauerbrey-type behavior.For adsorbed spheres in air, the sample's softness is synonymous to the inverse contact stiffness.
One needs to be careful, though, to not naively apply this argument to adsorbed spheres in a liquid.Exploiting this analogy is tempting because plots of Δf/n and ΔΓ/n versus n look similar to the corresponding plots obtained for planar layers.There is a linear dependence of Δf/n on n, the slope is positive, and the slope is large for soft contacts.Also, ΔΓ/n depends linearly on n.The slope is positive and is large for soft contacts.Graphs of this kind are shown in Section S4 (Supporting Information) to Ref. Decreasing Δf/n indicates increasing coverage, but there is no quantitative relation between coverage and −Δf/n.Because the coverage is not known in absolute terms, the acoustic ratio is plotted versus −Δf/n in the lower panels.These graphs allow for a comparison between experiment and simulation (see Figure 3).[22] One might derive apparent values, Jʹ app and Jʹʹ app , from the slopes and calibrate J app against the inverse contact stiffness, using a FreqD-LBM simulation.This representation was not chosen as the basis of the discussion here, because it is potentially misleading.In air, an increase in bandwidth can with confidence be attributed to a dissipative component of the contact stiffness.This is different in liquids.A perfectly rigid, structured sample (such as a rough metal surface) distorts the flow of the liquid such that the flow no longer is a plane wave.This distortion increases the bandwidth, as shown by Urbakh and Daikhin in Ref. [30] Rather than deriving an effective viscoelastic compliance from the QCM data, we plot ΔΓ/(−Δf) versus −Δf/n (see the bottom graphs in Figure 1) and search for agreement between experiment and simulation based on these plots.Expressed differently, the discussion hinges on the acoustic ratio, similar to what is done in the context of the ΔΓ−Δf extrapolation scheme. [43]his work differs from the previous demonstrations of FreqD-LBM [21,22] in how the particles were treated.For details see Section 3.1.The bulk of hard particles cannot be made part of the FreqD-LBM algorithm.To avoid this problem, the particle surfaces were configured to be oscillating boundaries.Demonstrating the success of this mode of FreqD-LBM is one of the novelties of this work.

Experimental Data to be Reproduced with the Simulation
Figure 1 shows experimental data obtained by the Krakow group, [1] which we attempt to reproduce.Polystyrene spheres were adsorbed to the surface of a resonator, which had been coated with SiO 2 .Details of the experimental procedures are reported in Ref. [1] Polystyrene is almost rigid at room temperature.The -potentials of the spheres were positive (between 70 and 80 mV).The ionic strength of NaCl electrolyte was 0.01 m.The SiO 2 surface was negatively charged at the pH of the experiment (pH 4,  ≈ −16 mV).The spheres' positive charge favors adsorption to the substrate, while at the same time preventing a contact between neighboring adsorbed particles.Apart from the sphere size (diameters of 26, 67, and 140 nm), the experimental parameters were similar in the different runs.
The top panels in Figure 1 show the adsorption kinetics.Data from the fundamental were discarded, because these often are contaminated by effects linked to compressional waves. [31]f/n is always negative.This result is nontrivial because large spheres may also increase the resonance frequency in case they are clamped in space by inertia. [32]ΔΓ/n is always positive, as it must be.Converting the diameter of the particles to a negative overtone-normalized frequency shift with the Sauerbrey equation results in values of −Δf/n = 150, 380, and 800 Hz, for the diameters of 26, 67, and 140 nm, respectively.The experimental values in all cases stay below these limits, indicating incomplete coverage.This can be expected, given the repulsive interaction between particles.
The bottom panels in Figure 1 show the same data as the top panels, changed to a representation which allows for easier comparison with the simulation.The simulation varies the coverage, but it does not attempt to reproduce the adsorption kinetics (i.e., the dependence of coverage on time).Also, the coverage itself cannot be easily inferred from −Δf/n because of liquid trapped between particles (increasing −Δf/n). [33]Substituting the coverage by −Δf/n avoids this problem.
Clearly, the outcomes of these three experiments were rather different.For the largest spheres with a diameter of 140 nm, the acoustic ratio was larger than unity.For the smallest spheres (d Sph = 26 nm), the acoustic ratio decreased with increasing coverage, while the decrease was moderate for the intermediate size (d Sph = 67 nm).
In work, we only attempt to reproduce ratios between ΔΓ and −Δf.Absolute values would also be interesting in the context of a calibration of −Δf/n against coverage.Such a calibration is outside the scope here.

General
In the first step, the simulation solves the continuum equations of hydrodynamics.In principle, this may be achieved with any of the methods of computational fluid dynamics.In passing, we note two analytical models.Tarnapolsky et al. have put forward a set of equations based on an equivalent mechanical circuit, which treats individual soft spheres. [34]This model is geared to bacterial adsorption.It takes both translation and rotation into account.The Krakow group has proposed models, which take the interaction between different particles into account. [35,36]The interaction is modeled with a screening function, which was originally developed to predict how particles adsorbed to a charged surface modify the streaming potential.A hybrid model using both numerical methods (the Finite Volume Method) and an advanced analytical formalism is reported in Ref. [37] Leshansky and coworkers have put forward an analytical treatment of particles above a QCM surface. [38]This model does not yet cover contact mechanics, but may be extended to do so.
The authors' group has used the Finite Element Method (FEM) in the past, but this code only produced robust results in 2D. [39][42] These algorithms were computationally expensive and the grid resolution was problematic.In previous work, we have shown that the frequency domain lattice Boltzmann method (FreqD-LBM) can amount to a detailed and realistic model for sufficiently soft samples.In this formalism, "particles" are domains with a (complex) viscosity larger than the bulk viscosity,  bulk .FreqD-LBM gave support to the ΔΓ-Δf extrapolation scheme. [43]It also explained a transient maximum in ΔΓ observed in QCM experiments on particle fouling. [22]he algebraic details of FreqD-LBM (first proposed by Shi and Sader in Ref. [44]) have been described in Refs.[21] and [22].They are repeated in the supporting information.The populations from standard LBM are replaced by complex amplitudes of oscillation.Because FreqD-LBM assumes that frequency-domain calculations are equivalent to time-domain calculations, the underlying equations must be linear in oscillation amplitude.Linearity is very well obeyed in experiment.The Reynolds number in QCM experiments is of the order of 10 -3 . [45]Nonlinear terms in the collision step therefore can be (and are) omitted.The liquid is quiescent at t < 0. At t = 0, a transverse periodic motion of the substrate is turned on.The substrate launches a wave into the sample.The simulation continues until a stationary state is reached.In the absence of particles, the "ring-in" process simulated in this way is realistic.The ring-in is not of much practical interest, though, because there is no sudden ring-in of this kind in experiment.The interpretation of the simulation results entirely rests on the periodic flow field eventually achieved.If the sample is a Newtonian liquid (no adsorbed particles), the displacement pattern eventually turns into a decaying shear wave.The wave's decay depth, , is given as  = (2/()) 1/2 . is the In the implementation chosen here, the FreqD-LBM simulation is limited to the bulk liquid.The particles are walls with a certain (local) oscillatory velocity, computed in each simulation step from the particle's velocity and its rate of rotation.In the steady state, the force and the torque exerted by the bulk liquid are equal and opposite to the force and the torque exerted by the contact (ignoring a small contribution from inertia).
dynamic viscosity and  is the density.This situation also carries the name "2nd Stokes problem".For 5 MHz resonators in water,  is in the range of 100 nm.
In Refs.[21] and [22] (treating soft particles), the spheres were modeled as domains with an increased (complex) viscosity.This variant of FreqD-LBM runs into problems, when the viscosity of the particle, | P |, is larger than ≈100  bulk .The relaxation rate in the collision step of the simulation then is so small that the populations traverse the domain in question with insufficient equilibration (Section S1.8, Supporting Information to Ref. [22]).This problem can be avoided by not letting the material inside the particles be part of the simulation.Rather, the surfaces of the particles are configured to be boundaries, which oscillate with some complex amplitude (Figure 2).The boundary conditions are implemented in the form of bounce back. [46]For the sake of a more accurate representation of the geometry, interpolation was applied, [47] following Section 11.2.2.1 in Ref. [48]  A moving boundary implies that the surface undergoes an oscillatory motion, similar to the bulk fluid.No update of the geometry (which would be required for steady motion) is needed.The problem is to find the complex amplitudes of oscillation for all nodes on the boundary.The condition leading to these amplitudes is a force balance between the liquid drag, on the one hand, and the forces exerted by the contact, on the other.Inertial forces play a role, in principle, but are negligible in magnitude.The latter statement can be corroborated with a calculation of the momentum relaxation time,  p = m/ fric with m the mass and  fric = 6r Sph the friction coefficient.With a density of 1 g cm −3 a sphere radius of r Sph = 100 nm, and a viscosity of  = 10 -3 Pa s, the momentum relaxation time comes out to be 2 ns.Inertial forces are so weak that the particles adopt the velocity of their environment on a time scale much smaller the period of oscillation.

Dynamics of the Sphere
In order to formulate the force balance condition, the particle's response to the forces exerted by the liquid must be specified.
FreqD-LBM flexible with regard to the details of this response.Here, the sphere as a whole was assumed to not undergo internal deformation.Deformation occurred at the contact, only, as expressed by the finite contact stiffness.The model might be extended by a shear deformation.A shear stress might be calculated from the forces acting at the surface and the shear deformation might be computed, using the material's compliance tensor (which might even include anisotropy).Similarly, a nontrivial hydrostatic pressure might be calculated and the particles might expand and contract in response to the oscillatory pressure.These are only meant to be examples.The dynamics of the particle is entirely separate from FreqD-LBM.Numerous other ways of treating the mechanics of the particle are conceivable.
In this work, only the contact was given a finite compliance.This is legitimate even for stiff particles because there is a stress concentration at the contact, leading to a deformation close to the contact. [49]The properties of the contact are parameterized by a shear stiffness,  sh , and a bending stiffness,  bend .The shear stiffness is the ratio of the integrated transverse force to the sphere's lateral displacement, where the latter is evaluated far away from the contact.The Mindlin model predicts the shear stiffness to be [50]  sh = 2G c r c (2) r c is the contact radius.In Mindlin theory, G c is related to the shear moduli and the Poisson numbers of the media making the contact. [51]In the simulation, G c was prescribed by the algorithm as |G c | = 1.5 GPa, the latter being a typical value for a glassy polymer (such as polystyrene).G c was assumed to be independent of frequency.A dissipative component of the contact's response was introduced in the form of a loss tangent (tan( L ) = G c ʹʹ/ G c ʹ).More specifically, we used G c ʹ = 1.5 GPa cos( L ) and G c ʹʹ = 1.5 GPa sin( L ).The parameter tan( L ) took values of 0.3, 1, and 3.The bending stiffness is proportional to the ratio of the integrated torque and the angle of rotation about the contact.Differing from other definitions of the bending stiffness in the literature, we divide this ratio by r Sph 2 . [52]r Sph is the sphere radius, which is assumed to be much larger than the contact radius, r c .Defined this way, the bending stiffness has the same dimensions as the shear stiffness.Following Dominik and Tielens [52] (see also Ref. [34]), we assume the bending stiffness to be related to the shear stiffness as These relations apply in the limit of r c << r Sph .The stiffness under twist (rotation about the surface normal) was chosen to be the same as the bending stiffness.

Further Details of the Algorithm
Periodic boundary conditions were applied at the side walls of the simulation box.At the bottom, the boundary condition was a transverse oscillation with unit amplitude.At the top, an impedance boundary condition was applied as described in Sec- tion S3.1 (Supporting Information) to Ref. [22] see also Section S5 (Supporting Information).The impedance boundary condition at the top is essential to the success of the simulation.It allows to place the upper boundary of the simulation box inside the range of the decaying shear wave.The boundary must be at some distance from the upper edge of the sample, but there is no need for the velocity to be zero at the boundary.Here, the height of the simulation volume was twice the particle diameter.
The simulation volume contained three particles placed randomly.Such simulations were repeated and the values of Δf and ΔΓ were averaged over five such runs.Multiple, randomly positioned particles are needed because a single particle would amount to a regular array, given the periodic boundary conditions.A regular array creates interference effects, which can be seen as periodic oscillations in the dependence of Δf/n and ΔΓ/n on coverage (Section S8, Supporting Information to Ref. [21]).The coverage varied logarithmically in ten steps between 5% and 26%.The coverage is defined as the shadow of spheres divided by the total area.In technical terms, the width of the simulation volume (always containing three particles) was adjusted accordingly.The grid resolution, Δx, was 1/6 of the particle diameter.Section S9 (Supporting Information) shows -using one example -how the derived values of ΔΓ/(−Δf) and the computation time depend on grid resolution.Δx = 1/6 d Sph is a reasonable compromise between accuracy and computational effort.The minimum center-to-center distance between particles was 4/3 of the particle diameter.The factor of 4/3 ensures that the particles never touch.The particles are not expected to touch in the experiment because they are positively charged.The degree of truncation of the sphere at the bottom was 5%, meaning that the sphere center was placed to a height of 95% of the sphere radius.This choice was motivated by estimations based on the DMT model. [49]The contact radius resulting from this choice, r c , is ≈30% of the sphere radius.r c /r Sph ≈ 30% appears as a large number, judging by experience from the macroscale.In the microscale, adhesive forces gain in importance.The (small) degree of truncation is not of much importance to the simulation, because the contact stiffness is chosen as 2r c G c   with   some numerical coefficient between 10 -3 and 5 (see Equation ( 2)).Because the contact stiffness is varied using the factor   , the degree of truncation only affects the geometry (to a small extent).(Figure 3) Initially, the particles are at rest.They start moving in response to the forces exerted by the liquid, once the movement of the substrate is turned on.A natural update step would be to let the particle obey Newton's second law (F = and to update the velocity following u → u + FDt/m.An analogous relation would apply to the rate of rotation and the torque.However, such updates quickly produce a divergence.The update must be slowed down, artificially.The update occurs as u → u +  FDt/m with  a numerical coefficient of the order of 10 -3 The ring-in process then no longer is realistic.This is no problem because the further analysis is based on the stationary state.Because of the slow update, the ring-in − including the update of the spheres' motion − takes longer than the ring-in of the liquid by about a factor of ten.The simulations reported here are computationally more expensive than the simulations on soft particles reported in Refs.[21] and [22].The update of the particle's motion was particularly problematic for small loss tangents.When the contacts were mostly elastic, the spheres oscillated and these oscillations decayed slowly.

Implementation
A typical simulation volume has n x × n y × n z = 15 × 15 × 15 grid points, depending on coverage.With no update of the sphere motion, a ring-in takes about 3 n y 2 steps, where y is the surface normal.This dependence on n y 2 is the consequence of the displacement pattern being a diffusive wave with the kinematic viscosity, , taking the role of the diffusivity.The characteristic time is some multiple of the square of the diffusion length.The factor of n y 2 needed for the ring-in lets the computation time scale as the 5th of the inverse grid resolution, in total.The update of the particle motion leads to a further factor of ≈10.
The code is written in Python and accelerated with Numba.The Numba keyword "parallel" was set to "True".The Message Passing Interface (MPI) was not employed.Parallelization with MPI would speed up the calculation.
Computations were performed on clusters of the Simulation Science Center Clausthal-Göttingen (SWZ).The clusters employed are either composed of 2 x AMD Epyc 7502 CPUs (64 physical cores, 128 logic cores) with 1024 GB RAM or 2 x AMD Epyc 7281 CPUs (32 physical cores, 64 logic cores) with 1024 GB RAM.Because the FreqD-LBM algorithm is more CPU intensive than RAM intensive, speed was limited by the CPU resources.

Results and Discussion
Figure 4 shows plots of ΔΓ/(−Δf) versus −Δf/n for a sphere diameter of 26 nm and a loss tangent of 0.3.Analogous plots for the different sphere sizes and loss tangents are shown in Section S12 (Supporting Information).The values of the parameter   were spaced logarithmically between 10 −3 and 5 (ten values in total).Only five diagrams are shown in Figure 4, where the values of   are indicated on the left-hand side.
The two columns to the left show Δf/n and ΔΓ/n versus coverage.Both −Δf/n and ΔΓ/n increase with coverage.Clearly, both −Δf/n and ΔΓ/n also increase with increasing contact stiffness (increase from top to bottom).The column on the right-hand side shows the same data represented as the acoustic ratio versus −Δf/n (cf. the bottom graphs in Figure 1).The experimental values are shown as dashed lines.
The plot of −Δf/n versus coverage (left in Figure 4) shows the characteristic curvature, indicative of saturation.For the conditions studied here, saturation sets in at a coverage of 0.3.For the soft particles studied in Ref. [21], this value was closer to 0.5.
The agreement between experiment and simulation was never quantitative in the narrow sense.There would be many possibilities to improve the agreement.Most importantly, a dependence of the contact stiffness on frequency is expected and might be included in the model.Also, a distribution of contact stiffnesses between different spheres can be expected.This might be included in the model as well, but extending the model in this way would lead to multiple configurations achieving similar agreement with the experimental data.The model therefore was not tuned to match the experiment more accurately.Even with these choices, certain values of   are closer to the experiment than others.For the data shown in Figure 4   = 0.07 or   = 0.044 lead to better agreement than   = 0.001 (at the top) or   ≥ 0.292 (the two panels at the bottom).
The agreement was quantified using two parameters called  2 Δf/n and  2 ratio .To derive these, −Δf/n and ΔΓ/(−Δf) were averaged over all coverages, n cov , and all overtones, n ovt , for the simulation and the experiment. 2 is defined to be the square of the difference between these values: Again, the values were averaged before calculating the square of the difference.
Figure 5 shows these square-deviations versus contact stiffness, where the contact stiffness was derived from the values of   as indicated on the left in Figure 4 with the relation  Sh = 2r c G c  .In the two upper graphs (sphere diameters of 26 and 67 nm), the procedure is successful.The contact stiffnesses giving the best match are 0.2 N m −1 and 1 N m −1 for d Sph = 26 and 67 nm, respectively.The agreement between simulation and experiment is best for the smallest loss tangent (tan( L ) = 0.3).These values of the contact stiffness are much below what Mindlin's theory combined with the DMT model would predict.The factor   amounts a few percent, only.The lowered value may be explained with a contact radius being smaller than what is predicted by the DMT model.It may also be explained with nanoscale roughness, impurities, electrostatic interactions, and the role of the ambient liquid.The DMT model (which predicts r c ∝ r Sph 2/3 ) combined with the Mindlin model (which predicts  Sh ∝ r c ) would suggest that the ratio of the contact stiffnesses obtained for particle diameters of 26 and 67 nm would be (67/26) 2/3 = 1.87.Judging from the above analysis, the ratio is 5 in the experiments by the Krakow group.
In order to gain further insight, system parameters (such as the surface charges) need to be varied and the dependence of the contact stiffness on these system parameters needs to be studied.Of course, combined experiments, for instance involving an AFM tip, would also help.As mentioned in the introduction, such measurements might also be exploited on a more practical level, adsorbing a layer of interest to the surface (DNA, brushes, responsive gels, …) and adsorbing spheres onto these layers.On the right, the data are represented in a way, which eliminates coverage.In this form, they can be compared with experiment (dashed lines).The same representation was chosen at the bottom of Figure 1.
Presumably, such a layer would take an influence on the contact stiffness.
For the largest sphere, the simulation does not lead to a clear minimum in  2 .The simulation has difficulties reproducing the large acoustic ratios found in these experiments.The most plau-sible explanation is that the grid resolution Δx = 1/6 d Sph is problematic for this case.A grid resolution of 24 nm results.FreqD-LBM only reproduces the Stokes equation if Δx << . takes values of 145, 112, and 95 nm for the overtones 3, 5, and 7. Presumably, the ratio Δx/ is not small enough for the largest spheres.

Conclusion
The stiffness of the contacts between particles and the surface of a QCM-D can be inferred from the comparison of −Δf/n and ΔΓ/n obtained during an adsorption experiment with equivalent values computed with a FreqD-LBM simulation.The simulation requires the particle surfaces to be configured as oscillating walls.The dynamics of the particles is modeled independently of the FreqD-LBM simulation.For the small particles, the simulation could be made to match the experiment.The derived contact stiffnesses can be interpreted in the frame of contact mechanics.They may also give access to the properties of adsorbed layers, located between the resonator surface and the spheres.

Figure 1 .
Figure 1.Experimental data to be reproduced with the simulation, adapted from Ref. [1].Polystyrene spheres with a positive surface potential were adsorbed to a negatively charged resonator surface.The sphere diameters, d Sph , are indicated at the top.The upper panels show the adsorption kinetics.Decreasing Δf/n indicates increasing coverage, but there is no quantitative relation between coverage and −Δf/n.Because the coverage is not known in absolute terms, the acoustic ratio is plotted versus −Δf/n in the lower panels.These graphs allow for a comparison between experiment and simulation (see Figure3).

Figure 2 .
Figure 2.In the implementation chosen here, the FreqD-LBM simulation is limited to the bulk liquid.The particles are walls with a certain (local) oscillatory velocity, computed in each simulation step from the particle's velocity and its rate of rotation.In the steady state, the force and the torque exerted by the bulk liquid are equal and opposite to the force and the torque exerted by the contact (ignoring a small contribution from inertia).

Figure 3 .
Figure 3.A cut through the simulation volume.Dashed red lines indicate the planes of the other cut.The coverage is 0.4 in this plot.

Figure 4 .
Figure 4.The two columns on the left and the center show the overtone normalized shifts of frequency and bandwidth, Δf/n and ΔΓ/n versus coverage.On the right, the data are represented in a way, which eliminates coverage.In this form, they can be compared with experiment (dashed lines).The same representation was chosen at the bottom of Figure1.

Figure 5 .
Figure 5. Evaluation of the agreement between the simulation and experiment expressed in terms of  2 Δf/n and  2 ratio as defined in the main text.Different colors denote different loss tangents as indicated in the legend.The parameter   from Figure 4 was converted to a contact stiffness with the relation  Sh = 2r c G c   .