Entropy production in an elementary, light driven micro-machine

We consider the basic, thermodynamic properties of an elementary micro-machine operating at colloidal length scales. In particular, we track and analyse the driven stochastic motion of a carefully designed micro-propeller rotating unevenly in an optical tweezers, in water. In this intermediate regime, the second law of macroscopic thermodynamics is satisfied only as an ensemble average, and individual trajectories can be temporarily associated with decreases in entropy. We show that our light driven micro-propeller satisfies an appropriate fluctuation theorem that constrains the probability with which these apparent violations of the second law occur. Implications for the development of more complex micro-machines are discussed.


INTRODUCTION
Advances in micro-fabrication techniques enable the manufacture of finely structured objects with ever greater precision [1]. As a consequence, increasingly sophisticated micro-machines are being developed [2][3][4] to perform a variety of previously unrealisable tasks in the colloidal regime. Examples include controlled transport, micro-fluidics, pumping and mixing [5][6][7][8][9], sensing [10], or even hydrodynamic manipulation [11]. More recently, self-organizing, dynamically reconfigurable and self-propelled micro-machines have been exhibited [12][13][14][15]. Quantitatively describing the behavior of machines at these length scales is not trivial. Classical thermodynamics was developed, in part, to help optimize the efficiency of conventional machines. This framework, however, is inadequate for the colloidal regime where the energy flows that drive the machines are comparable in size to stochastic thermal fluctuations. In this article, we analyze the influence of thermal fluctuations on an elementary micro-machine, in this case a light driven micro-propeller or light-mill, observing that they are constrained by an appropriate fluctuation theorem [16]. In doing so, we underscore the usefulness of applying the principles of stochastic thermodynamics [17] to the design of novel micro-machines.
As with biological analogues such as molecular motors [18], artificial micro-machines operate in an intermediate thermodynamic regime. Conventional machines, working at every-day length scales, conform to the laws of classical, macroscopic thermodynamics where physical processes are strictly irreversible and associated with increasing disorder, or entropy. In contrast, dynamical motion at microscopic length scales is governed by time symmetric equations of motion. The thermodynamics of micro-machines are neither purely reversible nor completely irreversible and the second law of thermodynamics, which prohibits entropy from decreasing with time, appears not as an absolute law, but as a limiting case, emerging only for sufficiently large systems, over sufficient time intervals [17,19]. For micro-machines, working in the colloidal regime, individual trajectories can be associated with decreasing entropy for finite periods of time. The relative probability of observing such a trajectory is constrained by the fluctuation theorem (FT), the term really applying to a family of theorems, and decreases exponentially with time, with the second law restored in the limit either of long times or large system sizes. For a general system the FT takes the following form: where Σ t is the entropy production of a trajectory integrated over time t (divided by Boltzmann's constant k B to make it dimensionless) and P(Σ t Σ) the probability density that Σ t takes the value Σ. Equation 1 therefore quantifies the relative probabilities of entropy producing and consuming trajectories. Since Σ t is extensive, it increases with the system size and averaging time. An immediate consequence of Eq. 1 is the second law inequality, 〈Σ t 〉 ≥ 0, where the average is taken over an ensemble of experiments with the same start time, averaged over a fixed time period t [19]. The second law inequality requires that the average entropy change is non-decreasing. Unlike the classical second law, it does not, however, prohibit decreases in entropy occurring in some of the realizations comprising the ensemble. Equation 1 is commonly referred to as the transient or detailed FT (DFT) to distinguish it from the integrated version, where P(Σ t < 0) 0 −∞ dΣ t P(Σ t ) etc. and 〈 . . . 〉 Σ t > 0 denotes an average over all trajectories for which Σ t > 0. Again, since Σ t is extensive, the right hand side of Eq. 2 approaches zero as the system size or time duration increase.
We wish to confirm appropriate versions of the fluctuation theorems, Eqs 1 and 2 for an elementary micro-machine. Previous tests of the FT involve dragging spherical beads through water with optical traps [20,21] under various conditions [22]. In these studies, a colloidal bead is dragged in a straight line, with a force that does not depend explicitly on position. By contrast, we use an autonomous micro-machine, with a more complex geometry and force profile. As will become clear, this increase in complexity gives rise to a number of experimental challenges, especially those relating to the design of our micro-machine and to the tracking of its motion.
Our chosen machine is a simple, light driven propeller or lightmill, which we fabricate using two-photon polymerization [23]. We track its motion as it rotates in a laser trap, showing that it conforms to an appropriate FT. The propeller has been carefully designed so that it rotates slowly about its axis, allowing its motion to be tracked with sufficient accuracy and resolution to observe trajectories with temporarily decreasing entropy. Small, natural asymmetries in the system cause the axial torque to vary strongly with orientation, so that the rotational motion is highly uneven. For the rotor used in this study, the ratio of the maximum to minimum torque is ∼ 2.6. Rather than being a defect of the study, these irregularities are an important part of it. The case of constant torque is comparatively trivial. In particular, diffusion of a particle exposed to a spatially invariant force or torque is very well understood [24]. In our case, the dependence of the torque on the orientation of the propeller provides a more demanding and interesting test of the FT.
In the following sections we review the necessary theoretical background, which closely follows the work of Speck, Seifert et al. [22], describe the design and fabrication of the rotor and the experimental procedure. We next present tests of the detailed and integrated forms of the rotational FT. We comment on the small discrepancies between theory and experiment and conclude with a discussion of the implications for the optimal design of future micro-machines.

THEORETICAL BACKGROUND
We consider the over-damped motion of a Brownian propeller or rotor, rotating about a fixed axis in water. This motion is described by the Langevin equation, where ϕ is the orientation of the rotor and ξ r is the rotational friction. The systematic, external torque, Γ(ϕ) varies with ϕ but is always of the same sign and Γ S (t) is the stochastic, Langevin torque with zero mean, uncorrelated and with variance determined by the fluctuation dissipation theorem, i.e., Thus, we focus on axial rotations, assuming that the motion of all other degrees of freedom are negligible or independent. Qualitative features of the motion follow from a consideration of the impulses, I, delivered to the rotor over a finite time interval, Δt, i.e., I Δt dtΓ(t). For short time intervals, the impulse from the systematic torque, I sys , is obviously proportional to Δt, i.e., I sys ≈ Γ(ϕ 0 )Δt, if ϕ 0 is the initial orientation and we assume the torque does not vary too much over Δt. By comparison, the impulse from the stochastic torque, I stoch , is a Gaussian random variable drawn from a distribution with a variance of 〈(I stoch ) 2 〉 2k B Tξ r Δt. The probability that the stochastic impulse exceeds the systematic impulse is therefore, In other words, for very small Δt → 0 we expect the total impulse to be dominated by the stochastic term, i.e., I ≡ I stoch + I sys ≈ I stoch . It may therefore act in the opposite sense to the systematic torque. As Δt increases, this becomes progressively less likely, until the total impulse is approximately equal to the contribution from the systematic torque. The angular displacements made by the rotor may be expected to behave similarly, i.e., for short times there is significant probability that the propeller rotates in the opposite sense to the applied torque while, for greater time intervals, the propeller is practically guaranteed to rotate with the systematic torque. Furthermore, if we consider the probability, p(ϕ, t)dϕ that, at time t, the propeller has an orientation in the interval [ϕ, ϕ + dϕ], it is clear that p(ϕ, t) must approach a steady state limit, p ss (ϕ), as the time interval over which data is collected increases, i.e., lim t → ∞ p(ϕ, t) p ss (ϕ). The time evolution of p(ϕ, t) is determined by the Fokker-Plank equation [25] with drift and diffusion coefficients derived from the Langevin equation, Eq. 3, where p(ϕ, t) is the probability distribution of the orientation and is the total probability current, the first and second terms of which correspond, respectively, to the drift and diffusion currents. As time progresses, p(ϕ, t) approaches a non-equilibrium steady state, p ss (ϕ) [26]. Since we are treating this as a one dimensional system, the associated steady state current, j ss and mean angular velocity, ω ss (ϕ) satisfy j ss p ss (ϕ)ω ss (ϕ), where j ss is independent of orientation in accordance with conservation of probability, i.e., As mentioned above, the torque, Γ(ϕ) does not change sign and the propeller rotates continuously. Obviously, this is a nonconservative system which is not at thermodynamic equilibrium. For such systems, we would usually expect thermodynamic properties (e.g., work done) to depend on the details of the stochastic trajectory of the propeller (ϕ(t) for t 0 ≤ t ≤ t 1 ), making accurate measurement challenging if not impossible. Fortunately, that is not the case here. Since the system is one dimensional we can define an effective potential from which the torque can be derived. If we write Γ(ϕ) T 0 + T 1 (ϕ), where T 0 is the average value of Γ(ϕ) in the range [0, 2π), then the corresponding potential is, Here, ϕ c is a continuously varying angular displacement that can become arbitrarily large and counts multiple complete rotations, i.e., it is the total angular displacement. The upper limit in the integral appearing in the final term of Eq. 9 can be changed from ϕ c , which measures the total angular displacement, to ϕ, the orientation of the propeller (i.e., ϕ ϕ c − 2πn for integer n ϕ c /2π ) since the integral of T 1 (ϕ ′ ) over a complete rotation is zero by construction. This step allows us to evaluate thermodynamic quantities in terms of the end points of trajectories, rather than having to measure the detailed motion. To illustrate, the dissipated heat is [22], where the integral in Eq. 10 is evaluated by writing the torque as the derivative of the potential in Eq. 9, and ΔV 0 V 0 (ϕ(t)) − V 0 (ϕ(0)) and ϕ again, is the orientation whereas ϕ c is the total angular displacement. Following Ref. 22, the total change in entropy along a trajectory is given by the sum of the change in the medium, and the change in the system entropy, where p ss (ϕ) is the steady state distribution of the propeller orientation. This latter quantity, the system entropy, is defined in analogy with the usual Gibbs entropy, with the integral over microstates replaced by an integral along the particle trajectory [27]. An additional term, related to the configurational entropy induced by the external field is significant for molecular systems [28], but negligible here. Since we have assumed that the sign of the torque is the same for all orientations (a condition satisfied for the propeller used in the experiments, see below), the total potential [Eq. 9] is monotonic. Trajectories associated with the consumption of medium entropy [Eq. 11] therefore relate to small rotations against the applied torque which, as discussed above, occur only over short time intervals. We note that Eq. 10 can be arrived at by various other routes. For instance, the torque can be separated into a part derived from a 2π periodic potential, V 0 say, and a constant, non-conservative term. By computing the work from an integral of the non-conservative torque, and the change in internal energy from V 0 , an application of the integrated first law reproduces Eq. 10 [22]. Of course, the force can be partitioned into conservative and nonconservative components in infinitely many ways, all leading to the same result. As stated above, the fundamental reason is that this is a one dimensional system so the force can be represented as the derivative of a scalar function, as in Eq. 9. The principles described above provide us with a simple way of testing the rotational fluctuation theorem for a driven micropropeller. First we measure the steady state orientation distribution, p ss (ϕ) and average rotational velocity, ω ss (ϕ). Next we evaluate the systematic torque from the steady state probability current, j ss , Eq. 8. This step requires the rotational drag, ξ r , which we estimate numerically [29]. We note that it is not possible to find an exact, independent value for ξ r since the precise dimensions of the micro-fabricated rotor are not available. The numerical estimate of ξ r is therefore treated as a guide and, since it is a scalar, it is straightforward to consider moderate variations in this parameter. Finally, Eqs 11 and 12 are used to evaluate entropy changes over time intervals of varying length and apply them in the integrated and detailed fluctuation theorems, Eqs 1 and 2. Achieving this experimentally is rather demanding, and the protocols and techniques we adopted are described in greater detail below.

PROPELLER DESIGN AND FABRICATION
Testing the rotational FTs described above places some demanding constraints on the design and behavior of our micro-propeller. It must exhibit three key characteristics. 1) The propeller should rotate about a single fixed axis, with minimal fluctuations in the orientation of the axis. 2) It should rotate slowly enough for the short term, entropy consuming trajectories to be accurately resolved by the available technology. Finally, 3) its shape should facilitate accurate tracking of its orientation. We discuss each of these issues in more detail below.
(1) Axis stability: Elongated, dielectric objects tend to align themselves with the axis of optical beams [30]. The stability of the rotation axis of a micro-propeller is therefore promoted by providing it with a long, sturdy central spindle.
(2) Rotation rate: Equation 6 allows us to find a crude estimate of the rotation rate required to observe entropy consuming trajectories, i.e., we require Δt where ω is the typical angular velocity at ϕ 0 . Taking Δt min to be the smallest time interval over which we can accurately measure gives ω 2 ≈ 4k B T/ξ r Δt min . Finally, if we approximate the rotational drag by that of a sphere, ξ r ≈ 8πμa 3 with an effective radius a 2 μm we get ω(0.1/ Δt min √ . In our case, Δt min 1 ms, suggesting that we require a rotation rate of <0.5 Hz. Most light-mills rotate at greater rates, typically in excess of several Hz [31][32][33], and more effort is devoted to increasing rotation frequency than decreasing it [34]. Designing a propeller that rotates this slowly without stalling completely is a greater experimental challenge than might have been expected. To meet it, we recall the basic symmetry properties of optical torque coupling [35]. An object with rotational symmetry and a mirror plane that includes the beam axis has a preferred orientation in a linearly polarized optical tweezer, while an object with chiral symmetry experiences a constant, orientation independent torque. We adopt a basic design for the propeller in which the torque it experiences should be controllable through a single parameter. This design consists of a propeller with two tiers of five cylindrical arms. When the tiers are aligned, the rotor has a mirror plane, and does not rotate. By gradually displacing the tiers with respect to one another, we can control the rotation speed. Natural imperfections in the beam and propeller result in a torque that, rather than being constant, is strongly dependent on orientation. This tendency becomes more conspicuous the closer the propeller is to having a mirror plane. (3) Tracking: Finally, it must be possible to track the orientation of the propeller with as much angular precision as possible.
To do this, we add microspheres to the ends of the rotor arms, and make use of the established methods for high precision sphere tracking [36]. The orientation of the propeller can then be accurately retrieved from measurements of the sphere positions.
With these principles in mind, we arrived at the following final design, shown in Figure 1.
We fabricate micro-propellers using two-photon polymerization. To implement the geometric design shown in Figure 1 (left), we use custom software (LabVIEW) to generate an instruction file to control the path of a laser beam in a two-photon polymerization machine (Photonic Professional, Nanoscribe). By specifying the beam path directly, rather than relying on slicing software to automatically generate a path from a CAD file, we retain fine control over the fabrication process. We deposit approximately 20 μL of photoresist material (IP-G 780, Nanoscribe) on a coverslip, and then heat it using a hotplate (100 C for 1 h), following the manufacturer's instructions for photoresist preparation. We then fabricate many (>50) copies of our micro-propeller design using the two-photon polymerization machine to scan a laser beam through the prepared photoresist according to our custom beam path file. The polymerized material is developed, and unpolymerized material washed away, by submersing the coverslip first in developing fluid (Microposit EC Solvent) for 20 min, followed by isopropanol for 5 min. We move the fabricated micro-propellers into suspension by placing a drop of water (Direct-Q 3 UV, Millipore) on top of them, and then apply mechanical agitation to detach them from the coverslip using a thin wire manipulated by a 3-axis translation stage (ULTRAlign Precision XYZ Linear Stage, Newport), similar to a previously employed method [37,38]. The suspension is transferred to a custom microscope sample chamber consisting of one microscope slide and three coverslips bonded together (UV Adhesive 81 or 68, Norland) using a pipette, and the remaining volume of the chamber is filled with additional water then sealed with adhesive. This process typically results in around 50-80% of the fabricated micro-propellers being transferred into the microscope sample chamber. The two-photon polymerization and transfer steps are all carried out in a clean room.

EXPERIMENT
We optically trap one of the micro-propellers inside the sample chamber using a custom-built holographic optical tweezers system, which is similar to a system described elsewhere [39]. The optical tweezers are based around a 1,070 nm laser (YLM-5-LP-SC, IPG Photonics) whose beam is expanded to fill a spatial light modulator, or SLM (XY Series HSP512-1064-DVI, Boulder Nonlinear Systems). The beam is then tightly focused using an oil immersion, 1.4 NA, 100× objective lens (Plan-Apochromat, Zeiss) to form optical trap(s). Relative translation of the sample and optics is achieved with a piezo-based translator (Mipos 140 PL, Piezosystem Jena) along the optical (z) axis, and with a servomotor-based stage (MS2000, ASI) in the orthogonal (xy) plane. Imaging of the sample is performed by illumination from a custom LED source (University of Glasgow). Light from the LED which is transmitted by the sample is then collected by the objective lens and directed to a camera (EoSens CL MC1362, Mikrotron) using a polarizing beam-splitter. Images from the camera are acquired on a computer through a PCI-e frame grabber (PCIe-1433, National Instruments). We use the "Red Tweezers" software [36] to generate and position optical traps through control of the SLM. A single micro-propeller is picked up with an optical trap, such that the long axis of the propeller's spindle aligns with the optical axis. The micropropeller is then translated to the center of the sample chamber in the xy plane, and to a distance of around 30 μm above the coverslip along the z axis using the microscope stages. Note that the latter distance is limited by the working distance of the objective lens. The micro-propeller is simultaneously confined in translational space and driven to rotate continuously about its spindle by the optical trap. We record a high speed (1 kHz) video of a trapped and rotating micro-propeller using a custom LabVIEW program to acquire images from the system's camera. The fast frame rate is achieved by i) reducing the active area of the camera to the minimum needed to view the propeller, ii) use of a signal generator (33220A, Agilent) to trigger the frame grabber, and iii) use of a circular buffer in software that losslessly transfers acquired images to an SSD drive (HyperX SH100S3/120G, Kingston). In our set up, the frame rate is limited by illumination intensity and control over the camera's exposure time.
Data analysis is performed offline, starting by extracting the position of the micro-propeller in each frame of the video using an additional custom LabVIEW program. We first find the Cartesian position of each of the five spheres on the micropropeller, using a symmetry transform tracking library, taken from the Red Tweezers software [36,40]. The propeller's orientation, ϕ, can then be defined by the Cartesian position of one of the spheres and the position of the propeller's center. The center is calculated in each frame as the average of the positions of all five spheres.
We record video of a rotating micro-propeller for a period of 15 min, equivalent to approximately 175 complete rotations, and extract the orientation in each frame as described above. The steady state probability distribution, p ss (ϕ) is constructed as a histogram by binning the available data. The bins have a width of one degree. We confirm that p ss (ϕ) is a true steady state by comparing with separate orientation distributions derived from the first and second halves of the total Brownian trail. These two distributions are the same, with the exception of a jagged, high frequency component. These sharp features are of very low amplitude in comparison to the overall shape of p ss (ϕ) and derive from multiple sources including the finiteness of the Brownian trail, tracking errors, center of mass motion and tilting. To prevent unrealistic contributions to the medium entropy introduced via the gradient of p ss (ϕ), we smooth the distribution function. To evaluate ω ss (ϕ) we form a histogram of the coarse-grained velocity, i.e., we take the numerical derivative of the orientation over a single time step and bin the result according to the average orientation over the time step. Figure 2 shows p ss (ϕ), ω ss (ϕ) and their product, j ss [see Eq . 7]. Smoothed data is shown in thin, colored lines with the raw data plotted on thicker gray lines. The small variation of j ss with orientation provides a measure of the errors in the experiment. Next, we calculate V(ϕ) using the integral in Eq. 9 and Γ ξ r _ ϕ (with ξ r 2.8 aN m s/rad). We then perform the arbitrary separation of V(ϕ) into a conservative potential, V 0 (ϕ), and our choice of T 0 (the average value of Γ(ϕ): T 0 -3.59 aN·m), see Figure 3.
Finally, we investigate total change in entropy over different time scales. We choose an interval, t, and split the time series of micro-propeller orientations from the experiment into T/t sections, each spanning time t. For each of these sections, the total entropy change is calculated from the sum of Eqs 11 and 12, using the end points of the trajectory and the previously calculated values for p ss (ϕ), T 0 , and V 0 (ϕ). We repeat this analysis for different intervals, t. Figure 4 shows histograms of entropy change, Σ t , over a sequence of time intervals, t. As suggested above, for the shorter time interval, t 1 ms, the distribution is strongly peaked close to zero, with a substantial fraction of trajectories FIGURE 2 | Steady state distributions of the orientation, p ss (ϕ), angular velocity, ω ss (ϕ) and probability current, j ss . The data are averaged over ≈ 175 revolutions, using bins of width 1°. In each case, the raw data (gray) is overlayed with the smoothed data. The scale used for j ss is numerically equivalent to that used for p ss (ϕ), with units Hz. showing negative changes in entropy. As the time interval lengthens, the distribution spreads out, and the peak shifts to higher positive values. At each stage the fraction of trajectories showing negative entropy change decreases. This process is quantified by the transient and integrated rotational FTs, which are tested in Figure 5. On the left hand side, Figure 5A we plot the logarithm of Eq. 1. In accordance with the FT, the graph is a linear one passing through the origin. The gradient of the best fit line is 1.14, compared with 1, as it would be for perfect agreement. The integrated form is tested by plotting the right and left hand sides of Eq. 2. The shaded region shows the effect of varying the value of ξ r by ± 15% when evaluating Γ(ϕ). Since the rotational friction scales with the cube of the length, this variation corresponds to a change in linear dimension of ≈ ± 5%. The strength of the agreement increases with the time interval, becoming very close for ta3 ms.

DISCUSSION
We have defined a version of the FT suitable for describing fluctuations of a rudimentary, light driven micro-machine. The scenario we have considered is somewhat contrived: we had to go to considerable lengths to ensure that we could reliably resolve the behavior of interest. In particular we had to design a propeller that would rotate slowly enough for the transient, entropy consuming trajectories to be accurately measured. The design also minimized motion of the center of mass, and angular fluctuations of the rotation axis, allowing us to eliminate five of the six physical degrees of freedom and focus on the remaining axial rotations. For the detailed FT, Eq. 1 and Figure 5A, the experimental data, which are plotted for the logarithm of Eq. 1 are represented by a straight line, indicating the exponential dependence expected from the FT, however the gradient is higher than expected. In the case of the integrated FT, Eq. 2 and Figure 5B, our measurements tend to under-estimate the entropy changes over very small time intervals, although the   Frontiers in Physics | www.frontiersin.org November 2020 | Volume 8 | Article 593122 6 agreement becomes more accurate for intervals of longer duration. There are a number of sources of error including image blur and fluctuations in the rotation axis of the micropropeller. Both of these factors are more significant for short time intervals. In the first case, each video image does not represent an instantaneous moment, but contains contributions from previous times. This has the effect of low pass filtering the measured rotations relative to actual rotations. Tilting of the rotation axis adds a small error to the measured displacement of each tracking sphere, when projected onto a fixed horizontal plane. In addition the effects of data smoothing (Figure 2) will be more significant for shorter intervals.
Nevertheless, the essential features are qualitatively reproduced. For short time intervals, entropy consuming trajectories occur with substantially higher probability, and this probability decreases exponentially as the duration of the time interval increases.
Although it was not straight forward to observe and quantify this behavior, we note that fluctuations influence the efficiency and precision of micro-machines whether they are measured or not. The concepts described in this article, and in the many articles devoted to both the fundamental theory, and to applications in biology, are directly relevant to the growing field of artificial micro-machines [2,3]. This will become increasingly more true as we try to design increasingly small machines. Indeed, the analysis of the micro-propeller is somewhat simplified by the fact that the underlying potential (Eq. 9; Figure 3) is many times greater than k B T. For this reason, the effect of diffusion in Eqs 8 and 11 is relatively weak in comparison to the drift terms produced by the systematic torques. Understanding the role of fluctuations in more complex micromachines, will present further challenges, especially for machines involving multiple degrees of freedom. Under these circumstances, the steady state distributions are themselves multi-dimensional, complicating the interpretation of force, locally. For instance, if we wish to design a micro-machine that applies a particular force in a particular configuration, we might also need to understand the morphology of p ss in multiple dimensions. Providing theoretical approximations for the effect of fluctuations on more complex micro-scale devices with many degrees of freedom, including synchronizing systems [41,42], presents an interesting challenge and one that could be technologically useful. One possibility would be to try to eliminate spurious degrees of freedom, by considering the cycle of the machine as a one dimensional curve embedded in a higher dimensional phase space. In fact, this is precisely what has been done for the rotor considered above. Ultimately, understanding the interplay between deterministic and fluctuating forces is essential for the design and optimization of complex, multi-dimensional, mesoscopic machines.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
SB designed and fabricated micro-propellers with support from DP, and performed all experiments, data analysis and figure preparation. SS supervised the project with theoretical support from MA. SS wrote the article with input from all authors. SB acknowledges funding from the EPSRC and DP thanks the Royal Academy of Engineering for financial support.