Formation of inverse Chladni patterns in liquids at microscale: roles of acoustic radiation and streaming-induced drag forces

While Chladni patterns in air over vibrating plates at macroscale have been well studied, inverse Chladni patterns in water at microscale have recently been reported. The underlying physics for the focusing of microparticles on the vibrating interface, however, is still unclear. In this paper, we present a quantitative three-dimensional study on the acoustophoretic motion of microparticles on a clamped vibrating circular plate in contact with water with emphasis on the roles of acoustic radiation and streaming-induced drag forces. The numerical simulations show good comparisons with experimental observations and basic theory. While we provide clear demonstrations of three-dimensional particle size-dependent microparticle trajectories in vibrating plate systems, we show that acoustic radiation forces are crucial for the formation of inverse Chladni patterns in liquids on both out-of-plane and in-plane microparticle movements. For out-of-plane microparticle acoustophoresis, out-of-plane acoustic radiation forces are the main driving force in the near-field, which prevent out-of-plane acoustic streaming vortices from dragging particles away from the vibrating interface. For in-plane acoustophoresis on the vibrating interface, acoustic streaming is not the only mechanism that carries microparticles to the vibrating antinodes forming inverse Chladni patterns: In-plane acoustic radiation forces could have a greater contribution. To facilitate the design of lab-on-a-chip devices for a wide range of applications, the effects of many key parameters, including the plate radius R and thickness h and the fluid viscosity μ, on the microparticle acoustophoresis are discussed, which show that the threshold in-plane and out-of-plane particle sizes balanced from the acoustic radiation and streaming-induced drag forces scale linearly with R and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt \mu$$\end{document}μ, but inversely with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt h$$\end{document}h. Electronic supplementary material The online version of this article (doi:10.1007/s10404-017-1888-5) contains supplementary material, which is available to authorized users.


Introduction
Arranging particles and cells into desired patterns for labon-a-chip biological applications using ultrasonic fields, i.e. acoustophoresis, by means of bulk and surface acoustic wave techniques, have attracted increasing interest in recent years (Bruus et al. 2011;Friend and Yeo 2011). When an ultrasonic standing/travelling wave is established in a micro-channel containing an aqueous suspension of particles, two main forces act on the particles: the acoustic radiation force and the streaming-induced drag force. In most bulk and surface micro-acoustofluidic manipulation devices, the latter is generally considered to be a disturbance because it places a practical lower limit on the particle size that can be manipulated by the former (Wiklund et al. 2012;Drinkwater 2016). Nevertheless, acoustic streaming flows have been applied to play an active role in the functioning of such systems. (Hammarstrom et al. 2012(Hammarstrom et al. , 2014Yazdi and Ardekani 2012;Antfolk et al. 2014;Devendran et al. 2014;Ohlin et al. 2015;Cheung et al. 2014;Huang et al. 2014;Patel et al. 2014;Destgeer et al. 2016; Abstract While Chladni patterns in air over vibrating plates at macroscale have been well studied, inverse Chladni patterns in water at microscale have recently been reported. The underlying physics for the focusing of microparticles on the vibrating interface, however, is still unclear. In this paper, we present a quantitative three-dimensional study on the acoustophoretic motion of microparticles on a clamped vibrating circular plate in contact with water with emphasis on the roles of acoustic radiation and streaminginduced drag forces. The numerical simulations show good comparisons with experimental observations and basic theory. While we provide clear demonstrations of threedimensional particle size-dependent microparticle trajectories in vibrating plate systems, we show that acoustic radiation forces are crucial for the formation of inverse Chladni patterns in liquids on both out-of-plane and inplane microparticle movements. For out-of-plane microparticle acoustophoresis, out-of-plane acoustic radiation forces are the main driving force in the near-field, which prevent out-of-plane acoustic streaming vortices from dragging particles away from the vibrating interface. For inplane acoustophoresis on the vibrating interface, acoustic streaming is not the only mechanism that carries microparticles to the vibrating antinodes forming inverse Chladni patterns: In-plane acoustic radiation forces could have a Electronic supplementary material The online version of this article (doi:10.1007/s10404-017-1888-5) contains supplementary material, which is available to authorized users.
The ability to use ultrasonic fields for manipulation of particles and fluids has a long history which can date back to many eminent scientists including Chladni (1787), Faraday (1831), Kundt and Lehmann (1874), Rayleigh (1883), King (1934), Gorkov (1962). As early as 1787, the German physicist Chladni (1787) observed that randomly distributed sand particles on a vibrating metal plate could group along the nodal lines forming a wide variety of symmetrical patterns. The various patterns formed at different modes of resonance were called Chladni figures. Chladni also reported that fine particles would move in the opposite direction, to the antinodes, which was further studied by Faraday (1831), who found that it was due to air currents in the vicinity of the plate, i.e. acoustic streaming. The latter phenomenon was revisited by Van Gerner et al. (2010 who showed that it will always occur when the acceleration of the resonating plate is lower than gravity acceleration. Zhou et al. (2016) recently proposed an approach which is able to control the motion of multiple objects simultaneously and independently on a Chladni plate.
Recently, Vuillermet et al. (2016) demonstrated that it is possible to form two-dimensional inverse Chladni patterns on a vibrating circular plate in water at microscale, which extended an earlier work from Dorrestijn et al. (2007), who showed formation of one-dimensional (1D) Chladni patterns on a vibrating cantilever submerged in water, where microparticles and nanoparticles were found to move to the antinodes and nodes of the vibrating interface, respectively. Both works have depicted the two-dimensional streaming field in the near-field and emphasized the effects of in-plane streaming flow on the collections of particles at vibrating antinodes or nodes. Practical manipulation on vibrating plates, however, is three-dimensional (3D) including out-of-plane and in-plane manipulation, and interestingly, in such systems, little work has been done on the impact of acoustic radiation forces, the main engine for particle and cell manipulation in other acoustofluidic manipulation devices. Unlike microparticle acoustophoresis in bulk and surface standing wave devices that have been well studied (Barnkob et al. 2012;Muller et al. 2012Muller et al. , 2013Lei et al. 2014;Hahn et al. 2015;Nama et al. 2015;Oberti et al. 2009), the literature is lacking a quantitative analysis of microparticle acoustophoresis over vibrating plate systems.
In this paper, we will show a detailed 3D study on the main forces for the formation of inverse Chladni patterns on a clamped vibrating circular plate in contact with water (see Fig. 1 for the configuration). Both out-of-plane and inplane microparticle acoustophoresis are discussed and the contributions of main driving forces are compared, which enables a clear presentation of the underlying physics of microparticle manipulation in such systems. The many key parameters, including the plate thickness and radius, the vibration amplitude and the fluid viscosity, on the microparticle acoustophoresis are discussed. We believe that this work could provide an excellent tool on analysing microparticle acoustophoresis in vibrating plate systems and on guiding device designs for the better control of patterning of microparticles at various sizes as well as for single particle and cell manipulation.

Numerical method
We use bold and normal-emphasis fonts to represent vector and scalar quantities, respectively. Here, we assume a homogeneous isotropic fluid, in which the continuity and momentum equations for the fluid motion are.
where ρ is the fluid density, t is time, u is the fluid velocity, p is the pressure and μ and μ b are, respectively, the dynamic and bulk viscosity coefficients of the fluid.
Taking the first and second order into account, we write the perturbation series of fluid density, pressure and velocity: (Bruus 2012) where the subscripts 0, 1 and 2 represent the static (absence of sound), first-order and second-order quantities, respectively. Substituting Eq. (2) into Eq. (1) and considering Repeating the above procedure, considering the equations to the second order and taking the time average of Eq. (1) using Eq. (2), the continuity and momentum equations for solving the second-order time-averaged acoustic streaming velocity can be turned into where the upper bar denotes a time-averaged value and F is the Reynolds stress force (Lighthill 1978). When modelling the steady-state streaming flows in most practical acoustofluidic manipulation devices, the inertial force u 2 · ∇u 2 is generally negligible compared to the viscosity force in such systems, which results in the creeping motion. The divergence-free velocity u M 2 = u 2 + ρ 1 u 1 /ρ 0 , derived from Eq. (4a), is the mass transport velocity of the acoustic streaming, which is generally closer to the velocity of tracer particles in a streaming flow than u 2 (Nyborg 1998).
In this work, only the boundary-driven streaming field was solved because an evanescent wave field is established (see below) such that the overall streaming field is dominated by the boundary-driven streaming. Moreover, as the inner streaming vortices are confined only at the thin viscous boundary layer [thickness of δ v ≈ 0.6 µm at 1 MHz in water (Bruus 2012)], for numerical efficiency, we solved only the 3D outer streaming fields using Nyborg's limiting velocity method (Nyborg 1958;Lee and Wang 1989) as those published previously (Lei et al. 2013(Lei et al. , 2014(Lei et al. , 2016. Although the inner streaming fields were not computed in this work, they can, of course, be known from the limiting velocity field.

Numerical model, results and discussion
To validate the numerical results, a clamped circular plate of radius R = 800 µm and thickness h = 5.9 µm was firstly considered, which has a same size to the one used in Vuillermet et al.'s experiments (2016 Table 1. The model configuration is shown in Fig. 3a, where a cylindrical fluid-channel-only model was considered. Cartesian (x, y, z) and cylindrical (r, θ , z) coordinates were used for the convenience of calculations. The finite element package COMSOL 5.2 (COMSOL Multiphysics 2015) was used to solve all equations. The modelled final particle (radius of 30 µm) positions driven by the main forces including acoustic radiation forces, streaming-induced drag forces and buoyancy forces at two vibrating modes are shown in Fig. 2a. It can be seen that the inverse Chladni patterns the  Vuillermet et al. (2016) Copyrighted by the American Physical Society. The particle properties used in simulations are included in Table 1 50 Page 4 of 15 microparticles form compare well with Vuillermet et al.'s (2016) experimental observations. In the following, we will show step by step why microparticles are gathered to the vibrating antinodes forming inverse Chladni patterns and the contributions of various driving forces on the acoustophoretic motion of microparticles at various sizes. It is noteworthy that we have previously applied a fluid-channel-only model to study the 3D transducerplane streaming fields in bulk acoustofluidic manipulation devices (Lei 2015), where the excitation of transducer was approximated by a Gaussian distribution of boundary vibration. The fluid-channel-only model applied in this  work has more merits because we can easily write down the displacement equation when the circular plate vibrates at a resonant mode (see Eq. (6) below), and thus, there is no need to make an approximation on the boundary vibrations as we did in the previous models (Lei et al. 2013(Lei et al. , 2016.

Resonant frequencies
Resonant frequencies at various modes were firstly modelled, which are shown in Table 2. For comparison, the modelled eigenfrequencies of first eight modes for another two cases, namely no load and load with air, are also presented. It can be seen that the resonant frequencies for vibrations in air and those in vacuum are very close; differences are small enough to be considered as numerical errors, suggesting that omitting the influence of air does not introduce any significant error on the resonant frequencies.
The resonant frequencies of vibration in contact with water, however, have been reduced at least by a factor of 3 for all the modes presented, which means that we have to consider the influence of external load introduced by the surrounding water. All the results shown in this paper are for the (4, 1) mode (δ v ≈ 1.84 µm) unless otherwise stated. The computations were performed on a Lenovo Y50 running Windows 8 (64-bit) equipped with 16 GB RAM and Intel(R) Core(TM) i7-4710HQ processor of clock frequency 2.5 GHz. The mesh constitution was chosen based on the method described in a previous work (Lei et al. 2013), which chooses the mesh size to obtain steady solutions, i.e. ensuring further refining of mesh does not change the solution significantly. This model resulted in 131,521 mesh elements, a peak RAM usage of 4.96 GB (at the acoustic step), and a running time of about 4 h for solving the steps described between Sects. 3.2 and 3.6 below.

First-order acoustic fields
The first-order acoustic fields were modelled using the COMSOL 'Pressure Acoustics, Frequency Domain' interface, which solves the harmonic, linearized acoustic problem, taking the form, where ω is the angular frequency and c is the speed of sound in the fluid. The acoustic fields in the model regime were created by a harmonic vibration of the bottom edge (i.e. the plate) coupled with radiation boundary conditions on all other edges. For comparison, we also tried adding perfect matching layers around the cylindrical domain to absorb all outgoing waves and found that the differences on all the (5) ∇ 2 p 1 + ω 2 c 2 p 1 = 0, modelled quantities between these two methods are within 3%. To give a clear presentation of results, we show here the results modelled form radiation boundary conditions. For a (m, n) vibrating mode, the plate displacement amplitude can be written as where J m (·) is the Bessel function of the first kind of order m and α mn is the nth zero of J m (·). The results presented in this paper were obtained at a vibration amplitude of 0.4 µm unless otherwise stated. The vibration amplitude has a limited effect on the shape of microparticle trajectories as both the acoustic radiation force and streaming-induced drag force scale with the square of the vibration amplitude (more discussions can be found below).
As shown in Fig. 3b, a standing wave field was established on the vibrating interface with acoustic pressure nodes and antinodes locating at plate displacement nodes and antinodes, respectively. The standing wave field is shown more clearly in Fig. 3d, where the in-plane circumferential acoustic pressure magnitudes are plotted. The out-of-plane acoustic pressure magnitudes over a vibrating antinode are plotted in Fig. 3c, which shows that the acoustic pressure magnitudes decay exponentially with the increase in distance from the vibrating interface. The reason is that the plate wave travels at the vibrating interface at a subsonic regime leading to an evanescent wave field: the plate wave velocity at substrate surface u = f r ≈ 55 m/s ≪ u l , where λ is the acoustic wavelength, f r is the resonant frequency and u l is the speed of sound in the liquid.

Acoustic radiation forces
The corresponding 3D acoustic radiation forces were solved from the Gorkov equation (Gorkov 1962), where E kin and E pot are the time-averaged kinematic and potential energy, ρ p and ρ f are, respectively, the density of particle and fluid, β p = 1/ ρ p c 2 p and β f = 1/ ρ f c 2 f are the compressibility of particle and fluid, and V 0 is the particle volume (see Table 1 for model properties). Equation (7) is valid for particles that are small compared to the acoustic wavelength λ in the limit r 0 / ≪ 1 (where r 0 is the radius of the particle) in an inviscid fluid in an arbitrary sound field. (Gorkov 1962) When a particle moves close to the vibrating plate, the acoustic radiation forces may oscillate weakly with a decrease in distance to the plate due to the multiple-scattering interaction and wall interference, while the force magnitudes will not be significantly affected (Wang and Dual 2012).
The modelled acoustic radiation force fields are shown in Fig. 4. As shown in Fig. 4c, the out-of-plane acoustic radiation forces also decrease exponentially with the increase in distance from the vibrating interface. In the near-field, at this vibrating amplitude, the out-of-plane acoustic radiation forces have a greater contribution on the sedimentation of microparticles than the buoyancy forces. With an increase in vibration amplitude, we can expect dominant out-of-plane acoustic radiation forces over buoyancy forces. Interestingly, as shown in Fig. 4b, the in-plane acoustic radiation forces carry microparticles away from the acoustic pressure nodes and converge at antinodes from all directions, in contrast with the conditions usually found in bulk and surface standing wave manipulation devices, where the acoustic radiation forces move most particles and cells of interest to the acoustic pressure nodes (Glynne-Jones et al. 2012). Examining Eq. (7), it can be seen that the acoustic radiation force is a gradient of the force potential, which contains a positive contribution from the kinematic energy (weighted by a function of the fluid and particle densities) and a negative contribution from the potential energy (weighted by a function of the fluid and particle compressibility). Comparing the contributions of these two terms in this model, it was found that the kinetic energy term dominates in the force potential, which drives microparticles to the vibrating antinodes.

Acoustic streaming fields
The 3D acoustic streaming field was modelled using Nyborg's limiting velocity method (Nyborg 1958;Lee and Wang 1989). It was shown that if the boundary has a radius of curvature that is much larger than the acoustic boundary layer, then the time-averaged velocity at the extremity of the inner streaming (the 'limiting velocity') can be approximated as a function of the local, first-order linear acoustic field. The outer streaming in the bulk of the fluid can then be predicted by a fluidic model that takes the limiting velocity as a boundary condition. The applicability and viability of the limiting velocity method have been further discussed recently (Lei et al. 2017). In Cartesian coordinates, the limiting velocity field at the driving boundaries (z = 0) can be written as (8a) where u L and v L are the x-and y-components of the limiting velocity field, u 1 , v 1 and w 1 are the x-, y-and z-components of the acoustic velocity vector u 1 , Re{·} denotes the real part of a complex value and * is the complex conjugate. A COMSOL 'Creeping Flow' interface was used to model the acoustic streaming field, which solves As only outer streaming fields are solved in this method, with the assumption of low velocity and incompressible flow, the first term in the left-hand side of Eq. (4a) is zero and thus u 2 = u M 2 (Hamilton et al. 2003). Then, as discussed by Lighthill (1978), the Reynolds stress in the bulk of the fluid can set up hydrostatic stresses, but in the absence of attenuation these will not create vortices, hence these terms are not included in Eq. (9b). The 3D outer acoustic streaming fields in the considered model regime were generated by the limiting velocity field on the vibrating interface (see Fig. 5a) along with no-slip boundary conditions (u 2 = 0) on all other edges.
The limiting velocity field and the 3D acoustic streaming fields are shown in Fig. 5. It can be seen that, similar to the distribution of in-plane acoustic radiation forces, the limiting velocities (i.e. the in-plane acoustic streaming velocity field) converge at the acoustic pressure antinodes from all directions leading to acoustic streaming vortices on out-of-planes perpendicular to the vibrating interface as those plotted in Fig. 5b, c, where, in order to give a clear demonstration of the 3D acoustic streaming fields, only the acoustic streaming vortices above one acoustic pressure antinode are plotted.

Acoustic streaming-induced drag forces
Based on the acoustic streaming velocity field, we can calculate the acoustic streaming-induced drag forces on microparticles from the stokes drag, where v is the particle velocity. Equation (10) is valid for particles sufficiently far from the channel walls (Happel 1965). Since microparticle acoustophoresis discussed in this work is closely associated with the vibrating plate, it is necessary to take into account the wall effect on the streaming-induced drag forces when a particle moves close to the bottom wall. When a sphere particle moves perpendicularly towards or in parallel to the vibrating plate, the streaming-induced drag force should be corrected by multiplying a wall-effect-correction factor χ or γ, respectively, which can be expressed as (Happel 1965) where H is the distance from the centre of the particle to the plate and A = 9/16, B = 1/8, C = 45/256 and D = 1/16. The 3D acoustic streaming-induced drag forces are shown in Fig. 6, where, for comparison, the buoyancy forces are also plotted. As shown in Fig. 6c, with the increase in distance from the vibrating interface, the out-of-plane streaming-induced drag forces rise rapidly to the maximum value in the near-field and then fall gradually to zero in the far-field. The wall effect can increase the maximum our-of-plane streaming-induced drag force by approximately a factor of 2 in this model. Also, it can be seen that, for a small vibration amplitude of w = 0.4 µm, the maximum out-of-plane streaminginduced drag force is larger than the buoyancy force on a particle with a radius of 30 µm. With an increase in vibration amplitude, we can expect even larger acoustic streaming-induced drag forces while the buoyancy forces remain the same. Therefore, it might be reasonable to say that introducing only the streaming effects is not enough to explain the sedimentation of microparticles, especially for those with r 0 < 30 µm, where the differences between the out-of-plane streaming-induced drag forces and the buoyancy forces are even larger, as plotted in Fig. 6d, because the former and the latter scale with the particle radius and particle volume, respectively.

Microparticle trajectories
From the acoustic radiation forces and streaming-induced drag forces that have been calculated, together with the buoyancy forces, microparticle (polystyrene beads) trajectories were modelled, following where m p is the particle mass, F B is the buoyancy, F G is the gravity and g is the gravity acceleration. In this work, Fig. 6 (Colour online) a 3D streaming-induced drag forces on a particle with a radius of 30 µm (|F d |, N); b in-plane |F d |; c out-of-plane |F d | [red arrow in (a)]; and d comparisons of maximum out-of-plane |F d | [peak in (c)] with the buoyancy forces for various particle sizes (radius of r 0 ). The inset in (c) shows the directions of the plotted forces above a vibrating antinode. F B and F G are the buoyancy and gravity, respectively it is assumed that all the forces, including acoustic radiation, streaming-induced drag and buoyancy forces, act on the centre of spherical particles (otherwise, integration of forces over the particle surface would be needed when the particles are close to the boundaries). It is noteworthy that, in addition to these main driving forces, a particle-particle interaction force was used in this model. The particle-particle interaction force can be expressed as where k s is the spring constant, r i is the position vector of the ith particle, and r e is the equilibrium position between particles. In this model, k s = 2.5 × 10 −4 N/m for polystyrene beads (Jensenius and Zocchi 1997) and r e was set as 2r 0 to avoid all particles being concentrated to a single point.
Here, a COMSOL 'Particle Tracing for Fluid Flow' interface was used to solve Eq. (12) to model the particle trajectories. The shape of the trajectories is independent of the pressure amplitude since both the acoustic radiation forces and steaming-induced drag forces scale with the square of pressure; results are presented here for an excitation amplitude of w = 0.4 µm. An array of tracer particles (given the properties of polystyrene beads of radius 30 µm) are seeded at t = 0. Acoustic radiation forces, streaminginduced drag forces and buoyancy forces act on the particles, resulting in the motion shown in Fig. 7. It can be seen that, in the considered model regime, particles with a radius of 30 µm first move towards the vibrating interface driven by the predominant out-of-plane forces and are then carried to their closest acoustic pressure antinodes by in-plane forces, resulting in spider-like trajectories and inverse Chladni patterns on the vibrating interface within seconds. Generally, particles closer to the vibrating interface take less time to settle for stronger driving forces. Particles with smaller sizes take longer to locate at the acoustic pressure antinodes for smaller driving forces and will follow the outof-plane streaming vortices leading to acoustic streamingdominated trajectories close to those shown in Fig. 5b, c while r 0 < 6.9 µm (see explanations below and videos in the Supplemental material).
Out-of-plane acoustophoresis. A single particle outof-plane acoustophoresis is directly acted upon by the acoustic radiation force, the buoyancy force and the acoustic streaming-induced drag force. The equation of motion for a spherical particle of out-of-plane velocity v out above an acoustic pressure antinode is then As we have seen above, particles are concentrated at the acoustic pressure antinodes, so we take here a particle staying above an acoustic pressure antinode to analyse the contributions of the many forces on the microparticle out-of-plane acoustophoresis. As shown in the inset of Fig. 8a, the streaming-induced drag forces, F out d , competes with other forces above an acoustic pressure antinode as the acoustic streaming flow drives particles away from the pressure antinode, while other forces bring particles to the pressure antinode. Based on the fact that there should be a threshold out-of-plane particle size, r out 0 : for r 0 > r out 0 , particles can be easily concentred to the (15) F out d ∝ r 0 and F out ac + F B + F G ∝ r 3 0 , acoustic pressure antinodes; while for r 0 < r out 0 , particles will follow the out-of-plane acoustic streaming vortices. We define the threshold particle radius r out 0 for crossover from these out-of-plane forces. The out-of-plane forces on particles at various sizes are plotted in Fig. 8a, which shows that, at a small vibration amplitude of w = 0.4 µm, the threshold particle size r out 0 ≈ 6.9 µm. Considering the walleffect-correction for the streaming-induced drag forces, r out 0 ≈ 9.1 µm. This threshold out-of-plane particle size may slightly vary with the vibration amplitude as F B + F G are independent of the vibration amplitude, while F out d and F out ac scale with the square of the vibration amplitude. As shown in Fig. 4c, the buoyancy force is approximately 1/20 of F out ac at w = 0.4 µm on the vibrating interface. With an increase in vibration amplitude, the contribution of buoyancy force will be even smaller on the microparticle acoustophoresis in the near-field. To calculate the limit value of r out 0 , we can set by ignoring the buoyancy forces, which gives Considering the wall-effect-correction for the streaminginduced drag forces, the limit value of r out 0 ≈ 9.4 µm.
(16) F out ac + F out d = 0 (17) In-plane microparticle acoustophoresis. For the inplane microparticle acoustophoresis, it is acted upon by the acoustic radiation force and the streaming-induced drag force. Similar to the analyses above, the equation of motion for a spherical particle of in-plane velocity v in is then As shown in Figs. 4b and 6b, both the in-plane acoustic radiation force, F in ac , and the streaming-induced drag force, F in d , move microparticles to the acoustic pressure antinodes (see also the inset in Fig. 8b). To evaluate the contributions of these two forces on the in-plane microparticle acoustophoresis, we compare their average values over the plate interface because considering the maximum force only may not be accurate. Since both of these in-plane forces point to the acoustic pressure antinodes, they jointly contribute to the focusing of microparticles to the acoustic pressure antinodes provided that the particle sizes are big enough to avoid being driven away from the vibrating interface by out-of-plane acoustic streaming vortices (as discussed in the previous step), which could provide evidence for the much larger particle velocities measured in experiments when compared with the predicted streaming velocities as shown in Vuillermet et al. (2016) work.
Although there is no threshold in-plane particle size for the reason that both the in-plane acoustic radiation force and streaming-induced drag force drive microparticles to the acoustic pressure antinodes, we can figure out the contribution of each force on the in-plane microparticle acoustophoresis. Again, based on the fact that we can expect a critical in-plane particle size, r in 0 : for r 0 > r in 0 , F in ac contribute more to the in-plane acoustophoresis; while for r 0 < r out 0 , F in d have a higher contribution. The value of r in 0 can be found from setting F in ac = F in d , which gives Considering the wall-effect-correction for the streaminginduced drag forces, r in 0 ≈ 27.6 µm. The in-plane forces on particles at various sizes are plotted in Fig. 8b. It is noteworthy that, different to the situation for r out 0 , r in 0 is independent of the vibration amplitude w because both F in ac and F in d scale with the square of w. Actually, it can be seen from Eqs. (17) and (20) that, ignoring the small effect of buoyancy forces in the nearfield, the relationships between the in-plane and outof-plane threshold particle sizes and the ratios of the corresponding streaming-induced drag force and acoustic radiation force are

Effects of key parameters on microparticle acoustophoresis
Having demonstrated the acoustophoresis of microparticles at various sizes for a particular plate (thickness of 5.9 µm and radius of 0.8 mm), in this section, we investigate the effects of many key parameters, including the plate radius and thickness and the fluid viscosity, on the performance of microparticle acoustophoresis in order to facilitate device design for a wide range of applications.
Effects of fluid viscosity. On the one hand, it can be seen from Eq. (8) that the magnitudes of limiting velocities (i.e. the strength of the outer streaming velocities) are independent of the fluid viscosity even though viscosity is the initial cause of acoustic streaming flows. Thus, with a change in fluid viscosity, the streaming-induced drag force, F d , scales linearly with µ, while F ac will remain the same. From Eq. (21), the following relationships are established, Therefore, to eliminate the 'side effect' of streaming flows on the microparticle manipulation, and we can conclude that lowering the fluid viscosity is a viable way to augment the weight of acoustic radiation force on microparticle acoustophoresis.
Effects of plate thickness and radius. To investigate the effects of plate thickness (h) and radius (R) on the microparticle acoustophoresis, we considered a series of h and R ranging from 2 to 14 µm and 0.3 to 1.4 mm, respectively. When one parameter was studied, the other parameter was kept the same. For each case, following the whole numerical procedure described in the sections above, we calculated the threshold in-plane and out-of-plane particle sizes, which are shown in Fig. 9. It can be seen that these two threshold particle sizes have similar variation tendencies: they grow with the increase in R and fall with the rise of h. Compare with basic theory. Turning to the theoretical aspect, as seen from Eq. (21), to determine how these two threshold particle sizes change with the many key parameters, we only need to figure out how the force ratio on the right-hand side varies with these parameters. If we define v rad as the contribution of acoustic radiation force on the particle velocity, considering Eqs. (15) and (19), we have Examining the acoustic field in the near-field, it can be seen from Fig. 3b that, if expanded in the radial direction, the acoustic pressure field (as plotted in Fig. 3d) can be approximated to a 1D standing wave on all circumferences for 0 < r ≪ R, in which the right-hand side of Eq. (23) has the following relation (Barnkob et al. 2012) where Φ ≈ 0.1685 in this work is the acoustic contrast factor and the thermoviscous effects are not included.
For a clamped circular plate with radius of R and thickness of h, the angular frequency for an unloaded case for each (m, n) mode follows (Leissa 1993) Φρ f ωr 2 0 , Fig. 9 Effects of plate radius on the threshold a in-plane particle sizes, r in 0 , and b outof-plane particle sizes, r out 0 (with wall effect) . For (a, b), the plate thickness is the same, h = 5.9 µm. For (c, d), the plate radius is the same, R = 0.8 mm 1 3 where E is the plate Young's modulus, ρ is the plate density and υ is the plate Poisson's ratio. Considering the surrounding water, for a given (m, n) mode, the angular frequency is reduced to where Γ mn is the non-dimensional added virtual mass incremental (NAVMI) factor, values of which can be found in Ref. (Amabili and Kwak 1996), Table 5, in the case of a clamped plate.
Combining Eqs. (21), (23), (24) and (26), the relationships between the threshold in-plane particle sizes and the many key parameters in a 1D standing wave field can be expressed as The calculated values of r in 0 using Eq. (27) and those obtained from our model for various R and h are shown in Fig. 10. It can be seen that the modelled r in 0 compare reasonably well with the calculated values under the 1D standing wave approximation. The differences between the calculated values and those modelled may be attributed to the reason that, compared to an approximated 1D standing wave, the acoustic field in the near-field is a more complex pattern. However, due to the complexity of the problem, the good comparisons between our model and the calculated values indicate that the approximated 1D standing wave may have captured the main features of (4, 1) mode and our model can be applied to study the basic physics of microparticle acoustophoresis on vibrating plate systems for even more complex vibrating modes.

Mode switching
Eigenfrequency studies show that two orthogonal vibrating patterns for each (m, n) vibrating mode could be excited at two adjacent frequencies (typically hundreds of Hz difference) provided that the modes are high enough (m ≥ 1).
As shown in Fig. 11, the phase angle between two adjacent acoustic pressure antinodes of these two orthogonal patterns is For this specific model, both the in-plane acoustic radiation force and streaming-induced drag force diverge from the vibrating nodes and converge at the vibrating antinodes, so when switching from one mode (e.g. mode 1 in Fig. 11) to the other orthogonal mode (e.g. mode 2 in Fig. 11) a particle tends to move from the vibrating antinode of the former to its closest antinode of the latter either clockwise or anticlockwise depending on the initial position of the particle (assuming the initial position of the particle slightly shifts from the vibrating antinode). The potential underlying mechanism for the circular manipulation of a single particle is schematically shown in Fig. 11. It can be seen that, for each mode switching, the particle can move by an angle, θ = π/2m, while its distance to the centre of the circular membrane will remain the same. To complete a full circle of (28) θ = π 2m .  2010) who showed that beads can be brought to any arbitrary point between the half and quarter-wave nodes when rapidly switching back and forth between half and quarter wavelength frequencies in bulk acoustofluidic devices.

Conclusions
We have investigated the 3D acoustophoretic motion of microparticles due to acoustic radiation, acoustic streaming, gravity and buoyancy over a clamped vibrating circular plate in contact with water. The underlying physics of microparticle acoustophoresis over vibrating plates has been studied in detail. Previous predominant analyses have emphasized the in-plane acoustic streaming flows on the formation of inverse Chladni patterns, which, according to this study, may not be complete. For in-plane microparticle acoustophoresis, both the in-plane acoustic radiation forces and the in-plane streaming-induced drag forces were shown to drive microparticles to their closest vibrating antinodes. For out-of-plane microparticle acoustophoresis above vibrating antinodes, in addition to the buoyancy forces, one has to consider the acoustic radiation forces in the near-field, which prevent the out-of-plane streaming vortices from dragging microparticles away from the vibrating interface. Based on the high efficiency of this numerical model, the threshold in-plane and out-of-plane particle sizes balanced from the acoustic radiation and streaminginduced drag force under all vibrating modes can be readily obtained. An important next step is to achieve a direct experimental verification of numerical modelling. Given a successful experimental verification, this 3D model could be extended to include the thermoviscous effects (Muller and Bruus 2014) to obtain more accurate results, but it would be very computationally expensive. According to a study by Rednikov and Sadhal (2011), the thermoviscous effects can increase the streaming velocities by 18% for water at 20 °C which, thus, will shift the threshold particle sizes.
The good comparisons between our modelling and experiments and basic theories indicate that our numerical model could be used together with high-precision experiments as a better research tool to study the many yet unsolved problems. For example, modelling suggests that mode switching between two adjacent frequencies may be used for circumferential manipulation of a single particle or a pair of particles, which might provide routes for the study of particle-particle and particle-wall interactions in acoustofluidics.
While we have shown here 3D particle size-dependent acoustophoresis over an ultrathin circular plate in water, we believe that this strategy could be applied to analyse 3D acoustophoretic motion of microparticles in other vibrating plate systems regardless of fluid medium and thickness, shape and material of plates. One particular application would be acoustophoretic handling of sub-micrometre particles, such as small cells, bacteria and viruses, whose movements are usually dominated by acoustic streaming flows. From the modelled results and the general scaling law given in Eq. (27), we can conclude that increasing plate thickness, decreasing the plate diameter and lowering the viscosity of the liquid are probably the most viable way to conduct such manipulation.
The above-mentioned applications demonstrate that our numerical model is timely and has a huge potential on studies of basic physical aspects of microparticle acoustophoresis in vibrating plate systems and the design of lab-on-achip devices. Fig. 11 (Colour online) A schematic representation of the underlying mechanism for the circular manipulation of a single particle by continuous mode switching between two (m, n) orthogonal modes. To complete a full circle of movement (i.e. θ = π/2m), 4 m times of mode switching are required