Oriented suspension mechanics with application to improving flow linear dichroism spectroscopy

Flow linear dichroism is a biophysical spectroscopic technique that exploits the shear-induced alignment of elongated particles in suspension. Motivated by the broad aim of optimizing the sensitivity of this technique, and more specifically by a hand-held synthetic biotechnology prototype for waterborne-pathogen detection, a model of steady and oscillating pressure-driven channel flow and orientation dynamics of a suspension of slender microscopic fibres is developed. The model couples the Fokker–Planck equation for Brownian suspensions with the narrow channel flow equations, the latter modified to incorporate mechanical anisotropy induced by the particles. The linear dichroism signal is estimated through integrating the perpendicular components of the distribution function via an appropriate formula which takes the biaxial nature of the orientation into account. For the specific application of pathogen detection via binding of M13 bacteriophage, it is found that increases in the channel depth are more significant in improving the linear dichroism signal than increases in the channel width. Increasing the channel depth to 2 mm and pressure gradient to 5 × 104 Pa m−1 essentially maximizes the alignment. Oscillating flow can produce nearly equal alignment to steady flow at appropriate frequencies, which has significant potential practical value in the analysis of small sample volumes.

Flow linear dichroism is a biophysical spectroscopic technique that exploits the shear-induced alignment of elongated particles in suspension. Motivated by the broad aim of optimizing the sensitivity of this technique, and more specifically by a hand-held synthetic biotechnology prototype for waterbornepathogen detection, a model of steady and oscillating pressure-driven channel flow and orientation dynamics of a suspension of slender microscopic fibres is developed. The model couples the Fokker-Planck equation for Brownian suspensions with the narrow channel flow equations, the latter modified to incorporate mechanical anisotropy induced by the particles. The linear dichroism signal is estimated through integrating the perpendicular components of the distribution function via an appropriate formula which takes the biaxial nature of the orientation into account. For the specific application of pathogen detection via binding of M13 bacteriophage, it is found that increases in the channel depth are more significant in improving the linear dichroism signal than increases in the channel width. Increasing the channel depth to 2 mm and pressure gradient to 5 × 10 4 Pa m −1 essentially maximizes the alignment. Oscillating flow can produce nearly equal alignment to steady flow at appropriate frequencies, which has significant potential practical value in the analysis of small sample volumes.
2019 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Introduction
Suspensions of particles in liquid or gas are found throughout the natural world and in many industrial processes; including blood, particulate-laden air, algae in open water or in bioconvection experiments, semen samples and industrial applications such as emulsions in food and cosmetics-these examples and more are elaborated in [1][2][3][4][5]. We will consider flowing suspensions of microscopic elongated rod-like particles, as a means to understand and improve flow linear dichroism spectroscopy, with particular focus on a specific technological application. The core idea of flow linear dichroism spectroscopy is that elongated particles undergo shear-induced rotation, which has the effect of concentrating their orientation in the direction of flow. The difference in absorbance of light polarized parallel, and perpendicular to, the flow direction can be used to reveal structural information; further information and wider context is given in [6][7][8][9]. The application of mathematical modelling to understand and quantify molecular orientation in flow linear dichroism was addressed by McLachlan et al. [10], building on classical oriented suspension mechanics [11,12]. This manuscript will generalize this work to the significantly more complex problem of a non-homogeneous shear environment of pressure-driven, and potentially time-varying, channel flow, and will moreover focus on a specific technological application. The methods and results will be adaptable to linear dichroism spectroscopy and beyond.
The specific application is a prototype hand-held device, developed by Linear Diagnostics Ltd (LD) designed to detect waterborne pathogens in fluids. The analyte is mixed with a reagent containing a synthetic biology micrometre-length fibre based on M13 bacteriophage, a filamentous virus known to infect Gram-negative bacteria (for example, Escherichia coli) [13]. M13 easily align in shear flow due to their relative rigidity and large aspect ratio, and are readily genetically engineered [14]. Microscopic elongated particles have been used to measure shear in microfluidic devices [15] and to determine wall shear stress in small biological structures [16]. In the absence of the pathogen, the fibres align in the flow field and therefore exhibit linear dichroism. The M13 are engineered with antibodies enabling binding to specific pathogens, therefore in the presence of the pathogen the alignment is disrupted, resulting in a reduced linear dichroism signal. The technique is described further in [14]. To improve the sensitivity of the device it is desirable to optimize the signal to noise ratio, which can be achieved by increasing the overall alignment in the system. To analyse the system theoretically involves solving a coupled pressure-driven steady or oscillating flow and alignment problem. While we will present results for a specific biotechnology application, the computational framework is generic to any channel flow of dilute suspensions.
The coupling between the flow and particle orientation is important to understand the dynamics of these suspensions. Jeffery [17] carried out a landmark study of the dynamics of suspended elongated particles in viscous flow, in particular, shear-induced alignment. Moreover, following seminal theoretical work [18][19][20][21], it has been established how the aspect ratio of particles affects emergent rheology. In the related context of suspensions of actively-swimming oriented particles, Pedley & Kessler [22] described a coupled model of orientation distribution, governed by a Fokker-Planck equation, and fluid flow, given by modifications to the Navier-Stokes equations proportional to the ensemble averages of orientation. These early studies form the theoretical framework for our model of homogeneous suspensions of semi-rigid M13 bacteriophage via the approximation of elongated, axisymmetric, rigid Brownian rods in the limit of infinite aspect ratio.
Other relevant theoretical analyses include work on homogeneous shear [12,23,24], timedependent shear [25], turbulent channel flow [26][27][28] and pressure-driven channel flow [29]. Including variations in particle concentration in channels has been investigated in a weak flow limit [30] and when the particle size is comparable to the channel width [31].
More recently, Taylor-Couette [32] and Rayleigh-Bénard [33] stability of perfectly aligned suspensions have been theoretically modelled via the transversely isotropic fluid of Ericksen [18,34]. Similar models have also been used to investigate fibrous biological systems such as growth in plant cell walls [35], propulsion in aligned cervical mucus [36,37] and the extracellular matrix [38]. This paper solves and compares two three-dimensional (3D) pressure-driven flow problems for an oriented suspension in a narrow rectangular channel geometry: steady flow, which models the continuously pumped loop used in existing spectroscopic devices, including the LD technology, and a potential novel oscillatory system, in which a smaller volume of analyte is oscillated back and forth. For practical particle concentrations, we have n * d a * 3 < 1, where n * d is particle number density and a * is the particle length, and hence the suspension is at the upper end of the dilute range [39]; this theory has previously been shown to hold for n * d a * 3 ≈ 1 [30]. The channel dimensions, pressure gradient, particle number density and frequency of oscillations are investigated as factors to improve alignment, and thus signal, in both flow types and to determine the viability of an oscillatory system for aligning particles.
The coupled orientation and flow model will be solved numerically by iterative coupling. Mathematical modelling of these suspensions is computationally challenging due to the additional independent variables associated with the particle degrees of freedom and the coupling between particle dynamics, velocity gradients and rheology. Rational simplification of the flow problem via lubrication theory, the application of a spectral method [12] and multicore parallelization of the array of spatially local orientation problems, will be shown to enable practical solution with workstation hardware.
The manuscript is organized as follows: the governing equations for the system, including the Navier-Stokes and Fokker-Planck equations, are summarized in §2. The steady flow model is presented in §3 and the oscillatory problem in §4. Results for both flow problems are presented in §5 and discussed and compared in §6.

Summary of equations governing dilute suspensions of elongated particles
Consider a 3D rectangular channel of width 2W * , depth 2h * and length-scale L * , where h * , W * L * . The axis origin is at the centre of the channel and it is assumed that the flow direction induced by pressure gradient G * , or the molecular orientation axis, is the x * -direction (figure 1a). The behaviour of dilute suspensions is approximated by altering the constitutive equation for the stress tensor σ * [40,41]; the motivating problem is at the upper end of the dilute regime (see §5). This constitutive equation, along with the coupled Fokker-Planck equation governing fibre behaviour [12], are briefly detailed below.
The dimensional Navier-Stokes and incompressibility equations are and where ρ * is density, u * (x * , t * ) is the velocity vector for spatial coordinate x * and time t * and ∇ * is differentiation in positional coordinates. The asterisk notation represents dimensional variables. The boundary conditions are given by no-slip and no-penetration at the walls, The boundary conditions for x * have not been defined here as they are specified for each flow system.
To determine the bulk stresses in the fluid, consider a homogeneous suspension of particles, dilute enough to ensure that particle-particle interactions are negligible [10,12]. Nitsche & Hinch [42] found that anisotropic translational diffusivity can produce migration and hence non-uniform concentration. In the present motivating problem, we estimate the timescale for translational diffusion as h 2 /D T ≈ (3 × 10 −4 ) 2 /10 −13 10 5 s and hence this effect will not be significant here. Furthermore, assuming the cell density is spatially homogeneous, the suspension   Figure 1. Schematic of the flow of a fibre-laden fluid, forced by pressure gradient G * with a distribution of orientation, in a thin rectangular channel. (a) Rectangular channel of width 2W * and depth 2h * and the x * -axis is the molecular orientation axis. The orientation parameter is calculated as the proportional difference in orientation parallel and perpendicular to the molecular orientation axis, highlighted by the arrows to the left of the channel. Note that for oscillating flow, the pressure gradient is replaced with G * e iω * t * for frequency of oscillations ω * . (b) The orientation angles θ ∈ [0, π ] and φ ∈ [0, 2π ) away from the coordinate axis describe the orientation of the fibres. configuration can be described by the orientation distribution function ψ(x * , p, t * ), satisfying the Fokker-Planck equation, where p is the unit orientation vector, p = (sin θ cos φ, sin θ sin φ, cos θ), (2.5) for θ ∈ [0, π ] and φ ∈ [0, 2π ) shown in figure 1b, ∇ p represents differentiation in orientation space and D * r is the rotational diffusion constant. The particle rotation due to shear is given by (2.6) where e * = (∇ * u * + ∇ * u * T )/2 is the rate of strain tensor, ω * = (∇ * u * − ∇ * u * T )/2 is the vorticity tensor ((∇ * u * ) ij = ∂u * i /∂x * j ) and α 0 = (r 2 − 1)/(r 2 + 1), for particle aspect ratio r, is a measure of slenderness [17]. The orientation distribution function must satisfy periodicity and the normalization condition, and s ψ dp = 1, where s is the surface of the unit sphere. Following Pedley & Kessler [22] (and the prior work of Batchelor [20] and Hinch & Leal [41]), the enhanced stress tensor for a suspension of inactive particles is where the Newtonian contribution is for pressure P * , identity matrix I and viscosity μ * . Here, σ * D is the extra stress due to the rotary diffusion of the particles and σ * P represents the interaction of the particles with the fluid, given by ψ dp (2.11) and σ * P = 4 μ * Φ α 2 e * : s p p p p ψ dp + α 3 e * · s p p ψ dp + s p p ψ dp · e * + α 4 e * s ψ dp + α 5 I e * : s p p ψ dp , (2.12) and where Φ is the volume fraction of particles. The constants α i (i = 2 . . . 5, r) relate to the aspect ratio of the particle and can be represented in terms of ellipsoidal integrals [17,20], which are detailed in equations (S1)-(S5) in electronic supplementary material, S1. For reference, Holloway et al. [43] have summarized this system of equations for particle suspensions in concise notation. All results in this paper include the full particle stresses σ D and σ P . Because the width and depth of the channel are comparable and much smaller than the channel length, we use a lubrication approximation, introducing the dimensionless parameter δ = h * /L * 1 to allow detailed analysis of the significant hydrodynamic interactions. The variables are scaled via where G * is a reference pressure gradient. For steady flow, G * is the pressure gradient in the channel and for oscillating flow it is equal to the amplitude of the oscillating pressure gradient. The full dimensionless models for steady and oscillatory flow are described in the relevant sections. Finally, the degree of orientation in a system can be described by an orientation parameter, S, calculated as = sin 2 θ cos 2 φ − sin 2 θ sin 2 φ , (2.14) for a 3D suspension where Note that this is different from the definition used by McLachlan et al. [10], where θ S = arccos(| sin θ cos φ|).

Steady flow
The first problem we consider is flow forced by a constant pressure gradient G * . §3a contains the Fokker-Planck equation and §3b contains the model for the flow. The full coupled problem is solved by iterating between the numerical methods described in §3a,b.

(a) Model for the orientation distribution function
The particles in the flow are small enough that each particle at a point in space is subject to an approximately linear velocity field in the local frame, and hence a constant shear rateγ * = (∇ * u * ) 2 + (∇ * u * T ) 2 . Therefore, the orientation distribution of the bacteriophage is formulated at each point in space independently; fluid flow influences orientation distribution only via the local Péclet number. We describe the model and its numerical discretization below. Applying lubrication theory, assuming α 0 = 1 and investigating steady-state behaviour, the Fokker-Planck equation reduces to where P e =γ * /6D * r is the local Péclet number and and where θ and φ describe the local coordinates at each point in space, which are calculated by rotating φ and θ away from the laboratory frame by an angle β = arccos(|∂u * /∂y * |/γ * ). The above equation is similar to the spatially uniform case presented by Strand et al. [12] (see also Bird et al. [45]), with the complication of requiring a local coordinate system. The Fokker-Planck equation (3.1) is spatially discretized via spherical harmonics [10,12,45]. Assume the solution takes the form where P m n are the associated Legendre polynomials. By substituting (3.4) into equation (3.1), Λ and Υ are given by Λ(P m n (cos θ ) cos(mφ )) = −n(n + 1)P m n (cos θ ) cos(mφ ), and Υ (P m n (cos θ ) sin(mφ )) = Strand et al. [12], are given in electronic supplementary material S2. The resulting system of equations for A m 0 n and A m 1 n , found by equating coefficients of sine and cosine, is

(b) Model for the velocity
Under the lubrication theory approximation, the flow problem can be expressed entirely as an elliptic partial differential equation for u(y, z) where the coefficients depend on moments of the orientation distribution ψ(φ, θ; y, z).
Scaling the system, via (2.13), neglecting terms O(δ) and smaller, the continuity equation (2.2) is unchanged and the Navier-Stokes equations (2.1) reduce to for global Péclet number P G = G * h * /D * r μ * . The y-and z-components of the dimensionless Navier-Stokes equations enforce P = P(x) and so the pressure gradient reduces to −1. The boundary where W = W * /h * . Equation (3.11), along with boundary conditions (3.12), are solved numerically using a finite difference scheme and constructing a matrix of equations from the resulting problem. The numerical solution to the Newtonian and suspension models are described in electronic supplementary material, S3, and the iterative process is described in electronic supplementary material, S4.
To optimize the efficiency of the above processes, the convergence of the spherical harmonic solution and the finite difference approximation are investigated for a range of spherical harmonic modes N and finite-difference grid spacings Y and Z. The full details and related figures are presented in electronic supplementary material, S3; we find that N = 10, Y = 150 and Z = 120 are sufficient for convergence. We will discuss the results in §5.

Oscillatory flow
The next flow system of interest consists of an oscillating pressure gradient where the length is still much larger than its width and depth but a smaller volume of fluid oscillates back and forth. The pressure gradient is now defined as for frequency of oscillations ω * . The dimensional scalings are given by equation (2.13), and the time is scaled against the frequency of oscillations as t * = t/ω * . The model for the time-dependent Fokker-Planck equation is given in §4a and the model for the velocity profile is given in §4b. The full coupled problem comes from iterating numerically between 4a,b over each time step.

(a) Model for the orientation distribution
Assuming the particles are small enough that at each point in the channel they are subjected to a spatially linear velocity field and thus a shear rate dependent only on timeγ * =γ * (t * ); the orientation distribution is again formulated at each point in space independently. Following §3a, the Fokker-Planck equation is where P e (t * ) =γ * /6D * r is the local Péclet number, t * = τ/6D * r and Λ and Υ are given by (3.2) and (3.3), respectively; note the introduction of a local time-scale τ .
The orientation distribution is discretized spatially in the basis of spherical harmonics, This is then substituted into (4.2), where Λ and Υ are given by (3.5)-(3.8). Equating coefficients of sine and cosine, the resulting system of ordinary differential equations for A m 0 n and A m 1 n is for integer N, where q = 0, 2, . . . , N and p = 0, 2, . . . , q. Here, A 0 0 0 (τ ) = 1/4π to satisfy the normalization condition, A 0 1 n = 0 for all n and A m 0 n = A m 1 n = 0 for odd n, m. The ordinary differential equations (4.4) and (4.5) are solved numerically using a second-order explicit Improved Euler method by first writing the equations in matrix form as described in §3a.

(b) Model for the velocity
Again using a lubrication theory approximation, the velocity problem can be expressed as a timedependent partial differential equation where coefficients of u(y, z, t) depend upon moments of the orientation distribution function ψ(φ, θ ; y, z, t). Scaling the Navier-Stokes and continuity equations, (2.1) and (2.2) respectively, retaining terms larger than O(δ), the continuity equation is unchanged and the dimensionless Navier-Stokes equations become for global Péclet number P G and Womersley number α 2 = ω * h * 2 ρ * /μ * . It is assumed the system is initially at rest, accounted for by altering the dimensionless pressure gradient to ∂P/∂x = exp(i(t + π/2)). Equation (4.6) is solved via an alternating direction implicit (ADI) method in which the time step is split in half and one spatial variable is discretized in each time step [46]; mixed derivatives are treated explicitly throughout this analysis. As in §3, a matrix equation is constructed for each time step and solved numerically; the details of this method can be found in electronic supplementary material, S5, and a description of the iterative procedure is described in electronic supplementary material, S6.
A numerical convergence study was undertaken, detailed in electronic supplementary material, S5. As in §3, Y = 150 and Z = 120 are satisfactory in the finite difference scheme and T = 1200 time steps ensure convergence.

Results
The velocity, orientation parameter and linear dichroism signal are compared for a range of particle number density, channel dimensions and pressure gradient. We will refer to the standard parameter set for these quantities as: channel half-depth h * = 3 × 10 −4 m, half-width W * = 1 × 10 −3 m, fixed pressure gradient G * ≈ 4.4 × 10 3 Pa m −1 , particle number density n * d = 7.33 × 10 17 phage m −3 and frequency of oscillations ω * = 10 rad s −1 (table 1(a)). Using a particle length of a * = 800 nm, the value n * d a * 3 = 0.3753 < 1 and so we argue that the approximation of a dilute regime is reasonable. We will also examine the effect of higher particle concentrations, which is beyond the strict regime of validity of dilute theory but with the aim of providing insight into this tractable regime. The dimensionless groups are detailed in table 1(c) and are subject to change depending on their composition of dimensional values; we do not directly vary the dimensionless groups. All other parameters are unchanged (table 1(b)).
The notation • and • y * represent the spatial average and the width average for steady flow, respectively, The linear dichroism signal is assumed proportional to the particle number density and is calculated as the root mean square (denoted as S rms ) of the orientation parameter across the   depth of the channel. The signal is given by where N * A is the Avogadro constant (≈ 6 × 10 23 phage mol −1 ), m * w is the molecular weight of M13 bacteriophage (m * w ≈ 2 × 10 7 g mol −1 from Linear Diagnostics Ltd) and * is the extinction coefficient (m 2 g −1 ), which is a measure of the absorption of light through the sample [7]. The extinction coefficient is taken as 0.38 m 2 g −1 at light wavelength 269 nm [47][48][49]. A number of the parameters in table 1 and their calculation are discussed in electronic supplementary material, S7. As a direct comparison with steady flow, the maximum value of the orientation parameter in the oscillating channel is presented.
The computational walltime for this problem was large, especially for oscillating flow. The numerical solution was parallelized using parfor in MATLAB  (a) Steady pressure gradient Figure 2 depicts the velocity profile, orientation parameter and linear dichroism signal for steady flow through the cross section of the channel, with the standard parameter set. As may be expected from the Newtonian solution [50], the velocity u * increases from zero at the channel walls to a maximum at the centre point ( figure 2a), corresponding to a minimum in the orientation parameter S (figure 2b). Associated with the large shear rate, the linear dichroism signal is largest near the channel walls and decreases to a local minimum at the centre (figure 2c); there is an increase at the middle width of the channel (y * = 0) where larger shear near the walls at z * = ±h * increases the orientation parameter. Local minima also exist at the channel corners. The orientation distribution function ψ is displayed for varying points in the channel cross section (figure 3). At the centre of the channel where the shear rate is zero (figure 3b) there are no biasing effects on the particles and so ψ is almost uniform in orientation space. Close to the walls at y * = ±W * the behaviour of ψ is similar ( figure 3c,d); the particles preferentially align in the flow direction, for orientation angles φ = 0, π and 2π and θ = π/2. Close to the corner of the channel this behaviour is more prominent. Near z * = h * (figure 3e), the M13 bacteriophage are more likely to be aligned at φ = 0 and π (note that the particles have no polarity so the two angles are equivalent) and θ = π/2. The flow profile accounting for a particle suspension at standard concentration is compared with the Newtonian profile u * N for the same flow set-up in figure 4. The actual difference shows that the profile for both flow systems is qualitatively the same; the suspension model produces a lower velocity however the change is small (figure 4a). To confirm that this change is small, the reduction in velocity relative to the magnitude of the Newtonian velocity is also calculated in figure 4b and the change is on the order of 10 −3 . This effect has been investigated for higher particle concentrations (electronic supplementary material, figure S4); above n * d ≈ 10 18 phage m −3 there is a greater reduction in the flow.
We compare the width-averaged orientation parameter calculated in this work (solid line) with equation (   In figure 6, the impact of the pressure gradient G * and particle number density n * d on the spatially averaged orientation parameter S and linear dichroism signal LD is examined. Increasing the number density reduces the orientation parameter; the particle effects become significant above n * d ≈ 1 × 10 18 phage m −3 . This effect is shown in more detail in electronic supplementary material, figure S4, in which the spatially averaged velocity profile and orientation parameter are calculated for varying number density. The velocity and orientation parameter are reduced approximately 30 and 19%, respectively, over three orders of magnitude. Increasing the number density significantly increases the linear dichroism signal due to their linear relationship (figure 6b).
The spatially averaged orientation parameter is calculated for a range of channel widths and depths in figure 7; the other parameters, including pressure gradient, remain fixed. The channel depth has a large impact on S; the orientation parameter increases as the depth is increased. The channel width has limited effect, unless the two dimensions are comparable (figure 7b). For channel widths below W * = 1 × 10 −3 m increasing the channel width increases the orientation parameter; this effect reverses as W * grows. Finally, the upper asymptote of the spatially averaged orientation parameter is investigated for large pressure gradient and channel depth ( figure 8). In both cases, S quickly tends towards a constant value, S ≈ 0.74.

(b) Oscillating pressure gradient
The velocity, orientation parameter and linear dichroism signal for oscillatory flow with the standard parameter set are shown in figures 9 and 10. Figure 9a,b shows how the velocity changes in time across the width and depth of the channel, respectively (away from edge effects in the other coordinate); the velocity oscillates and is zero at the walls as expected. A time-averaged velocity profile is plotted in figure 9c, this is qualitatively similar to the steady velocity profile in figure 2a but the magnitude of the velocity is reduced.
The orientation parameter is shown over time and for multiple values of z * in figures 10a for y * = 1.9 × 10 −3 m, and 10b for y * = 3.5 × 10 −4 m. Towards the middle width (figure 10a), the orientation parameter increases with the channel depth and in the centre of the channel (figure 10b) the opposite occurs. Figure 10c   dichroism signal in figure 10d. Note here that the peak values of the orientation parameter are similar for steady and oscillatory flow. Next, the spatially averaged orientation parameter S is compared for various pressure gradient and frequencies of oscillation ( figure 11). As the pressure gradient G * increases, the orientation parameter increases; there is an indication that the impact of G * on S is reducing as G * becomes large. The impact of increasing the frequency of oscillation ω * is to decrease S, a result seen throughout all figures. The relationship between the pressure gradient and frequencies of oscillation, for a choice of channel depth, is summarized in table 2.
The orientation parameter is calculated for increasing channel depth h * and for a number of frequencies of oscillation and channel width ( figure 12). Increasing W * , shown by line markers in figure 12a, decreases the orientation parameter in general. The orientation parameter has a nonmonotonic response to h * as the frequency of oscillations is increased; for small ω * , increasing h * increases the orientation parameter whereas for large ω * this relationship begins to reverse. This trend is highlighted in figure 12b, in which W * = 2 × 10 −3 m. The relationship between h * , W * , ω * and the orientation parameter S has been summarized in table 3; when ω * is large, a larger reduction in S is seen with increasing channel width.
To further investigate these properties, the orientation parameter is plotted for increasing channel width for a range of channel depth and frequencies of oscillations ( figure 13). As seen previously, the orientation parameter S increases with channel depth and decreases with increasing channel width. Increasing the frequency of oscillations reduces the orientation parameter in general; for larger channel depths, the decrease in the orientation parameter with ω * is greater.

Discussion
The coupled pressure-driven flow and orientation dynamics of a dilute suspension of elongated particles has been modelled, to provide a mechanistic underpinning for flow linear dichroism spectroscopy. The model was applied specifically to a prototype hand-held device for waterborne   pathogen detection, to explore the effect of changing the channel dimensions and pressure gradient, and to assess the feasibility of oscillatory flow. Mathematically, the problem involved coupling the Navier-Stokes equations, simplified via lubrication theory, to the Fokker-Planck equation for the shear-induced orientation dynamics. Whereas the flow problem involved two spatial independent variables and time, the orientation problem involved two angular independent variables and time, applied at each spatial location in the domain; the latter follows from the assumption that the particles are sufficiently small that they encounter a locally uniform shear rate. The flow and orientation problems were coupled via corrections to the fluid constitutive law involving moments of the orientation distribution, as described by established orientation Brownian suspension mechanics theory.
The system was solved numerically using spherical harmonics for the orientation distribution and a finite difference method for the fluid dynamics, with iterative coupling; the oscillatory problem was temporally discretized with a second-order explicit method for orientation dynamics and an alternating direction implicit method for the fluid dynamics. The oscillatory flow is rather more computationally expensive due to temporal dependence; a parameter set of 175 tuples was explored for the oscillatory flow, with each tuple requiring 14-20 h of walltime on five workstation cores. The dependence of particle orientation on channel width, depth, pressure gradient, particle number density and frequency of oscillation were reported and compared for both the steady and oscillatory systems.
To model linear dichroism, the analysis focused on a single orientation parameter, calculated as the average of the difference between the components of the particle alignment parallel to, and perpendicular to, the flow direction. It was noted that the earlier (homogeneous shear) analysis of McLachlan et al. [10] applied an alignment formula based on uniaxial orientation for elongated particles embedded in a membrane. However, because flow linear dichroism produces a biaxial orientation distribution-in other words, there is no rotational symmetry about the flow direction in general. The correct biaxial formula was compared with the uniaxial formula; it was found that the qualitative behaviour was similar, however, the uniaxial definition overestimates the degree of orientation.
In both steady and oscillatory flows, increasing the pressure gradient increases the shear rate and hence the degree of alignment in the system. For steady flow, a plateau was observed close to a -perhaps maximum achievable -value of S = 0.75, above G * = 5 × 10 4 Pa m −1 and h * = 2 × 10 −3 m, whereas for oscillating flow the mean orientation did not exceed S = 0.43 for the range of parameters analysed. This result does however suggest that oscillatory flow may be a viable method for producing practically useful levels of alignment, which in turn may be valuable in the analysis of small sample volumes. Alignment decays monotonically with driving frequency, however an increase from ω * = 10 rad s −1 to 40 rad s −1 results in less than a 25% reduction in alignment. Increasing channel height has a significantly greater impact than increasing channel width, although for oscillating flow there is evidence of a local maximum in alignment as channel height is increased. For a given limited sample volume, there is a clear trade-off between channel volume and oscillatory frequency; as channel volume is increased, a higher oscillatory frequency is needed to keep the sample underneath the detection window. These data reported should enable the optimal configuration for a given volume to be calculated.
Signal production was also modelled by multiplying the alignment parameter by the number density of particles; at lower volume fractions the effect of number density on flow and alignment is weak, so the relationship between number density and signal is predicted to be approximately linear. At higher volume fractions (n * d 10 18 ) the relationship is slightly weaker than linear due to flow and alignment attenuation. While increasing the number density may therefore appear to be a simple way to optimize signal, in practice there are trade-offs with scattering-which was not modelled in the present case-and cost of the reagent. Optimizing the remaining design parameters is therefore important.
Avenues for further work include experimental testing of the predictions of how signal should vary with channel dimensions and pressure gradient, and direct testing of oscillating flow linear dichroism spectroscopy. It will also be important to take account of scattering and hence produce a more accurate model of how number density affects signal production. Other avenues include incorporating flexibility of particles-since the persistence length of the M13 bacteriophage considered here is approximately 1.3 µm [51] and the fibres are of length 0.8 µm a rigid fibre approximation is reasonable, however, future work could incorporate the flexibility of the fibres. Depending on particle sizes and time scales, migration induced by anisotropic translational diffusivity can be important [42]; including spatial variations in the concentration would, therefore, be an interesting extension to the current problem.
While the present manuscript focused on flows of passive particle suspensions in thin rectangular channels, future work may include investigating active suspensions, which can display collective behaviour [52][53][54][55] and superfluidity [56]. The system considered here can be augmented by an extra stress term to account for self-motility [22], ψ dp, (6.1) which could model swimming elongated cells such as spermatozoa in pressure-driven flow or complex behaviours such as trapping in shear flow [57]. The dynamics of suspensions of elongated particles provides an enduring subject in fluid mechanics, with novel application areas continuing to be discovered. The challenging task of solving the complex models associated with non-homogeneous shear and time-dependent forcing is now becoming within reach of multicore computing hardware. We hope that the present manuscript provides both useful information in the design and optimization of the technological system considered, and an effective framework which can be applied and extended to other systems in biophysical spectroscopy and cell analysis.
Data accessibility. Two CSV files, containing the numerical results for each figure for steady and oscillatory flow, are given in the electronic supplementary material. Both CSV files contain dimensionless results and the dimensional results included in this manuscript.