Proof of principle study of ultrasonic particle manipulation by a circular array device.

A feasibility study of a circular ultrasonic array device for acoustic particle manipulation is presented. A general approach based on Green's function is developed to analyse the underlying properties of a circular acoustic array. It allows the size of a controllable device area as a function of the number of array elements to be established and the array excitation required to produce a desired field distribution to be determined. A set of quantitative parameters characterizing the complexity of the pressure landscape is suggested, and relation to the number of array elements is found. Next, a finite-element model of a physically realizable circular piezo-acoustic array device is employed to demonstrate that the trapping capability can be achieved in practice.

Introduction

The acoustic radiation force (King 1934) is the basis of ultrasonic methods of particle manipulation. The technique is complementary to optical and dielectric methods, in particular for larger biological cells (Hultström ; Bazou ) where optical methods reach their limits, and for particles that do not posses high optical contrast in a medium. Owing to the typical micrometre scale of the fields and nanoNewton force magnitudes, ultrasonic particle manipulation techniques are used in applications where large numbers of particles need to be handled, as in particle fractionation (Yasuda ; Ratier & Hoyos 2010), sorting (Johnson & Feke 1995) and medium exchange (Augustsson ) devices or for cell agglomeration that allows filtering of cells from fluids (Coakley 1997; Coakley ). Also, the acoustic radiation force was reported to be used in the production of meta-materials (Saito et al. 1998, 1999; Mitri ).

Of high interest are more general-purpose acoustic devices that allow particles to be manipulated around the acoustic cavity rather than driven to specified locations determined by a stationary acoustic field. A review in Courtney indicates four types of acoustic particle manipulation approaches. These are: mode switching, ‘acoustical tweezers’, linear arrays and counter-propagating waves.

The mode switching manipulation method is based on trapping the particles at the nodes of resonant standing waves and spatially shifting the nodes by means of switching between resonant frequencies (Trinh ; Min ; Glynne-Jones ). The acoustic tweezers method is similar in principle to the optical tweezer technique (Ashkin ) and operates by trapping the particles in the focal point of an acoustic beam (Wu & Du 1990; Wu 1991; Lee & Shung 2006; Lee ). Particle manipulation in this method is achieved by physical displacement of the transducers generating the field. A development of this approach based on manipulation by propagating Bessel beams was analysed in the work of Marston (2006) and Mitri (2008). In the method based on linear arrays, two-dimensional traps are formed by activating selected elements of the array opposing a passive reflector (Kozuka ; Démoré ). Here, the manipulation is achieved by selective switching of the transducers. The last approach involves generating a standing wave as a sum of two counter-propagating waves. This allows the locations of the field nodes to be changed by varying the relative phase between the two sources. The problem of generating counter-propagating standing waves without generating resonant standing waves (Kwiatkowski & Marston 1998) was resolved by using transducers acoustically matched to the fluid with an absorbing backing to prevent reflections leading to resonant modes (Courtney ). This method allows a manipulation of particles on a two-dimensional plane to be achieved using a pair of orthogonal phase-controlled counter-propagating waves (Courtney ).

The last approach allows multiple particle traps located at the nodes of a two-dimensional lattice to be moved en masse but does not allow the independent control of individual traps. Greater flexibility is therefore required to match the particle manipulation dexterity achieved using the optical tweezers (Ashkin ). Thus, the principles of the device design reported by Courtney are taken a step forward to produce a dextrous acoustic manipulation system. It is based on a circular ultrasonic array such as shown in figure 1. This conceptual device consists of a piezo-ceramic ring divided into a number of identical separated elements, an absorbing backing layer, a matching layer and a fluid chamber. The matching layer is used to optimize the coupling of the acoustic signal into the fluid chamber in which the particle trapping and manipulation take place. To generate an acoustic radiation force landscape f (King 1934) related to acoustic force potential U (Gorkov 1962) as f=−∇U, each transducer in the array is individually excited through a dedicated control channel. The device excitation is designed to create a force potential landscape with wells in which particles can be trapped, where the force potential U on a spherical non-elastic particle is given by where p is the total pressure field, a is the radius of the particle, ρ is the medium density, c is the speed of sound and the coefficients f1 and f2 are defined as where ρ0 and c0 are the density and the speed of sound in the particle, respectively (Gorkov 1962). Because the force potential field in equation (1.1) created by a pressure distribution J(r), where J(r) is a Bessel function of the first kind of order α, produces a potential with a local minimum at r=0 for orders α>0 (Marston 2006; Mitri 2008), such pressure distribution is capable of trapping particles with densities ρ0>ρ and would be a natural choice of the physically viable field in a circular geometry characteristic of the device considered here. For particles with small densities ρ0<ρ, a J0(r) field similar to the optical tweezers (Ashkin ) should be considered.

Figure 1.

Schematic of a 16 element circular piezo-electric transducer array device. Each element in the array is excited by an individual voltage source. Matching between the transducers and the fluid is achieved using a quarter-wavelength matching layer and the attenuating backing is used to absorb the waves reflected from the back face of the transducers. (Online version in colour.)

The purpose of the analysis presented in the manuscript was to demonstrate the feasibility of the circular ultrasonic array device to form and independently manipulate multiple acoustic radiation force traps for dextrous particle manipulation.

The analysis of the device is divided into two parts. In the first part, fundamental properties of a generic circular acoustic array device are examined. Here, the range of acoustic traps obtainable in the device is found, the dependance between the fraction of the device area in which control over the acoustic field can be achieved and the number of array elements is established, and the relation between the field at the boundary and the field inside the area enclosed by the boundary is determined. This relation allows the required excitation of the array elements in a practical device to be calculated and the complexity of the radiation force landscape to be characterized by a set of quantitative parameters. The maximum complexity attainable for a given number of array elements is found.

In the next part of the analysis, these general properties are applied to the more realistic device shown in figure 1. The major difference from the first part of the analysis is the model of the source of acoustic radiation at the device boundary. Whereas in the first part of the analysis, the boundary is formed by an array of non-reflecting matched sources, in the realistic device, the piezo-elements forming the boundary are mutually reflective because the transducer-to-fluid matching is not perfect. Therefore, finite-element (FE) modelling is required in the latter case to establish the acoustic field generated by an excited array element. This field, together with the relation between the fields on the boundary and in the area enclosed established in the previous part, is used to determine the excitation of the array that produces the required trapping field.

General formalism

Green's function approach

Control over the field distribution in the fluid is achieved by imposing the field on the fluid boundary using a transducer array. The required field on the boundary can be found using a reciprocity relation by starting with a physically viable field in the volume enclosed, and tracing it back to the boundary. However, the reciprocity cannot be implied straightforwardly for a discontinuous and open boundary. In the following discussion, an analysis is presented to establish the relation between the boundary and volume fields in such cases.

In general, the relation of the fields on the boundary and in the fluid volume enclosed can be found using an integral representation derived from Green's theorem for steady-state harmonic waves (Morse & Feshbach 1953, p. 803). For a fluid volume with no body forces, the pressure at an internal point P∈V can be written as (Morse & Feshbach 1953, p. 806) where G(r1,r0) is Green's function and n(r0) is the coordinate in the direction of the normal on a boundary S. This integral representation is not directly applicable because to calculate the pressure at an internal point, both the pressure p and its normal derivative on the boundary must be known. However, taking G(r1,r0) as a combination of fundamental solution of the Helmholtz equation due to a point source −δ(r1−r0) inside the volume V and a homogeneous solution, one of the terms in equation (2.1) can be dropped (Schmerr 1998, pp. 157–159).

The circular array device is described in two-dimensional plane by the coordinate system shown in figure 2, where r0=(r0,θ0) and r1=(r1,θ1) are defined in accordance with equation (2.1) as polar coordinates of the internal point P and the point source relative to the common origin O1, whereas r2=(r2,θ2) is the polar coordinate of an internal point P with respect to the origin O2 located on boundary S. The boundary S in equation (2.1) is defined by a circle of a fixed radius r0=R and the corresponding volume V enclosed by this circular boundary is an infinitely long cylinder of radius R. The fundamental solution of the Helmholtz equation in cylindrical coordinates due to a harmonic line source of wavelength λ located at r0 (figure 2) can be written using Graff's addition theorem as (Martin 2006, pp. 29–61) follows: where the wavenumber k=2π/λ. Green's function satisfying the boundary condition G(r1,r0)=0 is With this choice of Green's function, the first term in equation (2.1) vanishes and the pressure inside the volume p(r1) is uniquely determined by the pressure on the boundary p(R,θ0),

Figure 2.

Definition of coordinate systems used in the integral representation in equation (2.1).

The relation of the boundary and volume pressures implied by equation (2.4) is just an expression of the uniqueness theorem for the Helmholtz equation, stating that for a given boundary condition, the solution inside the volume enclosed is uniquely determined. Practically, this implies that by controlling the pressure distribution on a boundary, any pressure distribution satisfying the Helmholtz equation can be generated inside the volume enclosed. In the case of single-frequency devices, the acoustic landscape is constructed from linear combinations of Bessel functions of the first kind, which are the physically viable non-divergent solutions of the Helmholtz equation. Bessel function traps p(r)∝J(r) ei, of orders α>0 will be considered in the following because only such traps generate a trapping field at the centre r=0 for particles that are denser than the fluid (Marston 2006; Mitri 2008).

A Bessel pressure distribution of order α centred at a point r0 with r0=RT and θ0=θT (figure 3) can be expressed in terms of r1 using Graff's addition theorem as where βT=θT−π and p0 is the pressure amplitude. The corresponding pressure on the boundary is The fact that this boundary pressure gives exactly the pressure distribution as in equation (2.5) can be readily verified by substituting the last expression into the integral in equation (2.4) and using A practically important case of a discontinuous boundary will be considered next.

Figure 3.

Definition of coordinate systems used in the Graff's transformation equation (2.5) with a prototypical J1(kr2) trap located at (RT,θT).

The effects of a discontinuous boundary and the number of array elements

In the circular array device analysed here, the fluid boundary is controlled by a periodic array of separate piezo-elements. The effect of replacing the continuous circular boundary by a finite number N of sources equally distributed along the boundary can be qualitatively understood using the Nyquist theorem. Assuming a linear array and substituting the spacing of D=2πR/N between the elements into the Nyquist condition 1/D>2/λ gives the minimum number of elements to achieve a non-aliased signal Nmin=4πR/λ. If this requirement is violated, aliasing will occur at large angles relative to the normal to the inner surface of the boundary, i.e. in a ring adjacent to the inner boundary surface. To evaluate the effect of replacing the continuous boundary by a finite number N of sources in a circular array quantitatively, we start by representing the boundary by a piecewise constant function of constant pressure on segments of Δφ<Δ length, where Δ=2π/N, and with gaps of zero pressure in between. The integral in equation (2.7) is then replaced by the sum which, using the Poisson summation formula, can be expressed as and the pressure in the volume becomes To further simplify the analysis, we consider the case of infinitely small sources such that . The practical validity of this assumption is demonstrated towards the end of this section.

The sum in the last equation can be broken into , where the s=0 part gives the principal term pT(r1)=p0J(kr2) ei and the remainder contributes to the artefact field created owing to the boundary discontinuity. To determine the artefact field we consider the asymptotic behaviour of the Bessel function for large orders (Abramowitz & Stegun 1964, p. 365), from which one can conclude that while the first term in equation (2.11), J(kR)/J(kR), is limited for m<πeR/λ, the second term, J(kRT), becomes zero for |m+sN−α|>πe(RT/λ). Therefore, if there exist m such that J(kR)/J(kR) is limited ∀m: |m+sN−α|<πe(RT/λ), the artefact field can be estimated by the sum over a finite number of orders of m, where ΔM=πeRT/λ and A are the appropriate coefficients in equation (2.11).

When RT increases, the largest order for non-zero J(kRT) (e.g. N+α+πe(RT/λ) for s=−1) becomes greater than πeR/λ. In this case, only the lower order terms in m can be neglected such that in general the artefact field can be approximated by

In order to understand the structure of equation (2.14), consider first the RT=0 case, in which the only contribution comes from m=±N+α terms. In this case, which, with the help of the asymptotic behaviour equation (2.12), indicates that the artefact field introduces a disturbance outside a circle of radius r1=λ(N−α)/πe. The disturbance field is determined by J(kr1) and goes to zero for r1<λ(N−α)/πe. This case is depicted in figure 4, where p(r1) calculated using equation (2.7) for α=1, N=60 and a set of different values of RT is shown. Because ΔM increases with RT, the lowest order m contributing to the artefact field decreases as implied by equation (2.14). Thus, the disturbance-free region is reduced to r<λ(N−α)/πe−RT, as shown in figure 4.

Figure 4.

A J1 Bessel trap generated by a 60 element point-source array with R=10λ. In two of the cases shown, the trap centre distance from the centre of the array is RT<1/2λ(N/πe)∼3.5λ and in the other two cases, this condition is violated.

The maximum value of RT can be found by demanding that the trap lies inside the distortion-free region. This requirement implies that the useful area of the fluid chamber is defined by In the example shown in figure 4, the maximum value of RT corresponding to N=60 and R=10λ is RT∼3.5λ. The plots of p(r1), calculated according to equation (2.4), show that indeed, for RT<3.5λ, the trap field appears undistorted; for RT=4λ, the trap starts to overlap with the artefact field, and for RT=6λ, the trap is completely destroyed in the artefact field. This set of plots also illustrates how the undistorted region decreases with increasing RT.

Reversing the argument used to obtain the condition on the maximum useful region in equation (2.14), the minimum number of elements such that the useful region is at least as large as the inner radius of the fluid chamber R can be found to be This is a more stringent condition than the well-known Nyquist sampling condition N>4πR/λ on account of the cylindrical geometry of the device.

Up to this point in the analysis, the effect of the finite source size Δφ>0 has been neglected. This assumption is justified if sinc(mΔφ/2) varies slowly for . From here, the maximum size of the source can be estimated as Δl∼πR/20, which is readily satisfied if the number of sources is N>40.

This section can be summarized by indicating that as a result of a boundary discontinuity, an artefact field appears in addition to the principal trapping field. Because the artefact field is negligibly small inside the circular area defined by equation (2.14), it affects the principle field only outside this circle, whereas the field inside remains virtually undisturbed. This means that potentially any pressure field satisfying the Helmholtz equation can be generated by an array of N elements within the area defined by equation (2.14), and the finite number of array elements limits only the spatial extent of the controllable area but not the complexity of the pressure landscape. The complexity of the landscape is determined by the interaction of pressure fields of individual traps, which, in the case of a single operating frequency device considered here, are given by Bessel functions.

Multiple trapping

Selective manipulation of multiple particles is an essential capability of the device analysed here for which an ability to generate complex force landscapes with separated particle traps is required. The trapping complexity grows with the number of traps Np introduced; thus, given the number of traps Np, the minimum average spacing between the traps dmin is established to characterize the configuration complexity. For example, it is clear that for two particles, the minimum spacing will be λ, because at closer distance, two separate traps will blur into one.

The average spacing between the particles is defined using the area S occupied by the particles as where S is the area of the smallest rectangle enclosing all the particles, and dmin is defined as the smallest spacing d between the traps for which the force potential U contains Np individual and separate potential wells that allow Np particles to be trapped. This loose definition of dmin is better explained graphically. In figure 5, a potential landscape U given for configurations of Np=6 traps with different average spacings d is shown. Clearly, for d=1.5λ, the trapping ability of the field is absent but appears for the spacing of d=2.0λ; in that case, six separate closed traps become apparent. Thus, dmin=2.0λ for Np=6 configuration.

Figure 5.

Force potential U for a polystyrene bead with six traps generated by (b) a 60 element point-source array and (a) a continuous boundary. Potential landscapes with different average trap spacing d are shown.

The parameter dmin is independent of the number of array elements N. This can be further seen from figure 5. The plots shown in figure 5a correspond to an ideal case with a continuous boundary, and in figure 5b, the same configurations are obtained with a 60 element point-source array. One can see that in the case of an array with a finite number of elements, the field inside the circle defined by equation (2.14) is identical to the field generated by a continuous boundary. This indicates that dmin is independent of N; however, the maximum spatial extent of the configuration is limited according to the condition in equation (2.14).

Heuristically determined minimum average spacing dmin, minimum area Smin, minimum extent Rmin and the minimum number of array elements Nmin required are shown in table 1 for different Np. For example, for a configuration with Np=6 traps, a minimum of Nmin=48 array elements are required to contain the area of Smin=8λ2 lying within a circular controllable area of radius Rmin=2λ, where Nmin is found by substituting Rmin into equation (2.14). Thus, any of the parameters listed in the table can be regarded as a function characteristic of the configuration complexity with a variable Np. Also, one can see that these functions grow as a function of Np. For estimation of the complexity, an approximate formula for the minimum extent Rmin can be used,

Table 1.

The minimum average spacing dmin, the minimum area Smin of the configuration, the minimum lateral extent of the configuration Rmin and the minimum number of array elements required for generation of the configuration Nmin as a function of trap number Np.

Np	dmin	Smin	Rmin	Nmin
2	λ		0.5λ	9
4	λ	λ2	0.5	12
6	2λ	8λ2	2	48
9	5λ	100λ2		120
12	12.2λ	900λ2	15	362
16	20λ	3600λ2	30	724

Finite-element model of the device

In §2, the underlying properties of a circular acoustic array were derived by assuming some pre-existing boundary conditions and by analysing their effect on the field in the volume enclosed. However, the question of how the required boundary conditions can be realized in practice was not addressed and is the subject of this section.

In a practical device such as the one shown in figure 1, the field in the fluid volume is produced by a piezo-ceramic transducer array. The field generated by each element in the array is evaluated using an FE model. This calculated field is then used to establish the excitation of the array elements required to produce a boundary condition such as determined in §2.

The FE model of the device was obtained using a commercial FE package (PZFlex; Weidlinger Associates Inc.). The device consists of a piezo-ceramic ring (PZ26, Ferroperm piezoceramics A/S) of inner radius 5.35 mm and wall thickness 1 mm, divided into 16 elements, separated by air-filled gaps and with a backing layer and a quarter wave (341 μm at 2 MHz) matching layer (15% by volume alumina-loaded epoxy).

A two-dimensional model, assuming plane strain in the perpendicular direction, was used and solved with a time-stepping algorithm. In order to reduce the size of the model, only half the system was modelled with a symmetrical boundary condition at the mid-plane x=0 and, rather than model a backing layer thick enough to absorb all incoming radiation, backing material (10% by volume tungsten-loaded epoxy) was modelled extending the modelled area out to a square of width 4 mm wider than the outer diameter of the piezo-electric array, and absorbing boundary conditions applied to remove reflections from the boundaries of the model. In a practical device (Courtney ), the backing layer of 9 mm thickness was used, absorbing 99 per cent of the incident energy.

Since all the array elements are assumed to be identical, it is sufficient to calculate the system response to a single excited element, e.g. the bottom element in figure 1. This element was excited with a single cycle at 2 MHz and the model run until pressure in the device had decayed to negligible levels. Steady-state responses (in terms of pressure complex amplitudes) were calculated by performing Fourier transforms of the response at each node of the model, extracting the complex amplitude at 2 MHz and normalizing to the applied voltage amplitude at the corresponding frequency.

Given the response of the system to excitation of one element (figure 6a), the response to excitation of multiple elements can be calculated as the linear superposition of the single-element response, rotated appropriately and scaled with a complex amplitude to allow different excitation amplitudes and phases to be applied to each transducer. The total field in the fluid chamber is thus given by where p(r) is the total pressure in the fluid chamber, p(r) is the pressure response of the system to the excitation of an element i and V is the excitation amplitude of the element.

Figure 6.

FE-based calculation of a force potential U in a 16 element device. (a) Single device excited at the bottom of the array. (b) All the elements excited to generate a trap at RT∼0.937λ distance from the centre.

According to §2, in order to generate the required pressure landscape p(r), it is sufficient and necessary to impose the corresponding boundary pressure p(r=R,θ) at some R. However, given only an N element device, it is only possible to satisfy this boundary condition approximately. To determine the optimal amplitudes (V) giving the best approximation, the quadratic average error function between the required boundary pressure p(R,θ) and the boundary pressure generated is defined as and the values of V are found by solving the set of 2×N linear equations derived from the condition The force potential distribution with a trap at RT∼0.937λ, corresponding to the maximum displacement distance according to equation (2.11), obtained using the outlined optimization procedure, is shown in figure 6b. Here, similar to the idealized case analysed in §2, the field inside the area defined by equation (2.11) is distortion free, and the field outside this area is dominated by the artefact field contribution.

This analysis demonstrates that it is possible to impose the boundary conditions necessary for creating pressure traps using a physically realizable ultrasonic array (figure 1) based on the fields predicted by an FE model. However, it is clear that the successful operation of a physical device using this method depends on how accurately the FE model represents the performance of each array element. For example, a real device is three dimensional, there will be inter-element variability and there may be other manufacturing imperfections leading to deviations between a numerical model and the true performance. The robustness of the proposed method depends on its insensitivity to such deviations, which are likely to impact on the physical design of the device. This is the subject of ongoing research.

Conclusions

A circular transducer array device (figure 1) for manipulation of particles was analysed. The purpose of the analysis was to examine the ability of the device to produce pressure landscapes required for particle manipulation. This proof of principle analysis was carried out in two steps. First, a general theoretical approach based on Green's function integral representation of volume pressure as a function of the boundary pressure was developed. This formalism has been used to establish the upper limit of controllable area of the operational volume as a function of the number of array elements and to evaluate the complexity price function as a function of the number of traps. Using the conclusions drawn from the formal approach, an FE model of a practically realizable device was used to obtain a preliminary assessment of its feasibility. Although a further in-depth modelling of the practical device is essential, the preliminary results presented indicate the feasibility of using such a device for particle manipulation.