The aerodynamics and power requirements of forward flapping flight in the mango stem borer beetle ( Batocera rufomaculata )

The need for long dispersal flights can drive selection for behavioral, physiological and biomechanical mechanisms to reduce the energy spent flying. However, some energy loss during the transfer of momentum from the wing to the fluid is inevitable, and inherent to the fluid-wing interaction. Here, we analyzed these losses during the forward flight of the mango stem borer ( Batocera rufomaculata ). This relatively large beetle can disperse substantial distances in search of new host trees, and laboratory experiments have demonstrated continuous tethered flights that can last for up to an hour. We flew the beetles tethered in a wind tunnel and used high-speed videography to estimate the aerodynamic power from their flapping kinematics and particle image velocimetry (PIV) to evaluate drag and kinetic energy from their unsteady wakes. To account for tethering effects, we measured the forces applied by the beetles on the tether arm holding them in place. The drag of the flying beetle over the flapping cycle, estimated from the flow fields in the unsteady wake, showed good agreement with direct measurement of mean horizontal force. Both measurements showed that total drag during flight is ~5-fold higher than the parasite drag on the body. The aerodynamic power estimated from both the motion of the wings, using a quasi-steady blade-element model, and the kinetic energy in the wake, gave mean values of flight-muscle mass-specific power of 87 and 65 W kg muscle -1 , respectively. A comparison of the two values suggests that ~25% of the energy is lost within the fluid due to turbulence and heat. The muscle mass-specific power found here is low relative to the maximal power output reported for insect flight muscles. This can be attributed to reduced weight support during tethered flight or to operation at submaximal output that may ensure a supply of metabolic substrates to the flight muscles, thus delaying their fatigue during long-distance flights. ~30 mW, respectively, in two hawkmoths flying at 2.5 ms -1 . With a mean body mass of 2 gr and flight varying between 21% and 34% of the hawkmoths’ body mass, the muscle mass-specific power of the moths was between 29.4 and 47.6 W kg -1 muscle for the PIV measurement, and 44.1 -71.4 W kg -1 muscle for the blade-element, respectively. These values are somewhat lower than our 71 and 84 W kg -1 derived from the PIV and blade-element approaches. The difference may correspond to the lower wingbeat frequency and wing loading of the moths as well as their more streamlined body. Using a blade-element approach, Dudley and Ellington (1990a) estimated body mass-specific power of forward flying bumblebees. Their results indicate 20-35 W kg -1 for flight at 2.5 ms -1 . With a flight-muscle ratio of 0.37 out of the total body mass in Bombus terrestris (Polidori and Nieves-Aldrey 2015), this gives a muscle mass-specific power of 54-95 W kg -1 muscle. These values are somewhat higher than ours but they include the inertial power of the wing mass under the assumption of perfect elastic energy storage; whereas our estimates include the inertial power of the added mass but not the inertial power of the mass of the wing itself. Ellington (1985) reported aerodynamic power values between 167-186 W kg -1 muscle for hovering honey bees, bumblebees and drone flies. Hence the muscle mass-specific power found here is within the range expected for flying insects. Considering the fact that weight support of our beetle was below the level required for free-flight, the submaximal power output of the flight muscle is to be expected. Submaximal muscle power output of the flight muscles might also delay fatigue during prolonged flight.


Introduction
Flight can greatly improve an animal's dispersal capacity (Johnson, 1969). It allows terrestrial organisms to quickly move substantial distances across direct routes that are not limited by obstacles, surface texture or topography. However, flapping flight is an energetically demanding form of locomotion, requiring high levels of power output that typically exceed the levels required for running and swimming in similarly sized animals (Butler, 2016). Consequently, animals engaged in long flights (e.g. migrations) often tend to resort to power and energy saving strategies that conduce to improved flight range and time in air, for a given amount of fuel (Hedenstrom and Alerstam, 1998;Pennycuick, 1969;Rayner, 1999;Tucker, 1973).
Being the smallest active flyers, insect flight speed is typically lower compared to that of birds and bats and their flapping frequency is higher (a consequence of having smaller wings). As a result, most of the wing speed required for aerodynamic force production is provided by the work of flight muscles, rather than by the forward speed of the body. The unique flapping kinematics and the low advance-ratio (forward flight speed/wingtip flapping speed) of insects compared to larger flying vertebrates is bound to affect energy expenditure during long distance flights. However, our understanding of the aerodynamics, power requirements and energetics of insects during prolonged forward flight is rather poor. To some extent this is because insect flight utilizes unique unsteady aerodynamics (reviewed by Sane, 2003) and most of our understanding of these unsteady mechanisms, and our ability to model them, pertain to hovering and low-speed flight. For instance, theoretical models for bird flight predict a "U-shaped" speed-power curve (Pennycuick, 1969;Rayner, 1999): i.e. the power requirements are high in both slow and fast flight and have a minimum at intermediate flight speeds. In contrast, empirical measurements on bumblebees found no significant change in oxygen consumption when flight speed increased from 0 -4 ms -1 (Ellington et al. 1990). A blade-element analysis of the wing and body kinematics of forward-flying bumblebees (Dudley and Ellington, 1990a) and computational fluid dynamics simulations (Wu and Sun, 2005) confirmed a low minimum in the power-speed curve that does not greatly differ from the power required for hovering flight. A similar result was found for fruit flies, using numerical simulations (Sun and Wu, 2003). In contrast, the power requirement of hawkmoths reveals a U-shaped power-speed relationship (Warfvinge et al., 2017; Downloaded from https://academic.oup.com/iob/advance-article/doi/10.1093/iob/obaa026/5902830 by guest on 18 September 2020 Willmott and Ellington, 1997). Thus, it seems that the power requirement for insect forward flight is hard to predict based on simple models, necessitating direct measurements of its value at the specific flight speed.
To study the aerodynamics of insect forward flight researchers have used computational fluid dynamics (e.g., Harbig et al., 2014;Lee et al., 2018;Sun and Wu, 2003;Wu and Sun, 2005), quasi-steady modelling based on flapping kinematics (Dudley and Ellington, 1990b;Dudley and Ellington, 1990a;Jensen, 1956;Weis-fogh, 1964), qualitative flow visualization using smoke released near free-flying insects (Srygley and Thomas, 2002;Thomas et al., 2004) and Particle Image Velocimetry (PIV) studies of free-flying or tethered insects in a wind tunnel (Bomphrey et al., 2005;Henningsson and Bomphrey, 2011;Johansson et al., 2012;Johansson et al., 2013;Warfvinge et al., 2017;Young et al., 2009). The PIV technique, in particular, enables quantitative measurement of the unsteady flow-field in the wake of the flying insect. Since the flying insect transfers momentum from the flapping wings to the surrounding air, its unsteady wake can be used to link between the flapping kinematic of the wings and the aerodynamic forces and flows generated during flight.
The mango stem borer (Batocera rufomaculata De Geer) is a large beetle (body mass 2-7 gr) from Southeast Asia that has invaded the Mediterranean and several Caribbean islands (https://www.biolib.cz/en/taxon/id244539/). The females lay their eggs under the bark of the host tree (in Israel, mostly Ficus) and the developing larvae burrow into the stem of the tree, feeding on its sap and tissues. The emerging adults then burrow their way out of the tree stem and disperse to other trees by flight (Palaniswami et al., 1977). A previous study showed that these insects are capable of prolonged flight in flight-mills, lasting up to an hour and covering a distance of several km during a single flight (Brown et al., 2017). Interestingly, the beetles in the study that had a smaller body mass (~2 gr) tended to fly longer durations in flight-mills compared to larger individuals (with body mass ~4-6 gr). The flight speed of the smaller beetles was lower but they flew distances that were ~50% longer than larger beetles before reaching exhaustion (Brown et al., 2017). Intrigued by the prolonged flight capabilities of these beetles, we sought to identify their aerodynamic efficiency and power output during forward flight.
Specifically, we evaluated whether prolonged flight at the reported flight-mill speeds (2-2.5 ms -1 ) Downloaded from https://academic.oup.com/iob/advance-article/doi/10.1093/iob/obaa026/5902830 by guest on 18 September 2020 relies on high aerodynamic efficiency. To do so we flew the beetles tethered in a wind tunnel and estimated the power output using the observed kinematics and a quasi-steady blade-element model. We then compared the results with estimates derived from PIV measurements in the unsteady wake behind the flying insect.

Insects
Beetles for the study came from a population maintained at Tel Aviv University. Field-caught males and females were placed together in a large metal cage containing a small fig tree. The laid eggs were collected and incubated at 26° C in the dark. The hatched larvae were housed individually in plastic vials filled with moist peat. They were fed with an artificial agar-based medium (Surulivelu et al., 1977)

Wind tunnel
Experiments were carried out in the insect-flight wind-tunnel at Tel Aviv University. The threemeter-long wind-tunnel has a rectangular (0.2 × 0.21 m) cross-section and enables testing at wind speeds ranging from 0.11 to 7 ms -1 . The working section of the tunnel is 0.5 m long and equipped with transparent walls for visualizing the insect and flow. Wind speed is controlled by adjusting the voltage of a DC blower (RadiCal, ebm-papst, Germany) pulling air out at the downwind end of the wind-tunnel. The air passes through a series of 1/8 inch honeycombs (PAMG -XR1 5052, Plascore, US) to reach a uniform straight airflow. The wind speed throughout the working section was calibrated and verified before each experiment using a hotwire anemometer (SwemaAir40, Swema, Sweden). A small hole in the side of the working section allows a horizontal stainless steel rod (diameter 3 mm) to reach the center of the working section. At its tip is a vertical tether to which the insect is mounted using a drop of hot glue applied to the dorsal prothorax. Although tethered, the beetles were able to actively adjust their body angle by moving the head, meso-, meta-thorax and abdomen relative to the fixed prothorax.
During flight, the beetles tilted (pitch) their longitudinal axis by 24° ± 7° relative to the oncoming wind. This body angle is significantly lower (t-test, p<0.001) than the body angle when the wind was turned on but the beetles were not flying (39° ± 8°).

Direct force measurements
The free end of the horizontal rod holding the tether passed through the wind-tunnel wall into two low friction ball-bearings that allowed it to rotate about its length. The free end terminated in a 90° lever arm that converted the rotational motion of the rod to vertical force applied on a force transducer (0.5N, accuracy ± 0.05%, model ULC, Interface, USA,). The lever arm magnified the force applied by the tethered beetle through a mechanical advantage of 1:21. Both the force transducer and ball-bearings support were attached to a Perspex plate mounted on top of a second force-transducer (3N, accuracy ± 0.015%, Zimech LB9H) that measured the weight of the entire tethering setup. In this arrangement of force transducers and levers, the horizontal (up/down-wind) forces generated by the insect are measured with the first force transducer and the vertical forces (lift/weight) are measured with the second one. The mechanical advantage of the lever arm and the higher sensitivity of the first force transducer result in the horizontal forces having a negligible effect on the measurement of the vertical force and vice versa. Measurement data from both force transducers were continuously recorded at 100 Hz into the hard drive of a desktop computer via an analog-to-digital interface. The measured data were low-pass filtered with a cut-off frequency of 10 Hz and averaged for the duration of the measurement (5-6 flapping cycles =~0.2 sec). The entire force measurement system was calibrated and tested by placing measured weights at the tip of the tether (static) and by measuring the drag of a smooth plastic ball (diameter 4 cm) at a series of wind speeds and comparing the drag coefficient to published data on spheres. This yielded an error estimate of ± 4.1% for the horizontal forces and 0.7% for the vertical forces. For each beetle, forces in the wind-tunnel were measured three times: 1) with no wind and the beetle not flapping (i.e., the vertical force transducer measured the weight of the beetle and the horizontal force transducer measurement was null; 2) with a wind speed of 2.5 m s -1 and the beetle not flapping (body weight minus body lift measured in the vertical axis and parasite drag in the horizontal axis); and 3) with the wind on and the beetle flapping. We allowed the beetle to fly for two minutes to warm up and reach steady flight and then moved on to record the force, video and PIV data (see below) from the flying beetle.
By subtracting force measurement 1 from 2 we obtained the body lift and parasite drag for the given flight speed. Subtracting measurement 2 from measurement 3 gave the net horizontal and vertical thrust generated by the wings of the flying insect. The measurement of net horizontal thrust represents an underestimate of the true horizontal thrust provided by the beetle, because of the inability to separate drag from the measurement of steady-state thrust (Pennycuick, 1969;Rayner, 1999). In other words, the forces measured on the tether during flight are the reaction forces from the subtraction of forward thrust and drag. Since our force transducer only measures the difference between the two, the actual thrust must be higher than the measured force in order for the thrust generated by the wings to overcome the profile drag.

High-speed videography
Two high-speed cameras (Fastcam SA3_120K, Photron, USA) equipped with 50 mm lenses (Nikkor, Nikon, USA) were mounted perpendicular to the working section walls and captured the flapping motion of the insect inside the wind-tunnel at 2000 fps. We used infra-red back light illumination, which enabled us to carry out the experiments under normal visible light levels.
The cameras were calibrated using direct linear transformations (Hedrick, 2008) allowing us to reconstruct the 3D flapping from three visible landmarks on the right wing (Fig. 1). The 3D position data of wing landmarks were converted to three time-varying flapping angles relative to the body and an average stroke plane, as described in Meresman and Ribak, (2017) and Supplementary Material S1. In the experiments described here, we first used the high-speed cameras to analyze the wing flapping kinematics of ten beetles flying at 2.5 ms -1 . We then performed wake flow measurements using 2D-PIV (described below) on another six beetles at the same flight speed and same room temperature conditions.

Particle image velocimetry
PIV was utilized to measure the flow-field behind the tethered beetle during flapping flight. The PIV field of view captured the streamwise-normal plane, thereby measuring the streamwise and normal velocities. A double-head Nd:YAG laser was used to illuminate the flow field. The laser wavelength was 532 nm, emitting 200mJ/pulse at 15 Hz (Evergreen, Quantel). The light sheet was formed using a 15 mm cylindrical lens followed by a 500 mm spherical lens. At mid downstroke, the light sheet intersected the beetle's left wing at 2/3 the wing length. The time between the laser pulses was set to 200 sec to allow sufficient pixel displacement per interrogation window. The camera used was 8MP (3320x2476 pixle 2 ) CCD double exposure camera operating at 4 Hz with a 12-bit dynamic range capturing a field of view of 220x141 mm 2 .
The particles used to seed the flow were olive oil droplets with an average dimeter of 1 m, generated using a Laskin nozzle. The PIV analysis was performed using an interrogation window of 64 × 64 pixle 2 with a 50% overlap (Insight 4G, TSI). The obtained velocities were post processed, including local vector validation based on a median test from the surrounding 5 × 5 vectors. Holes in the vector maps were filled with the local mean value calculated from surrounding measured vectors. The number of erroneous vectors per map was less than 5%. A total of 3,000 vectors were yielded per maps and about 300 vector maps were acquired per beetle (Table 1) to enable statistical convergence of the data. The spanwise vorticity was calculated using least-square fitting to the acquired 2D velocity fields (Raffel et al., 2007).

Data analysis: sectional drag based on PIV
We measured the mean drag from all six beetles using PIV maps. In three of the beetles (beetles 4-6 in Table 1) we also explored the variation of the flow field within the flapping cycle. To do so we first divided the flapping cycle into eight distinct phases as described in Table 2. We then binned the PIV maps according to phase and calculated averages per phase for each of the three beetles.
The sectional drag was estimated for the beetles during flapping flight from the momentum deficit developed at the wake behind the wing as in Ben-Gida et al. (2013) and Nafi et al. (2020), with the principle originally formulated by Goett (1939). The drag force is derived from the momentum equation. The streamwise x-component of the momentum equation can be re-written to express the drag force per unit span: where is the velocity vector in Cartesian coordinate system, p is the pressure and ρ is the density. Defining S to be sufficiently far from the wing, where the pressure is assumed to be p ∞ (constant), and after substituting the continuity equation into Eq. 1, the drag force per unit span can be written as follows (Ben-Gida et al., 2013): , is referred to as the velocity deficit drag (Goett, 1939), while the unsteady drag component D 1 ' is added due to the flapping motion. The steady drag term can be obtained from the near wake velocity field. We assume that the unsteady drag, given that the wing flapping motion is close to harmonic (Ellington, 1999;Weis-Fogh, 1973), is approximately zero over a wingbeat cycle. Here, U ∞ is the mean undisturbed streamwise velocity in the wind tunnel, and h and l are the vertical and horizontal extent of the computed velocity field in the wake, respectively. The measured sectional drag (based on the momentum deficit) is proportional to the sectional profile drag (Goett, 1939). To obtain the drag per beetle we multiplied the sectional drag by the wing span.
Data analysis: energy expenditure from flapping kinematics: We calculated the power required for moving the wing through air at the observed kinematics using a quasi-steady blade-element model (Sane and Dickinson, 2002) adjusted for forward flight (Dickson and Dickinson 2004). We found the instantaneous aerodynamic force generated during flapping as the sum of three forces: where F t , F r and F a are the translational, rotational and acceleration (added mass) forces, respectively. They are calculated as: [Eq. 4] where and are the wing area and length, respectively, U is the wingtip speed, and are the non-dimensional first and second moment of wing area. and are the non-dimensional lift and drag coefficients defined in the Supplementary Material S1 and is the nondimensional tip velocity ratio (Dickson and Dickenson, 2004) that relates forward flight speed to wing flapping speed at the wingtip.
The rotational forces were calculated as in Sane and Dickinson (2002), assuming the axis of rotation is located at 0.25 the wing chord from the wing's leading edge: [Eq. 5] where and are the mean and local (radial) non-dimensional chords, respectively, and is the instantaneous angle of attack.
The instantaneous force due to the acceleration of the added mass (Sane and Dickinson, 2002) is: where is the instantaneous wing flapping angle and the single and double dots over angle symbols in Eq. 5 and 6 denote first and second time derivatives.
From the resultant instantaneous force, F aero , (Eq. 3), which is assumed to act perpendicular to the wing's chord, we calculated the mechanical power from the scalar product of the force and wing velocity vector.
[Eq. 7] Here, we assumed that the resultant forces act at a point located at 0.56 the distance from the hinge to the tip of the wing, and used the instantaneous velocity of the wing (relative to air) at this point. The instantaneous power was converted to work per cycle by integrating the power over the flapping cycle period. The contribution of the wing mass to inertial work and power was not included because 1) this extra work can be recycled in oscillatory flapping using elastic energy storage (Dudley and Ellington, 1990a;Ellington, 1985); and 2) not including the inertia of wing mass left us with the aerodynamic power that is comparable to the changes in kinetic energy of the flow, derived from the PIV measurements (see below). Power values (P) were divided by the flight-muscle mass of the beetles to account for the intraspecific variation in body mass in our beetles. The flight muscle mass was derived from body mass using the allometric equation in Brown et al. (2017). Thus, power is hereafter reported as muscle mass-specific power, denoted by, P*.

Data analysis: estimation of induced power and efficiency
In order to evaluate the calculated power expenditure of the beetles against a standard baseline, we compared the power estimations derived in our experiments with the theoretical minimal power output required for forward flight. This power output can be formulated as: where is the parasite power from the drag on the body. It was found by multiplying the parasite drag measured by the force transducer by the wind speed in the wind tunnel.
is the profile drag that results from the drag on the flapping wings. Pennycuick (1969) proposed that: with for a pigeon. Tucker (1973) (Warfvinge et al., 2017). For lack of better data on insects we calculated for 3.5. The induced power, , is required for weight support in the air. Pennycuick (1969) used the actuator disc principle (von Miss, 1945) to formulate the equation for the theoretical (minimal) induced power, relevant for medium to high flight speeds: where is the forward flight speed and W is the weight of the flyer. Since our tethered beetles supported only a part of their body weight during flight (see Results), we used the vertical force measured by the force transducer instead of total body weight.
The aerodynamic propulsion efficiency ( can be calculated as the ratio of the power output ( and the power invested moving the wing through air (P): [Eq.

11]
For P in the denominator, we used values from our two independent approaches: Power found from the blade-element analysis and the power based on the change in kinetic energy in the wake (see below). The efficiency is if 100% of the power spent moving the wings is converted to thrust for moving forward at the flight speed while supporting the body weight in air. It is less than 1 and greater than zero if some of that energy is dissipated in the process.

Kinetic energy of the surrounding fluid
We utilized a control volume (CV) approach to estimate energy conservation and transport within the wind tunnel. The control volume was defined as an imaginary box that spans from the incoming flow towards the beetle and continues towards the wake, where a distinct momentum deficit is observed (by the PIV). Its size was defined as the volume of the wind tunnel section equal to twice the length of the PIV field of view multiplied by the vertical length of the PIV field of view and by wind tunnel width. The tethered insect was held at the center of the control volume. Since the beetle was tethered, we can assume that there was no change in potential energy within the control volume since there was no mass change in altitude and the air mass moving can be assumed to be negligible given its density. Therefore, any changes in kinetic energy within the control volume are due primarily to the interaction of the beetle's wings and body with the oncoming flow.
The kinetic energy at the inlet side of the control volume (upstream section) was constant because the mean velocity profile at the upstream section (prior to interaction with the beetle) is uniform: . The kinetic energy of the flow at the outlet side of the control volume is a function of the momentum deficit; , where h is the wind tunnel height. It is noteworthy to mention that the kinetic energy at the outlet is based on the PIV measurement performed at the 2/3 wingspan plane. Since the kinetic energy is an integral property of the CV, we are using the mean velocity profiles at the wake to estimates its value. It is acknowledged that the flow at the wake is turbulent and 3D by nature, hence, the instantaneous velocity profiles across the wing span will be different. However, the average profiles at the midspan region encompass the momentum arriving from the wing, mixing at the wake, interacting and dissipating. The resultant mean profiles comprise of all these phenomena; thus, we suggest that at the mid span region, the mean profile will appear somewhat similar. The difference between the two: input and output, assuming steady-state conditions, is equal to the amount of work done within the control volume that relates to the aerodynamic work of the flapping wings and drag of the body. Hence, the energy expenditure of the beetle during flight can be evaluated from a comparison between the energy exerted on the flow using the PIV measurements and the work of the flapping wing and from the body (parasite) drag.
Data are reported below as mean ± Standard deviation.

Results:
The results are divided into two sub-sections; flapping kinematics and force-transducer data measured directly from the insect motion; and wake flow data measured from the PIV. All three measurements (kinematics, force and PIV) were independent in the sense that there was no interaction ('cross-talk) between them.

Direct measurements (beetle):
During flight, the beetle's stroke plane was tilted from the horizontal plane by β = 52° ± 6.2° (n=10). Figure 2 shows the three time-varying flapping angles relative to this stroke plane. The mean flapping amplitude was 161° ± 13.14° and the mean flapping frequency was 35 ± 1.9 Hz, resulting in a mean advance-ratio (forward flight speed/mean wingtip speed) of 0.30 ± 0.03. The upstroke and downstroke each comprised half of the flapping cycle (Table 2).

Direct force measurements
With the tethered beetles not flying and the wind speed in the wind tunnel set to 2.5 ms -1 the force transducer measured the parasite drag of the beetles as 0.6 ± 0.35 mN. When the beetles were flying at the same wind speed, they applied a net forward force of 2.3 ± 1.86 mN on the tether, implying that their horizontal thrust (at least 3.0 ± 1.6 mN) was in excess of the drag on the body and wings. In contrast, on the vertical axis, the beetles supported only 55% (± 19%) of their body weight on average (Fig. 3). The time-averaged resultant force measured with both force transducers ( , where and are the vertical and horizontal forces, respectively) showed a good correlation (r = 0.87) to the resultant force calculated from the blade-element analysis and averaged over the flapping cycle (Fig. 4). The estimates from the blade-element were on average 18% ± 25.8% lower than the measured forces.

Estimated aerodynamic work
Based on the blade-element analysis and the flight muscle mass from the allometric equation in Brown et al. (2017), the mean muscle mass-specific power of the flying beetles (excluding the energy spent on accelerating the mass of the wing) was 86.9 ± 18.76 W kg -1 muscle (Fig. 3).
Substituting the estimated power in Eq. 11 yielded a propulsion efficiency that varied substantially among the ten beetles. For 2 the average propulsion efficiency was 0.42 ± 0.21 (Fig. 5). The estimated work was on average 2.4 ± 1.06 mJ for a single flapping cycle.

Wake characteristics
The instantaneous spanwise vorticity maps in the wake of the beetle measured using PIV are shown in Figure 6. These maps are representative ones, demonstrating the flow evolution at the wake region that results from the flapping wing motion. One can observe the unorganized pattern of the vorticity fields, indicating strong unsteady activity due to the flapping motion and phasedependent instantaneous regions associated with the downwash and upwash of fluid in the normal direction. Figure 6a corresponds to the early stages of the downstroke (DD in Table 2

Drag trends over the wingbeat cycle using PIV
The sectional drag, estimated based on the momentum deficit in the wake, corresponded to the wingbeat phase. Figure  conditions, when the body is stationary (here, tethered), one may assume that these will be equal.
The mean (steady) sectional drag as estimated by the PIV and multiplied by the wing span gave an estimate of the total drag of the flying beetle. With few exceptions, the drag values were in agreement with the horizontal thrust measured by the force transducers (Table 1). The discrepancies were on average ~30%.  Table 3 reflects the averaged kinetic energy over the wingbeat cycle. We can observe a trend similar to the drag over the wingbeat phase where during both the downstroke and upstroke, the change in kinetic energy values fluctuate between 1.2 and 1.4  10 -3 Joules. However, these variations are not large, mainly due to the energy being a scalar-integral property; i.e., the velocity profiles are averaged twice to allow a 1-number representation of the energy. Converting to power (muscle mass-specific power = 64.9 ± 7.8 W kg -1 muscle) and substituting in Equation 11, the calculated propulsion efficiency is between 0.57 ± 0.08 and 0.99 ± 0.13 for k =1 and 2.5, respectively. Calculation of propulsion efficiency for k>2.5 yielded values exceeding 1.0. For k=2 (Pennycuick, 1969) the propulsion efficiency becomes 0.85 ± 0.11.

Discussion
Forward flight speed Brown et al. (2017) reported flight speeds for B. rufomaculata flying in a flight-mill, in which beetles with a body mass similar to those used here (2 gr), had a mean flight speed of 2.0 ms -1 ; while larger beetles with a mean body mass of 4.5 gr flew at a mean flight speed of 2.5 ms -1 . Our beetles were flown at a wind speed of 2.5 ms -1 and generated a net horizontal force on the tether, implying that they were generating horizontal thrust in excess of the drag on the body and wings.
Therefore, had they not been tethered, they would have accelerated to higher flight speeds. Thus, the 2.0 and 2.5 ms -1 flight speeds observed by Brown et al. (2017)  are not available for B. rufomaculata, preventing us from differentiating between these two options. Failure to provide enough vertical force to support body weight is common in tethered flight preparations (Jensen, 1956;Ribak et al., 2017;Shkarayev and Kumar, 2015). The force transducers were specifically added to our research setup to account for such tethering biases in our calculations of propulsion power and efficiency. Had we used the weight of the beetles in Eq.
10 instead of the vertical force they actually produced, the resulting propulsion efficiency would have been erroneously high. Constraints due to tethering limit our ability to directly relate our results to free flying-beetles. Nevertheless, it would not have been possible to measure forces directly and perform accurate flow measurements over time had the beetles in our study were not tethered. Despite its limitations, the tethered flight setup described here does provide an efficient means to relate wing flapping kinematics with aerodynamics enabling the measurement of profile drag and flight efficiency.

Wake-wing interaction
The interaction between the beetle and the wake flow is described here qualitatively via the instantaneous wake flow vorticity-velocity maps and quantitatively via the beetle's aerodynamic loads. The flow maps (Fig. 6) show fluctuations in the flow that are associated with the flapping cycle. Since the PIV system sampling rate was lower than the wingbeat frequency, these flow maps are not consecutive ones, although they do correspond to the different wingbeat phases.
The fluctuations are based on applying Reynolds decomposition to the velocity field, in which we subtract the ensembled average for each experiment set from the instantaneous maps. The resulting fluctuating velocity field represents the turbulence portion of the flow (Tennekes and Lumley, 1972). Johansson et al. (2012) noted that the wake of tethered beetles is rather complex, partly due to the presence of the elytra interfering with the flow over the hindwing. Here, we identified whether the wing generated upwash or downwash motion based on the direction of the velocity vectors coupled with concentrated regions of vorticity . This identification is somewhat subjective but provides a qualitative link between the wing motion and the momentum change in the wake. The downstroke was associated with downwash and advection of the concentrated vorticity downstream. The upstroke was associated with upwash, suggesting that it was also active in generating aerodynamic forces. The dorsal and ventral stroke reversals (mainly the VSR) showed transitions associated with multiple concentrated vorticity regions (positive and negative). The appearance of multiple rapid momentum changes in individual fluctuating maps resulted from the high wingbeat frequency that caused the flow to change direction over a short space. It should be noted that there was some delay between the wing motion and the arrival of change of momentum into the PIV field of view (which was located a few chord lengths downstream). This can be clearly observed for the VSR transition, where the wing moves from downstroke to upstroke.
The wing motion generated momentum changes in the wake velocity field -both streamwise and normal. These momentum changes were identified by the appearance of concentrated regions of vorticity within the field of view, in an unorganized manner. As a consequence, the aerodynamic loads exerted by the beetle also changed and were manifested by this interaction. The sectional drag values changed over the cycle as a function of the phases, fluctuating between 0.04 to 0.08 Nm -1 . The variations that followed the wingbeat kinematic phases demonstrate that the sectional drag peaked slightly after and before the dorsal and ventral transition phases (DSR and VSR), respectively. In addition, for two out of the three beetles, the drag during the upstroke was slightly higher than during the downstroke. The peak in drag during the stroke reversals probably resulted from the wing rotation. Rapid wing rotation at stroke reversal is associated with peaks in the aerodynamic forces, despite the flapping wing having low translation speeds (Dickinson, 1994;Dickinson et al., 1999;Sane and Dickinson, 2002). The forward speed of the body during wing supination has an effect similar to advanced wing rotation, a condition associated with peaks of drag at the VSR (Sane and Dickinson, 2002). At the DSR, as the wing pronates the forward speed of the body has the opposite effect (delayed rotation), which can lead to a large peak in instantaneous drag at the beginning of the downstroke (Sane and Dickinson, 2002).
These unsteady effects of wing rotation fit quite nicely with the bimodal shape of the sectional drag curve during the downstroke (Fig 7). The higher mean drag during the upstroke, compared to the downstroke, is probably due to the higher angle of attack during this phase (Fig. 3). The increase in angle of attack is due to the forward flight speed interacting with the tilted stroke plane, leading to a steeper wingtip path relative to air during the upstroke (Ellington, 1984).

Comparing the two calculation approaches
In the experiments and analyses that followed, we utilized two independent methods to measure the force and power requirements for forward flight. Both methods gave similar thrust-drag values and provided an estimate of the power and mechanical work values that were in the same order of magnitude and within the scope of insect flight muscle physiology. With respect to forces, the estimates from the blade-element analysis were comparable to the direct measurements with the force transducers (Fig. 4). The indirect drag measurements from the flow fields (obtained from the PIV data) were comparable to the direct horizontal force measurements from the force transducer. The discrepancies were about 30% on average. Further inspection revealed that in four of the experiment sets we obtained similar drag values using both techniques, while in the two other experiments, the values differed (Table 1). Therefore, one can conclude that estimating the drag using the momentum deficit at the wake region provides a good estimate of the aerodynamic forces exerted by the beetle during flight, which can also be used for free-flight experiments, assuming that the wake flow is measured (Ben-Gida et al., 2013).
The measured and calculated forces using the two techniques enabled us to estimate the propulsion efficiency of the beetle. Our calculation of propulsion efficiency was based on the minimal power output, which depends on the unknown parameter k in the calculation of profile power (Eq. 9). To account for the uncertainty in the latter we calculated for a range of k values ( Fig. 3), which indeed have a large effect on the estimated efficiency. We note that values for k>2.5 give unrealistically high efficiencies with efficiency (using either the blade-element approach or the flow field data) exceeding 1.0 in some beetles. The power output for k =2 (Pennycuick, 1969) gives a mean efficiency of 0.42, based on the blade-element. This value implies that 58% of the aerodynamic power invested is wasted during conversion to the minimal forward thrust that propels the body forward while keeping the insect in air. This efficiency is much lower than the 0.81 propulsion efficiency of swimming bottlenose dolphins (Fish, 1993) and 0.89 of Manta rays (Fish et al., 2016) which use lift-based propulsion but need to provide little body support; and higher than the 0.33 of swimming muskrats, which use drag-based paddling for surface swimming (Fish, 1984). Since the efficiency is calculated relative to the minimal power output, replacing the power output with a more realistic (higher) estimate by including additional morphological and kinematic parameters, would yield higher power output and propulsion efficiencies (Eq. 11). Specifically, the induced power found using Eq. 10 is the theoretical minimum for fixed wing flyers with elliptical (ideal) lift distribution over the span.
This calculation of induced power does not take into account the flapping amplitude, stroke plane angle, downwash speed and span efficiency (Henningsson and Bomphrey, 2011;Henningsson and Bomphrey, 2013), which would result in a higher induced power and, therefore, a higher estimate of propulsion efficiencies. Nevertheless, we suggest that comparing the power from the blade-element and flow fields to a standard minimum has a more meaningful interpretation, in pooling all of the above into the term "propulsion efficiency" and thus providing the means to compare very different flyers against a standardized baseline.
It is important to note that although the PIV and blade-element techniques calculate aerodynamic of free-flying hawkmoths using PIV and blade-element modelling, taking a slightly different approach to that presented here. They found an aerodynamic power expenditure of ~20 mW and up to ~30 mW, respectively, in two hawkmoths flying at 2.5 ms -1 . With a mean body mass of 2 gr and flight muscle varying between 21% and 34% of the hawkmoths' body mass, the muscle mass-specific power of the moths was between 29.4 and 47.6 W kg -1 muscle for the PIV measurement, and 44.1 -71.4 W kg -1 muscle for the blade-element, respectively. These values are somewhat lower than our 71 and 84 W kg -1 derived from the PIV and blade-element approaches.
The difference may correspond to the lower wingbeat frequency and wing loading of the moths as well as their more streamlined body. Using a blade-element approach, Dudley and Ellington (1990a) estimated body mass-specific power of forward flying bumblebees. Their results indicate 20-35 W kg -1 for flight at 2.5 ms -1 . With a flight-muscle ratio of 0.37 out of the total body mass in Bombus terrestris (Polidori and Nieves-Aldrey 2015), this gives a muscle mass-specific power of 54-95 W kg -1 muscle. These values are somewhat higher than ours but they include the inertial power of the wing mass under the assumption of perfect elastic energy storage; whereas our estimates include the inertial power of the added mass but not the inertial power of the mass of the wing itself. Ellington (1985) reported aerodynamic power values between 167-186 W kg -1 muscle for hovering honey bees, bumblebees and drone flies. Hence the muscle mass-specific power found here is within the range expected for flying insects. Considering the fact that weight support of our beetle was below the level required for free-flight, the submaximal power output of the flight muscle is to be expected. Submaximal muscle power output of the flight muscles might also delay fatigue during prolonged flight.
Our study highlights some of the sources of mechanical energy loss associated with the aerodynamics of forward flight in a large beetle. The ability of B. rufomaculata to undertake prolonged flights at a speed of ~50 body lengths per second underscores the remarkable dispersal capacity and endurance of these large insects.      Table 2 where the cycle starts with DSR (a) and ends with UD (h). b-d and f-h are the downstroke and upstroke, respectively, and e is the transition (VSR) between them. The maps are available at higher resolution in Supplementary Material S2 Figure 7 Sectional drag (D', drag per unit span) calculated from the PIV over one wingbeat cycle. The sectional drag was estimated and averaged for each wingbeat phase in each of beetles 4-6 (color coded lines). The horizontal axis are the different wingbeat phases as depicted in Table 2.