Olivier Macherey
email: macherey@lma.cnrs-mrs.fr

Sensitivity to temporal fine structure probed with pulse-spreading harmonic complexes

 provided evidence that NH subjects are sensitive to TFS at very high frequencies in stimulus conditions where spectral cues should not be present. This is important because it implies that phase locking may still provide useful information above 4000 Hz where it is commonly believed to be absent. Here, we examine sensitivity to TFS using different stimuli than the frequency-shifted complexes previously used. Pulse-spreading harmonic complexes (PSHCs) are spectrally-dense complexes with a low f0 (here 2 Hz). The phase of the components is set so that the envelope rate can be manipulated independently from the f0.

Here, we use a special case of PSHCs for which the TFS peaks were shifted relative to the envelope peaks from one envelope period to the next. For upward PSHCs (up-PSHC), the TFS peaks were progressively advanced, yielding a rising pitch percept across the stimulus, while for downward PSHCs (down-PSHC), the peaks were progressively delayed, yielding a falling pitch percept. In Experiment 1, subjects heard up-PSHC, down-PSHC and regular inharmonic complexes bandpassfiltered with a lower cut-off (Fc) of 2000 Hz and were asked to identify whether the pitch was rising, falling or flat in a one interval, three alternative forced choice task. The envelope rate was 98, 200 or 450 pps. For the two highest rates, most subjects correctly identified the direction of the pitch change while for the lowest rate, the stimuli could not be discriminated. This shows that, for a given Fc, there is an envelope rate below which the auditory system is insensitive to changes in TFS. The lower rate limit of TFS sensitivity was measured in Experiment 2 using an adaptive three interval, two alternative forced choice task where subjects heard two down-PSHC and one up-PSHC and had to pick the latter. The stimuli were bandpass filtered in seven frequency regions with Fc varying from 250 to 11200 Hz and were presented in the presence of threshold-equalising noise. The lower limit significantly increased from 42 pps to 691 pps with increases in Fc. The ratio between Fc and the envelope rate at threshold also increased from 6 to 18 when Fc increased from 250 to 2000 Hz, then remained flat with further increases in Fc or slightly decreased at the highest Fc. In the lower part of the spectrum (<2000 Hz), this lower limit resembles the lower limit of melodic pitch. The results in the upper part of the spectrum suggest that TFS cues are still available in high frequency regions (up to 10 kHz) when the rank of the lowest harmonic present in the passband is 18, and all harmonics are presumably unresolved.

Introduction

The spectral decomposition performed by the cochlea resembles that achieved by a bank of overlapping bandpass filters: in response to a sound, the vibration pattern of a particular point on the basilar membrane corresponds to a bandpass-filtered version of the original sound waveform. This filtered waveform can be decomposed into a slowly varying envelope superimposed on a carrier, referred to as the temporal fine structure (TFS), which oscillates at a rate close to the characteristic frequency of this particular cochlear location. Auditory nerve fibers can, to some extent, phase-lock to this TFS and there is evidence that this temporal neural information is used by the binaural system for frequencies up to about 1400 Hz [START_REF] Zwislocki | Just Noticeable differences in dichotic phase[END_REF]. The exact limit of phase locking in humans and whether or not it is important for pitch perception remain, however, a matter of debate. Most animal species exhibit a steep decrease in phase locking for frequencies above about 1000 or 2000 Hz (Koppl, 1997). The finding that pure tone frequency discrimination in humans also deteriorates for frequencies above 2000 Hz has often been cited as evidence that frequency discrimination depends on phase locking information (e.g., [START_REF] Moore | Frequency difference limens at high frequencies: evidence for a transition from a temporal to a place code[END_REF]. However, [START_REF] Oxenham | Pitch perception beyond the traditional existence region of pitch[END_REF] reported that harmonic complex tones only containing components above 4000 Hz can still convey a clear, musical pitch percept. They concluded either that phase locking information is not needed for complex pitch perception or that it is still usable by humans at these very high frequencies. [START_REF] Moore | Sensitivity of the human auditory system to temporal fine structure at high frequencies[END_REF] investigated sensitivity to TFS at very high frequencies (above 8 kHz) by asking their subjects to discriminate between bandpass filtered harmonic and inharmonic complex tones. In each trial, the inharmonic complex contained the same components as the harmonic complex except the components were shifted by a fixed amount in Hz. These two stimulus waveforms have identical envelopes but different TFS. [START_REF] Moore | Sensitivity of the human auditory system to temporal fine structure at high frequencies[END_REF] showed that subjects could discriminate between these two stimuli when the lowest harmonic present in the passband was the 12 th . According to usual definitions of resolvability [START_REF] Moore | Resolvability of components in complex tones and implications for theories of pitch perception[END_REF], this would imply that the components were unresolved and that the subjects could not use differences in excitation patterns to perform the task. They, therefore, concluded that TFS sensitivity was mediated by phase locking information which in turn implied that phase locking cues could still be used by their subjects at frequencies higher than 8 kHz. This is controversial because as previously mentioned, animal experiments show poor phase locking at these frequencies and there are currently no human single-unit data available in the literature. Note, however, that significant phase locking up to 7 kHz in chinchillas were obtained by Kale and Heinz when averaging a lot of data to lower the noise floor of their recordings, thereby suggesting there may still be usable temporal information even at these high frequencies [START_REF] Verschooten | The upper frequency limit for the use of phase locking to code temporal fine structure in humans: A compilation of viewpoints[END_REF]. A recent opinion paper on the use of phase locking cues in humans illustrates this controversy: estimates of the highest frequency at which phase-locking information could be used by humans ranged across "experts" in this field from 1.5 kHz to 10 kHz [START_REF] Verschooten | The upper frequency limit for the use of phase locking to code temporal fine structure in humans: A compilation of viewpoints[END_REF].

One draw-back of using frequency-shifted inharmonic complex tones to investigate TFS sensitivity is that these stimuli produce ambiguous pitch percepts and that subjects were not explicitly asked to compare the pitches of inharmonic and harmonic complex tones; any audible difference between the harmonic and inharmonic tones could be used as a cue. Furthermore, harmonic and inharmonic complexes produce different envelopes at the output of a given auditory filter which can potentially provide another cue to discriminate them. This led [START_REF] Moore | Sensitivity of the human auditory system to temporal fine structure at high frequencies[END_REF] to randomize the starting phase of their components, which in turn affected the shape of the stimulus envelope. Both this pitch ambiguity and the phase randomization may have an impact on measures of TFS sensitivity and might impair subjects' performance.

Here, we reexamine the TFS sensitivity of normal-hearing subjects using pulse spreading harmonic complexes (PSHCs; [START_REF] Hilkhuysen | Optimizing pulse-spreading harmonic complexes to minimize intrinsic modulations after auditory filtering[END_REF]. PSHCs are harmonic complexes whose phase relationship is tailored to manipulate the envelope rate independently from the f0.

In Experiment 1, we show that it is possible to produce two PSHCs that have the same long-term spectrum, and identical envelopes but are perceived as either rising or falling in pitch across their duration, thereby avoiding the pitch ambiguity of inharmonic complexes. Acoustically, these stimuli only differ in their TFS. In Experiment 2, we use these stimuli to measure the lower limit of TFS sensitivity in different frequency regions and discuss the putative mechanisms underlying this limit.

Experiment 1: Pitch contour identification

Methods

Subjects and Equipment

Six normal-hearing subjects with ages ranging from 26 to 57 years (mean of 38.7) took part.

Subject S4 was the author. They were seated in a double-walled sound proof booth in front of a computer screen. Stimuli were presented over an Etymotic ER2 earphone positioned inside one of their ears. The earphone was connected to a RME BabyFace Pro sound card. The experiment took place in the presence of a background noise that was played from a Zoom H4N recorder connected to the sound card.

Generation of PSHCs

The generation of PSHCs involves dividing the harmonics of a complex into several groups (called sub-complexes) and summing them separately before adding them back together with particular delays so that the overall envelope repetition rate remains regular across the stimulus duration.

Following [START_REF] Hilkhuysen | Optimizing pulse-spreading harmonic complexes to minimize intrinsic modulations after auditory filtering[END_REF], a PSHC is defined by

𝑠(𝑡/𝑘) = ∑ sin (2𝜋𝑓 0 𝑖𝑡 + 𝜑(𝑖, 𝑘)) 𝑁 𝑖=𝑀 (1)
where f0 is the fundamental frequency, M is the lowest harmonic present, i is the harmonic number, N is the highest harmonic present, and t is the time. The starting phase of harmonic i is given by

𝜑(𝑖, 𝑘) = 2𝜋 𝑖 𝑘² 𝑟 𝑗 ( 2
)
where k is the PSHC order which determines the envelope repetition rate (equal to k²f0) and j is defined as

𝑗 = (𝑖mod𝑘) + 1 (3)
The mod operator indicates the remainder after division of i by k. Equation (3) assigns the i th harmonic to the j th sub-complex so that there are in total k sub-complexes. In our previous publications [START_REF] Hilkhuysen | Optimizing pulse-spreading harmonic complexes to minimize intrinsic modulations after auditory filtering[END_REF]Mesnildrey et al., 2015), rj was a random element without replacement from {1, 2, …, k} whose role was to spread the sub-complexes evenly across the stimulus period, presenting them in random order. Here, we consider two special cases of PSHCs where the sub-complexes are presented either in descending order (down-PSHC) or in ascending order (up-PSHC). For down-PSHC, rj is given by

𝑟 𝑗 = 𝑗 (4)
while for up-PSHC, it is given by

𝑟 𝑗 = 𝑘 -𝑗 + 1 (5)

Stimuli

Stimuli consisted of up-PSHCs, down-PSHCs and inharmonic (INHARM) complexes synthesized in Matlab with a sampling rate of 96 kHz. Three versions of each of these stimuli were generated with envelope rates equal to 98 pulses per second (pps), 200 pps and 450 pps.

For the INHARM stimuli, the frequency components were randomly drawn with the constraint that the frequency separation between consecutive components was constant and equal to the envelope rate for each condition. For up-and down-PSHCs, the fundamental frequency was 2 Hz and the PSHC order k was 7, 10 and 15 for the three envelope rates, respectively.

All stimuli were bandpass filtered using a zero-phase 6 th -order Butterworth filter (36 dB/oct. attenuation slopes) with -3-dB cut-off frequencies of 2000 and 2540 Hz, yielding a 2-ERBN wide passband [START_REF] Glasberg | Derivation of auditory filter shapes from notched-noise data[END_REF]. The stimuli had a total duration of 500 ms (corresponding to a full period of the PSHCs) and were shaped with 20-ms raised-cosine onset and offset ramps.

The PSHCs display gliding spectral ripples, as shown in Figure 1. For each presentation of the upand down-PSHCs, the ripple starting phase varied randomly so that subjects could not identify the stimulus solely based on its onset or offset. The stimuli had an overall level 40 dB higher than the threshold in quiet of a pure tone with frequency corresponding to the middle of the passband.

A continuous background of threshold-equalising noise (TEN) was presented at a level of 25 dB/ ERBN (same SNR as [START_REF] Moore | Sensitivity of the human auditory system to temporal fine structure at high frequencies[END_REF]).

The time waveforms and spectrograms of the three stimuli are shown in Figure 1 for an envelope rate of 450 pps. The blue curves in the top panels represent the time waveforms while the red curves show the Hilbert envelopes. At a given rate, the envelopes of the down-and up-PSHCs are identical but are slightly narrower and have a larger peak amplitude than the envelope of the INHARM complex (when compared at equal RMS). The spectrograms in the bottom panels show 2-seconds of the three stimuli. The down-and up-PSHCs exhibit linear spectro-temporal ripples.

The ripple density is determined by the envelope rate which corresponds to the frequency separation between adjacent ripples, while the ripple velocity is determined by the f0. Given that the PSHCs contain only harmonics of 2 Hz, they are both periodic, exactly repeating every 500 ms.

The shallow-slope filters applied to the stimuli make them similar to continuous Shepard tones (referred to as "Risset glissando") except that in our case the frequency separation between adjacent ripples is equal to a constant value in Hz whereas the components in Shepard tones are uniformly spaced on an octave scale [START_REF] Shepard | Circularity in judgments of relative pitch[END_REF][START_REF] Risset | Pitch control and pitch paradoxes demonstrated with computer-synthesized sounds[END_REF].

Procedure

The whole experiment was completed in one session lasting roughly an hour.

To determine the appropriate level of the complexes, the detection threshold was measured for a 2254-Hz (geometric center of the passband), 500-ms pure tone using a three interval, three alternative forced choice adaptive procedure with a 2-down, 1-up rule (c.f. [START_REF] Moore | Frequency difference limens at high frequencies: evidence for a transition from a temporal to a place code[END_REF]. The signal level was adapted in steps of 8 dB until the second reversal, then in steps of 4 dB until the fourth reversal and finally in steps of 2 dB until the tenth reversal when the procedure stopped. The threshold was calculated as the mean signal level at the last six reversals. Two threshold estimates per subject were collected and their threshold was assumed to be the mean of these two estimates. As previously mentioned, the overall level of the complexes and the level of the background TEN used in the main part of the experiment were set at 40 dB and 25 dB/ERBN above this pure tone threshold, respectively.

Experiment 1a used an identification task. On each trial, subjects heard one of the 9 stimuli (3 stimulus types x 3 envelope rates) and had to indicate whether they perceived the pitch as rising, falling or staying constant across the 500-ms stimulus. Feedback was not provided. All stimuli and envelope rates were presented in random order in blocks of 90 trials. Each subject performed a total of 5 blocks, yielding 50 trials per stimulus in total.

In Experiment 1b, the same subjects took part in a three-interval two-alternative, "odd-man-out" task to investigate whether they could discriminate between the up-PSHC, down-PSHC and INHARM stimuli at the lowest envelope rate (98 pps). They were presented with a reference stimulus in interval 1 followed by the same reference and a target stimulus randomly presented in intervals 2 and 3. There were three different reference-target trials (up-PSHC vs. down-PSHC, up-PSHC vs. INHARM and down-PSHC vs. INHARM). The trials were presented in blocks of 105, consisting of 75 trials at 98 pps (25 trials for each condition) and 30 easier "catch" trials (5 trials per condition for the two higher envelope rates). All subjects performed two blocks in total.

Feedback was provided.

Results

The top panel in Figure 2 shows the percent correct scores obtained in Expt. 1a for the six subjects and the three envelope rates. Five out of six and six out of six subjects had scores higher than 92%

for envelope rates of 200 and 450 pps, respectively, showing that their perception was, at these rates, consistent with expectations based on the time-frequency representations of the stimuli (Figure 1). For the lowest rate (98 pps), however, all subjects performed at chance. Further inspection of the data revealed that they perceived the stimuli to remain constant in pitch in 90% of trials (i.e., for the three stimulus types), indicating that both the down-and up-PSHCs were also perceived as constant in pitch.

The results of Expt. 1b (Figure 2, bottom panel) show that only S5 performed significantly above chance for discriminating the inharmonic complex from each of the two PSHCs (up vs. flat and down vs. flat conditions). However, none of the subjects could discriminate between the up-and down-PSHCs. This shows that at 98 pps, subjects were not sensitive to TFS. It is possible that S5 used remaining envelope differences between the inharmonic tone and the PSHCs to perform the task. The observation that he could not discriminate between the two PSHCs which had identical envelopes would argue in favour of this explanation.

Discussion

The results of Experiment 1 show that, for a given passband, there is an envelope rate below which down-and up-PSHCs cannot be discriminated and are not being perceived as changing in pitch. Below we describe the cues that the subjects may have used to perform the pitch contour identification task of Expt. 1a and discuss the possible reasons why their performance got poorer at the lowest envelope rate. These cues are excitation pattern cues and phase locking cues.

The top panels of Figure 3 show simulations of basilar membrane vibrations using a Gammatone filterbank [START_REF] Slaney | Auditory toolbox[END_REF] in response to 2-seconds of 450-pps envelope rate, down-and up-PSHCs in the absence of the background TEN. Each line represents the evolution of the RMS value of a certain auditory filter output (whose characteristic frequency is indicated on the y axis) computed across successive periods of the envelope. There is a clear amplitude modulation at the output of all filters at a frequency equal to the fundamental frequency (2 Hz). The starting phase of this modulation is, however, different across filters so that at any given time, some filters have a high amplitude and others have a low amplitude. This is illustrated in the bottom panels of Subjects may either selectively attend a particular peak in the excitation pattern and follow its course across the stimulus duration or they may perceive a continuously changing residue pitch based on the continuous analysis of the pattern of peaks present in the excitation pattern (e.g., using a pattern-matching mechanism based on central templates; c.f. [START_REF] Terhardt | Pitch, consonance, and harmony[END_REF]. When the envelope rate is decreased, this has the effect of increasing the ripple density and, therefore, reducing the size of the peaks in the excitation pattern. This is similar to what happens when decreasing the f0 of a harmonic complex filtered in a fixed frequency region: the peaks in the excitation pattern get smaller and smaller as the harmonics become less and less resolved. This may, therefore, explain why the identification task was impossible at 98 pps. The second potential mechanism is temporal. The top panels of Figure 4 show the output of an auditory filter centred on 2254 Hz over 2 seconds (4 periods of the f0) for both PSHCs. The similarity between the outputs obtained for the two PSHCs suggest that there is no information in the envelope at the output of a given auditory filter that can allow the identification of the pitch contour. The middle panels of Figure 4 show a few envelope periods of these same signals. It can be observed that the TFS changes from one envelope period to the next. To quantify this change and to determine the time intervals conveyed by the TFS peaks, the auditory filter output waveform of the first envelope period of the stimulus was cross-correlated with the output waveform of the second envelope period. The lags corresponding to the first three maxima of the cross-correlation function were stored. Then, the second period was cross-correlated with the third period and so on until the end of the stimulus. The bottom panels of Figure 4 show the reciprocal of the lags corresponding to the three highest peaks of these cross-correlation functions as a function of time. These values can be viewed as the possible "pitches" conveyed by the TFS peaks at the output of this specific filter at any given time. This analysis shows that the TFS information within a single auditory filter provides a cue that may allow subjects to perceive the pitch contour of down-and up-PSHCs because the reciprocal of the distance between the TFS peaks progressively decreases for down-PSHC and progressively increases for up-PSHC.

However, these TFS variations are still present at low envelope rates (not shown here), including at 98 pps where subjects could not do the task. Therefore, to explain why performance is getting worse at low rates, we also need to assume that there is a frequency-dependent temporal limitation on the time intervals that can be analyzed by the central auditory system to extract pitch, as already proposed in previous publications [START_REF] Moore | An Introduction to the Psychology of Hearing[END_REF][START_REF] Pressnitzer | The lower limit of melodic pitch[END_REF]. Experiment 2: Measuring the lower limit of TFS sensitivity

Methods

Subjects

Ten subjects with ages ranging from 21 to 41 (mean of 25.5) took part in Experiment 2, including one subject (the author) who had participated in Experiment 1. All subjects were musicians and/or Master's students in Acoustics. They were selected after verifying they had normal absolute hearing thresholds up to 14 kHz. Pure tone detection thresholds at frequencies of 0.25, 0.5, 1, 2, 4, 8, 11.2 and 14 kHz were measured using the same adaptive task as for Expt. 1.

Stimuli

The stimuli were down-and up-PSHCs of various envelope rates filtered into seven different frequency regions. The filtering was identical to that for Expt. 1 except that Fc was equal to 0.25, 0.5, 1, 2, 4, 8 or 11.2 kHz depending on condition. The upper cut-off frequency was set so that the passband width would always be equal to 2 ERBN. The overall presentation level of the complexes was again 40 dB higher than the pure tone detection threshold for a frequency corresponding to the centre of the passband and the level of the TEN noise was 25 dB/ERBN above the detection threshold. Because we did not measure the thresholds in the middle of each passband, this reference threshold was obtained by interpolation from the threshold measurements made for the different values of Fc.

Procedure

The main task was a 3I2AFC, 2-down, 1-up adaptive task. In the first interval, the subjects always heard the down-PSHC. In intervals 2 and 3, the up-PSHC and down-PSHC were presented in random order. Subjects had to indicate which stimulus was different than the other two. The three stimuli presented in each trial had the same envelope rate, so subjects could only use differences in TFS to do the task. The starting envelope rate was equal to 0.2Fc. Following two correct responses, the envelope rate was divided by (1.05) α while following one incorrect response, it was multiplied by (1.05) α . α was equal to 4 until the second reversal, then it was switched to 2 until the fourth reversal and finally to 1 until the tenth reversal at which point the procedure stopped.

While we aimed to use the same f0 of 2 Hz for all stimuli, this was not always possible because, for a given f0, the envelope rate can only be equal to k²f0 where k is an integer and can therefore only take a finite number of values. Therefore, the f0 value closest to 2 Hz was selected for each desired envelope rate. This means that the f0 could slightly vary from trial to trial (from 1.7 to 2.4 Hz in the lowest frequency region and substantially less in higher frequency regions). Note, however, that the f0 was constant within a trial and could not be used as a cue to perform the task.

A measurement block involved collection of seven discrimination threshold measurements, one for each frequency region, presented in random order. Subjects performed between 5 and 6 blocks in total and the first block was considered as training and was not taken into account in the analyses. This experiment was completed in two sessions of 2 hours.

Results

Figure 5A shows the discrimination thresholds (i.e. reflecting the lower limit of TFS sensitivity) obtained as a function of Fc. The lower limit of TFS significantly increased with increases in Fc, ranging from 42 pps for Fc=0.25 kHz to 691 pps for Fc=11.2 kHz. Except for one subject whose threshold was lower at 0.5 than at 0.25 kHz, all subjects showed a monotonic increase as a function of Fc. The other functions illustrate transition rates from good to poor frequency discrimination for harmonic complex tones, replotted from [START_REF] Krumbholz | The lower limit of pitch as determined by rate discrimination[END_REF].

The solid blue line connecting the filled squares in Figure 5B shows the same data expressed as the ratio between Fc and the lower limit of TFS. We refer to this dimensionless variable as the rank of the lowest virtual harmonic present in the passband "NLOW" at threshold as it represents the rank of the lowest-frequency ripple present in the passband. A one-way ANOVA showed a significant effect of Fc on NLOW (F(6,54)=43, p<0.001). Paired t-tests (without correction for multiple comparisons) performed on the data obtained for successive values of Fc showed that NLOW increased significantly from 0.25 to 0.5 kHz (p=0.001), from 0.5 kHz to 1 kHz (p=0.003) and from 1 kHz to 2 kHz (p<0.001), remained approximately constant from 2 to 8 kHz and finally decreased from 8 to 11.2 kHz (p=0.042). Interestingly, the across-subject mean NLOW value was higher than 18 for Fc=2, 4 and 8 kHz and was still higher than 16 for Fc=11.2 kHz. These large NLOW values suggest that the ripples were spectrally unresolved in the passband according to usual definitions of resolvability [START_REF] Moore | Resolvability of components in complex tones and implications for theories of pitch perception[END_REF].

Discussion

As previously explained, our subjects may have used either excitation pattern or phase locking cues to perform the tasks in the two experiments. If subjects based their judgments on excitation pattern cues, we would expect the lower limit of the use of TFS to depend on the bandwidths of the auditory filters. The dashed line of Figure 5B illustrates the limit of harmonic resolvability as defined by [START_REF] Moore | Audibility of partials in inharmonic complex tones[END_REF]; components separated by less than 1.25ERBN are assumed to be unresolved. This limit falls below the lower limit of TFS sensitivity for all regions except the lowest, suggesting that, at threshold, the ripples of our stimuli were unresolved. The dash-dotted line also shows the same "resolvability" limit, but using ERB values proposed by [START_REF] Shera | Revised estimates of human cochlear tuning from otoacoustic and behavioral measurements[END_REF], which are smaller at high frequencies. These ERB values likely represent a conservative definition of resolvability given that they were obtained at lower levels than used here [START_REF] Oxenham | Level dependence of auditory filters in nonsimultaneous masking as a function of frequency[END_REF]. For Fc<4 kHz, the lower limit of TFS still falls above this line, confirming that the ripples at threshold were unresolved even assuming narrow auditory-filter bandwidths. For Fc>4 kHz, it cannot be excluded that the ripples were partially resolved. However, if the subjects based their judgments on excitation pattern cues, we would expect NLOW for the lower limit of TFS sensitivity to continue to increase up to the highest Fc. This was not the case as NLOW remained mostly constant above 2 kHz. Also, it should be noted that the random fluctuations produced by the background TEN should reduce the ability to use excitation-pattern cues [START_REF] Jackson | The role of excitation-pattern, temporal-fine-structure, and envelope cues in the discrimination of complex tones[END_REF]. Thus, unless our subjects managed to use residual excitation pattern cues in the lower skirt of the filter, these data suggest that the lower limit of the use of TFS is related to the use of temporal (phase-locking) cues rather than excitation-pattern cues. Moore et al. (2009) measured sensitivity to TFS for several f0 values in different frequency regions. They showed that for stimuli bandpass filtered so that the lowest component in the passband, NLOW, was harmonic 13, subjects could discriminate harmonic and inharmonic tones for f0=50, 100 and 200 Hz. At these f0, harmonic 13 corresponds to Fc=1.3, 2.6 and 5.2 kHz, respectively. Our results are consistent with their data because at these Fc values, our subjects could always discriminate between down-and up-PSHCs when NLOW was the 14 th or lower (see Figure 5B). They also showed that subjects required the presence of lower harmonics to perform the same task for f0=35Hz. Similarly, Figure 5B shows that our NLOW was lower in low frequency regions. [START_REF] Tarnowska | The effect of musicianship, contralateral noise, and ear of presentation on the detection of changes in temporal fine structure[END_REF] showed above-chance performance for musicians on a variant of the TFS task used in Moore et al. (2009) for Fc=3.8 kHz and NLOW=19 which is also consistent with the present data. Finally, [START_REF] Moore | Sensitivity of the human auditory system to temporal fine structure at high frequencies[END_REF] extended these measurements to higher fundamental frequencies (f0=800 and 1000 Hz) and showed sensitivity to TFS when NLOW was equal to 12 although only 6 out of their 8 subjects could perform the task at the highest f0. Their stimuli are equivalent to Fc=9.6 and 12 kHz, respectively, and can be compared to our highestfrequency region (Fc=11.2 kHz) for which subjects could discriminate down-and up-PSHCs for NLOW =16 and lower. It is possible that the poorer performance [START_REF] Moore | Sensitivity of the human auditory system to temporal fine structure at high frequencies[END_REF] observed for Fc=12 kHz than for Fc=9.6 kHz is related to the degradation of phase locking cues at very high frequencies. The trend of decreased performance we observed in the highest frequency region may have the same origin.

In summary, the finding that subjects were still sensitive to TFS when NLOW was equal to 18 for Fc=2, 4 and 8 kHz suggests that the ripples were unresolved and that subjects used temporal cues to perform the task. These results, therefore, confirm the findings of [START_REF] Moore | Sensitivity of the human auditory system to temporal fine structure at high frequencies[END_REF] and add further evidence that phase locking cues can be used by humans at very high frequencies.

Finally, it has been shown that the frequency discrimination of harmonic complex tones depends on the rank of the lowest harmonic present in the passband. When the frequency region is fixed, f0 difference limens (DLs) are poorer at low than at high f0 [START_REF] Ritsma | Frequency selectivity and the tonal residue[END_REF][START_REF] Cullen | Rate discrimination of high-pass-filtered pulse trains[END_REF][START_REF] Krumbholz | The lower limit of pitch as determined by rate discrimination[END_REF]. The transition from good to poor performance, however, depends on frequency region. [START_REF] Krumbholz | The lower limit of pitch as determined by rate discrimination[END_REF] compiled the results of several of these studies (c.f. their Figures 7 and10) and defined the transition rate as the rate at which the f0 DL equals 2.5% except for the [START_REF] Ritsma | Frequency selectivity and the tonal residue[END_REF] data for which the criterion was set to 1%. These transition rates, expressed in terms of the rank of the lowest harmonic number in the passband are replotted in Figure 5B. The upward triangles also show the lower limit of melodic pitch defined as the rate at which subjects could no longer recognize that the f0 of a complex tone in a melody changed by a semitone [START_REF] Pressnitzer | The lower limit of melodic pitch[END_REF]. The functions showing these transition rates as a function of Fc follow a similar pattern as the lower limit of TFS, first showing an increase followed by a decrease or a flat portion above 1 kHz. Although speculative, it is possible that the transition from good to poor f0 discrimination of complex tones with increasing harmonic rank reflects the decreased sensitivity to TFS in high-frequency regions.

Conclusions

These experiments show that it is possible to generate pairs of harmonic complexes that share the same long-term spectrum and the same temporal envelope but have different TFS and elicit different pitch percepts. These down-and up-PSHCs add to the arsenal of tools for studying sensitivity to TFS. The lower limit of TFS sensitivity (expressed as the envelope rate below which subjects were insensitive to TFS) was shown to increase with increasing frequency region. The observation that subjects were still sensitive to TFS in high frequency regions (> 4 kHz) when the stimuli likely contained only unresolved components suggests that phase locking cues remain usable at high frequencies, confirming the results of previous studies.

Figure 1 -

 1 Figure 1 -Top: In blue: temporal waveforms of three stimuli used in Experiment 1, all at a rate of 450 pps. The red line shows the Hilbert envelope. Only 3 envelope periods are shown. Bottom: Spectrograms showing 2 seconds of the same three signals.

Figure 2 -

 2 Figure 2 -Top: Results of Expt. 1a showing the percentage of correct identification at each of the three envelope rates for each subject. Bottom: Results of Expt. 1b showing the percentage of correct discrimination of the three pairs of stimuli for each subject.

Figure 3

 3 Figure3which show excitation patterns calculated across a single envelope period at three successive times across the stimulus. These excitation patterns represent the RMS outputs of the Gammatone filters expressed on a dB scale. All excitation patterns show peaks separated by a frequency equal to the envelope rate. However, the peaks drift during the stimulus period. For the down-PSHC, the peaks shift progressively downwards across the stimulus duration while for the up-PSHC, they shift upwards. The perception of down-PSHCs as falling in pitch and of up-PSHCs as rising in pitch may reflect the use of excitation pattern cues in two different ways:

Figure 3 -

 3 Figure 3 -Top: RMS energy at the output of several Gammatone filters in response to downand up-PSHCs with an envelope rate of 450 pps. The RMS energy was calculated across successive envelope periods. The red function corresponds to the filter located in the center of the passband (2254 Hz). Bottom: Excitation patterns at three successive times (t1, t2 and t3) during a period of the stimulus for the same signals.

Figure 4 -

 4 Figure 4 -Top panels: Output of a gammatone filter centered on 2254 Hz (center of the passband) for the two PSHCs. Middle panels: Short portions of the same output showing the TFS. Bottom panels, Reciprocal of the lags corresponding to the three highest peaks of the cross-correlation function between one envelope period and the next. Blue, red and lack symbols are for the highest, second-highest and third-highest peaks, respectively.

Figure 5 -

 5 Figure 5 -A: Boxplots showing the distributions of the lower limit of TFS sensitivity for each Fc. Each boxplot shows the median, 25 th and 75 th percentiles and extreme values. Outliers (red '+' symbols) correspond to values that are more than 1.5 times the interquartile range away from the bottom or top of the box. The solid line shows the mean across subjects. B, The same data as in panel A are plotted in blue as a function of the lowest virtual harmonic in the passband (see main text). The dashed and dashed-dotted lines show the limit of harmonic resolvability for two sets of auditory filters (see main text for the definition of resolvability).

Acknowledgements

I thank Bob Carlyon and Brian Moore for their comments on a previous version of this manuscript.


Juan D'adamo

Manuel Collaud

Roberto Sosa

Ramiro Godoy-Diana

Wake and aeroelasticity of a flexible pitching foil

A flexible foil undergoing pitching oscillations is studied experimentally in a wind tunnel with different imposed free stream velocities. The chord-based Reynolds number is in the range 1600-4000, such that the dynamics of the system is governed by inertial forces and the wake behind the foil exhibits the reverse Bénard-von Kármán vortex street characteristic of flapping-based propulsion. Particle Image Velocimetry (PIV) measurements are performed to examine the flow around the foil, whilst the deformation of the foil is also tracked. The first natural frequency of vibration of the foil is within the range of flapping frequencies explored, determining a strongly-coupled dynamics between the elastic foil deformation and the vortex shedding. Cluster-based reduced order modelling is applied on the PIV data in order to identify the coherent flow structures. Analysing the foil kinematics and using a control-volume calculation of the average drag forces from the corresponding velocity fields, we determine the optimal flapping configurations for thrust generation. We show that propulsive force peaks occur at dimensionless frequencies shifted with respect to the elastic resonances that are marked by maximum trailing edge oscillation amplitudes. The thrust peaks are better explained by a wake resonance, which we examine using the tools of classic hydrodynamic stability on the mean propulsive jet profiles.

I. INTRODUCTION

An oscillating elastic plate constitutes probably the simplest model of a bio-inspired propulsive mechanism: the back and forth motion of a structure subject to deformation as it interacts with the surrounding fluid is the basic fluid-dynamical problem behind animal swimming and flying [see e.g. 1, 2, and references therein], regardless of the broad collection of different geometrical and material properties, as well as the specificities of different kinematics with different degrees of complexity involved in insect [START_REF] Wootton | Functional morphology of insect wings[END_REF][START_REF] Bomphrey | Insect and insect-inspired aerodynamics: unsteadiness, structural mechanics and flight control[END_REF], bird [START_REF] Biewener | In vivo strain in the humerus of pigeons (Columba livia) during flight[END_REF][START_REF] Tobalske | Biomechanics of bird flight[END_REF] and bat [START_REF] Song | Aeromechanics of membrane wings with implications for animal flight[END_REF] wings [START_REF] Chin | Flapping wing aerodynamics: from insects to vertebrates[END_REF], and fish fins and tails [START_REF] Mchenry | Mechanical control of swimming speed: stiffness and axial wave form in undulating fish models[END_REF][START_REF] Blake | Fish functional design and swimming performance[END_REF].

From an engineering perspective, flapping-based propulsion has come back to the center stage in the past two decades especially with the development of micro-air-vehicles (MAVs)-fuelled by miniaturisation-see e.g. [START_REF] Whitney | Conceptual design of flapping-wing micro air vehicles[END_REF][START_REF] De Croon | The delfly[END_REF][START_REF] Chen | A Review of Design and Fabrication of the Bionic Flapping Wing Micro Air Vehicles[END_REF] for a review. With respect to conventional aerodynamics, the main difference is rooted in the unsteady mechanisms of force production by a flapping wing. Average thrust and lift are the outcome of the periodic flapping motion of a structure that accelerates and decelerates, as opposed to the case of fixed-wing aircraft, where thrust is produced by a jet engine or propeller and lift is the result of the flow around the static wing. Moreover, these unsteady mechanisms are intimately linked to the problem of vortex shedding, especially the forced vortex shedding that occurs on the time scale clocked by the flapping characteristic frequency. The distinctive feature of the wake of a flapping foil is a sequence of counter-rotating vortices that bears resemblance to the vortex street behind a bluff body but where, in the regimes of interest for propulsion, the sign of vorticity is inverted: the reverse Bénard-von Kármán (BvK) vortex street [START_REF] Freymuth | Propulsive vortical signature of plunging and pitching airfoils[END_REF][START_REF] Anderson | Oscillating foils of high propulsive efficiency[END_REF][START_REF] Godoy-Diana | Transitions in the wake of a flapping foil[END_REF]. The average flow field behind a flapping foil producing thrust is a jet flow. In the case of a body in self-propulsion, as the body accelerates towards a cruising regime, the jet profile behind the flapping body tends to a momentum-less wake profile (typically an inner jet surrounded by a drag wake) where the average thrust is balanced by the global drag [START_REF] Afanasyev | Wakes behind towed and self-propelled bodies: asymptotic theory[END_REF][START_REF] Arbie | Stability of momentumless wakes[END_REF]. In the case of flapping flight, the aerodynamic force produced by the flapping wings is directed obliquely downwards, and it is partly used for thrust and partly as lift to counter the weight of the flyer, the limit case being that of hovering, where all the thrust is directed downwards and there is no cruising speed.

In addition to the imposed flapping motion, the other main ingredient governing the dynamics of wings, fins and bio-inspired appendages is their structural deformation response [see e.g. [START_REF] Marais | Stabilizing effect of flexibility in the wake of a flapping foil[END_REF][START_REF] Huera-Huarte | On the impulse produced by chord-wise flexible pitching foils in a quiescent fluid[END_REF]. The latter is determined by the material and geometric properties of the body and it is the basic solid mechanics problem of forced vibrations. Without the surrounding fluid, i.e. for a body flapping in vacuum, the whole problem is defined by the interplay between the solid inertial forces that drive the deformation and the elastic restoring force. Adding the fluid to this forced vibration problem means, from the point of view of the structure, new fluid forces (viscous and inertial). But from the point of view of the fluid there is a whole new set of equations of motion to be considered, where the elastic solid determines the boundary conditions. Physically, the flapping foil is thus a fluid-structure interaction problem governed by a few dimensionless parameters that measure the respective importance of inertial, viscous and elastic forces in both the solid and the fluid. We can be more specific considering a flat plate of chord length L, span H, mass per unit surface µ, bending stiffness B = EI as a model, immersed in a fluid of density ρ and kinematic viscosity ν, and actuated in pitching oscillations of amplitude θ R and frequency f f . The main dimensionless parameters are then: the mass ratio M = ρ ρ s , which compares the density of the fluid to that of the solid; the Reynolds number Re = U L ν, which compares fluid inertial vs. viscous forces-where U is a characteristic velocity that can be the average flying or swimming speed in a cruising regime, or the flapping characteristic velocity defined in terms of f f and θ R ); and the dimensionless frequency f + = f f f n , that compares the forcing frequency f f to the natural frequency of the first mode of vibration of the structure f n and thus contains the information about the bending stiffness.

The parameter space defined by M, Re and f + encompasses a wide variety of problems, ranging from relatively stiff flapping wings of some insects [START_REF] Combes | Flexural stiffness in insect wings i. scaling and the influence of wing venation[END_REF] to the very flexible parapodia of swimming snails [START_REF] Mohaghar | Characteristics of swimming shelled antarctic pteropods (limacina helicina antarctica) at intermediate reynolds number regime[END_REF]. In this paper we focus on an aerodynamic problem with M ∼ 10 -3 and Re ∼ 10 3 , and we examine the case of a rectangular foil of span-to-chord aspect ratio of ∼ 4.6 bending chord-wise mostly on its first mode of deformation as in [START_REF] Thiria | How wing compliance drives the efficiency of self-propelled flapping flyers[END_REF][START_REF] Ramananarivo | Rather than resonance, flapping wing flyers may play on aerodynamics to improve performance[END_REF]. We use two-dimensional (2D) Particle Image Velocimetry (PIV) on the symmetry plane at mid-span to examine the flow in the wake and cluster-based reduced order modelling to identify the coherent flow structures. Simultaneously, the elastic deformation of the foil is tracked. The first natural frequency of vibration of the foil is within the range of frequencies of the forcing pitching oscillations, which determines a strong coupling between the structural dynamics of the elastic foil and the vortex shedding in the wake. We use a control-volume calculation based on the PIV data in order to compute the average drag force for each flapping configuration, which allows us to determine optimal points in the parameter space in terms of thrust generation. We show that propulsive force peaks occur at dimensionless frequencies shifted with respect to the elastic resonances, which are marked by maximum trailing edge oscillation amplitudes. The thrust peaks are better explained by a wake resonance [START_REF] Triantafyllou | Optimal thrust development in oscillating foils with application to fish propulsion[END_REF][START_REF] Moored | Hydrodynamic wake resonance as an underlying principle of efficient unsteady propulsion[END_REF][START_REF] Moored | Linear instability mechanisms leading to optimally efficient locomotion with flexible propulsors[END_REF], which we examine using the tools of classic hydrodynamic stability on the mean propulsive jet profiles [START_REF] Triantafyllou | On the formation of vortex streets behind stationary cylinders[END_REF].

II. CONFIGURATION AND EXPERIMENTAL SETUP

The flow configuration is described in Fig. 1a for a closed loop wind tunnel that has a test section of dimensions L x = L y = L z = 180 mm. A flexible foil of chord length L = 35 mm and span s = 160 mm undergoes an oscillatory rotation, animated by a DC motor through a crank mechanism. The parameters of the movement are a forcing angular frequency ω f = 2πf f and a fixed angular amplitude θ R = 10 ○ . We choose an aspect ratio s ∶ L ≃ 4.6 that favors 2D flow structures even though the wake is 3D. Concerning flow blockage, for rigid motions, the frontal length of the foil L y ∶ 2L sin(θ) ≃ 14.8 assures that the wind tunnel walls should not affect the wake flow dynamics.

As described by Figure 1b, the total displacement of the foil w(x, t) can be written as the sum of a rigid motion xθ and an elastic deformation w, i.e. w = w + xθ, as illustrated on Fig. 1b. Considering forced oscillations included as θ = θ 0 sin(ω f t), it is possible to characterize the dynamics using an Euler-Bernoulli beam framework for the displacements w with times and spatial derivatives correspondingly ẇ, w ′ . We can write the conservation of momentum for the elastic foil as:

µ ẅ + Bw iv + C Dp L ẇ + w ′ U ( ẇ + w ′ U ) = 0 ( 1
)
where µ is the mass per unit length, B = EI is the flexural rigidity composed of the Young modulus E and the moment of inertia I. C Dp is the drag coefficient for the flow normal to a plate, a resistive contribution as proposed by [START_REF] Taylor | Analysis of the swimming of long and narrow animals[END_REF] and later revisited by [START_REF] Ramananarivo | Rather than resonance, flapping wing flyers may play on aerodynamics to improve performance[END_REF][START_REF] Eloy | The origin of hysteresis in the flag instability[END_REF][START_REF] Piñeirua | Modelling of an actuated elastic swimmer[END_REF]. In our case, we have contributions from the flapping motion ẇ and from the projection of the free stream flow w ′ U on the same direction. Reynolds numbers based on the chord length are in the range Re ∼ [1600, 4000], corresponding to the working velocities of the wind tunnel of U ∈ [0.7, 1.8]m/s. Because viscous friction and added mass have small contributions, fluid forces on the foil are in practice given by the last term of Eq. 1. The expression can be further simplified if we consider that the first mode of elastic deformation is dominant. It is thus sufficient to study the projection of the trailing edge displacements δ = w(L, t).

We have chosen the physical parameters of the foil and the kinematic parameter space in order to define a problem with strong coupling between the fluid and solid elastic dynamics. On the one hand, considering a frontal area of the flapping object given by a characteristic length times the span s, a vortex shedding frequency can be roughly estimated using the Strouhal number

St = f vs U = f vs 2δ R U (2)
and a reference value for bluff bodies St ≈ 0.2. For the range of wind tunnel speeds of the present experiment, the vortex shedding frequencies-which are here driven by the frequency of the forcing pitching oscillation (i.e. f vs ≡ f f ) are thus in the range [10 . . . 30]Hz. On the other hand, considering the foil as a cantilever beam, its elastic natural frequency is f n = k 2 B µ 2π, where the wavenumber k is such that kL = 1.875 for the first elastic bending mode. For a 50µm polyethylene foil, the elastic natural frequency is f n ≃ 12.27 Hz. This value was confirmed experimentally through image processing of free damped vibrations of the foil, which gave f n = 11.73 Hz.

For image acquisition, we used a monochromatic high-speed video camera (Weinberger SpeedCam in a mode of 400 frames/s with a spatial resolution of 512×512 pixels 2 ). We also performed Particle Image Velocimetry measurements using a LaVision system, composed of an ImagerPro 1600 × 1200 CCD camera with 14-bit dynamic range capable of recording double-frame pairs of images at 14 Hz and a two-rod Nd-YAG (15mJ) pulsed laser. The field of view was 170mm×60mm (about 5L× 2L), with a spatial resolution of 0.8mm (0.023L) and the number of snapshots for each case was 680, in order to assure statistical convergence. A typical instantaneous vorticity field ω(x, t) is presented in Fig. 1c, where its contour levels follow the foil deformation. Two rows of eddies exhibit the characteristic arrangement of the reverse BvK vortex street, with positive (negative) eddies placed in regions where y > 0 (y < 0). This configuration produces a net momentum injection in the wake, typical of propulsion produced by flapping motion, with a jet-like structure noticeable by averaging in time the corresponding component of the velocity field, u x (x).

Assuming ergodicity, time-averaged velocity fields ⟨u(x)⟩ and fluctuations intensity of the j velocity component ⟨u ′ j 2 (x)⟩ were determined respectively as:

⟨u(x)⟩ = 1 N N i=1 u(x, t i ) ⟨u ′ j 2 (x)⟩ = 1 N N i=1 (u j (x, t i ) -⟨u j ⟩) 2 (3)
Then, with this definition, mean net momentum induced by flapping, u = u x -U is illustrated as contour levels on Fig. 1d.

The PIV measurements in the present experiment are not time resolved, since the maximum sampling frequency of the system f S ∼ f nat ∼ f vs . However, the flow temporal dynamics was accurately recovered using a clustering method. The algorithm is based on the k-means algorithm (see Appendix A 2 for details) and was used to classify the acquired PIV fields and obtain a number of different representative states. As clustering is similar to phase averaging, the velocity and vorticity fields are more clear than direct snapshots, and the whole period of flapping can be reconstructed.

III. RESULTS

A. Flexible foil dynamics

We proceed to characterize the response of the elastic foil by varying the forcing frequency f f or its non-

dimensionalized form f + = f f f n in a range f + ∈ [0.5, 1.8]. We define the reduced velocity U R = U (2πf n L)
that expresses the coupling between the free flow U and the elastic foil dynamics.

The results are presented in Fig. 2a, where the non-dimensional displacements of the trailing edge of the foil δ + are obtained by comparing the total to the rigid displacements δ + = δ δ R . Similar curves are obtained for each reduced velocity U R , and the data points are fitted with a third order polynomial around the maximum of displacement δ + max . The corresponding maxima (f + max , δ + max ) are plotted on Fig. 2b their relationship, together with a second-order polynomial fit. In the inset of the same Figure, a linear dependence between the reduced velocity U R and the maximum of amplitude δ + max is displayed. It can be seen that the maximum deformation δ + max decreases and f + max shifts with U R . Modeling the foil dynamics as a nonlinear oscillator with damping contributions from ẇ and U R , it is reasonable to find such behavior. This could be done by reducing eq.( 1) to a system of equations regarding trailing edge dynamics. Through these results it is possible the approximation of displacements for different flow conditions by interpolation, criterion supported by Fig. 2c, where a scaling can be performed showing that the regime is similar for each reduced velocity, Another way to scrutinize the kinematics is to compare the changes on the Strouhal number defined in eq.2 with respect to a different Strouhal number defined using the elastic displacement 2δ, St * = 2δf f U . Fig. 2d shows a family of similar curves for the different reduced velocities. More insight on the fluid-structure coupling can be gained recalling the different terms in eq. ( 1) for the Euler-Bernoulli model of the elastic displacements of the foil. We define a Cauchy number

C * Y = ( ẇ + U w ′ ) 2 (ω n L) 2 (4)
that compares the kinetic energy of the fluid velocities (we can identify the flapping velocity ẇ that is responsible for a major component of damping drag, and also the contribution from the free stream velocity coupled by the angle, U w ′ ) to that of the solid elastic vibrations given by (ω n L) 2 , for each reduced velocity U R . Both ẇ and w ′ are evaluated at the trailing edge of the foil, therefore, C * Y (U, f f , w ′ ) is a function that changes for every measured case. The inset of Fig. 2d shows that the Cauchy number C * Y provides a satisfactory scaling for the Strouhal numbers relationship.

B. Flow field results

The flow dynamics in the wake of the foil was inspected for different parameters (f + , U R ) through PIV measurements. A representative state of the flow was already presented in Fig. 1c (c) for the case of U R = 0.29 and f + = 1.47, where the well-known reverse BvK vortex street regime is established. A first characterization of the induced velocity field can be obtained by representing contour lines of the difference between time mean averaged and free flow as shown in Fig. 1d. A jet-like flow structure is the result in each forced case, and we can compare them conveniently by plotting a profile at y = 0 as displayed on Fig. 3a. For U R = 0.29 an increment in the overall magnitude of the velocity profiles is observed up to f + = 1.47. This is expected as ω f = 2πf f L gets higher, but effective flapping amplitudes reduce for f + > 1.3 recalling Fig. 2a. Maxima for the curves are summarized in Fig. 3b, where it is clear that for increasing U R , induced motion diminishes monotonically. On the other hand, flapping also produces transversal velocity fluctuations ⟨u ′ y 2 ⟩ 1 2 which are displayed in Fig. 3c. These augment with f + and attain a maximum for a value in the range 1.30

< f + < 1.63].
Furthermore, Fig. 3d outlines the evolution of transversal velocity fluctuations maxima for each forced case. We shall analyze how these components constrain effective streamwise momentum production or, in a more general context, how they are related to the determination of drag forces. We will also later discuss this aspect in terms of the wake receptivity or its hydrodynamical stability properties.

C. Drag forces

Time-averaged aerodynamic forces, particularly drag/thrust in the streamwise direction can be estimated from PIV 2D-velocity field data. There are different approaches to achieve this goal in the literature (see for instance [START_REF] Kurtulus | Unsteady aerodynamic forces estimation on a square cylinder by tr-piv[END_REF]). Here we start by considering the general framework depicted in Fig. ,4a that consists in integrating both the momentum budget in a control volume V enclosing the flapping foil and the pressure forces on the corresponding boundary surface S . The expression for the force reads:

F(t) = - D Dt V ρu dv + S (-pI + T) ⋅ n ds, (5)
(y = 0)⟩ 1 2 for each (f + , U R ) case.
FIG. 3: Forced wake flow characteristics.

where p is the pressure field, I the unit tensor, T = µ(∇u + ∇ T u) is the viscous stress tensor (with viscosity µ) and n is a unit vector orthogonal to the boundary S .

Neglecting the effect of inhomogeneities in the flow along the spanwise direction, the projection of the time-averaged Eq. ( 5) in the streamwise direction yields

F D = +∞ ∞ u x (1 -u x )dy + 1 ρ +∞ ∞ ∆pdy - +∞ ∞ ⟨u ′ x 2 ⟩dy (6)
where, for the sake of simpler notation F D ≡ ⟨F D ⟩, u x ≡ ⟨u x ⟩ and in the equation the integral is evaluated on a fixed wake coordinate x = x wake . Viscous terms are therefore very small compared with convective and pressure

-2 0 2 4 6 -2 0 2 x wake Trefftz u x (x, t) U S x L (a) -1 0 1 -1 -0.5 0 0.5 y L ⟨u x ⟩ (1 -⟨u x ⟩) ⟨u ′ y ⟩ 2 -⟨u ′ x ⟩ 2 (b)
FIG. 4: (a) Schematic representation of a control volume approach to determine forces from velocity fields. (b) Momentum fluxes for drag estimation in equation ( 6) for U R = 0.29, f + = 1.47 at x wake = 2.5L. ∆p = (p(x wake , y) -p ∞ ) terms. The main difficulty of this expression is the estimation of the pressure field from the 2D velocity field. This is achieved either by means of the Poisson equation (see e.g. [START_REF] Fujisawa | Evaluation of pressure field and fluid forces on a circular cylinder with and without rotational oscillation using velocity data from piv measurement[END_REF]), or by integrating the Navier-Stokes (NS) equation along the control surface [START_REF] Unal | Force prediction by piv imaging: a momentum-based approach[END_REF]. The time-averaged pressure gradient on x = x wake is given by:

1 ρ ∂⟨p⟩ ∂y = ⟨u x ⟩ ∂⟨u y ⟩ ∂x + ⟨u y ⟩ ∂⟨u y ⟩ ∂y + ∂⟨u ′ x u ′ y ⟩ ∂x + ∂⟨u ′ y u ′ y ⟩ ∂y (7)
where, for wakes, the last term is much larger than the rest, so ∆p ≈ ρ⟨u ′ y 2 ⟩ accounts for the change of pressure in the wake due to transverse velocity fluctuations gradients. Figure 4b presents an example of how the terms in Eq. ( 6) can be evaluated at a position x wake and integrated in order to obtain drag forces. It can be seen that the transverse velocity fluctuations (orange circles) bring a contribution of the order of the momentum flux in the stream direction (blue squares).

A possible approach to explicitly compute the force was proposed by [START_REF] Noca | Measuring instantaneous fluid dynamic forces on bodies, using only velocity fields and their derivatives[END_REF] and later revisited by [START_REF] Wang | A lift formula applied to low-reynolds-number unsteady flows[END_REF], where pressure terms can be replaced by products of velocity and vorticity impulse contributions (for a formulation in the case of a cylinder wake, see e.g. [START_REF] Adamo | Spatio-temporal spectral analysis of a forced cylinder wake[END_REF]).

However, in order to obtain accurate results, the evaluation of the flow field at a position x wake far enough from the body is needed, which is out of reach of the PIV measurement window of the present experiments.

To circumvent this problem, we used an alternative method following the study presented in [START_REF] Hall | Minimum induced power requirements for flapping flight[END_REF] (see also [START_REF] Minotti | Determination of the instantaneous forces on flapping wings from a localized fluid velocity field[END_REF]). Through a control volume scheme, under the hypothesis of periodicity, for a wake flow expressed as u = U i + ∇φ = (U + u)i + vj + wk, a superposition of a potential flow ∇φ and a far field velocity U , it is possible to estimate timeaveraged forces on a flapping body. We recall, that in the context of wake flows, the Reynolds numbers of these experiments Re ∼ 10 3 allow us to consider that starting at the very near wake, inertia is dominant with respect to viscous forces.

A light loading hypothesis applies when the unsteady wake and its associated potential are simply convected in the x-direction with speed ∼ U , given that ∇φ ≪ U . To leading order, the potential is φ = φ (x -U t, y, z). Furthermore, the unsteady flow is periodic in the x-direction so averaging over one temporal period is equivalent to averaging over one spatial period. It results that the momentum balance reduces to:

F = - 1 T ∭ V ρ (ui + vj + wk) dV ( 8
)
where V is the Trefftz volume, which encloses an x-periodic control volume in the wake, characterised by a length Trefftz = U T = U f f , as showed in Fig. 4a. In non-dimensional form, 8) cannot be performed.

˜ Trefftz = Trefftz L = U T L = U f f L .
Recalling that St = f f 2δ R U and δ R = L sin(θ R ), it follows that ˜ Trefftz = U f f L δ R δ R = U f f δ R sin(θ R ) = 2 sin(θ R ) St ( 9
)
We can evaluate the integral in Eq. ( 8) on the volume defined by this length and divide by the period T or multiply equivalently by f f or St in non-dimensional form. The result of this computation expressed as a drag coefficient C D = 2F x (ρU 2 2δ R ), considering homogeneity in the spanwise direction, is shown in Fig. 5. The computed forces are plotted as a function of the forcing frequency for all the different values of the reduced velocity U R on Fig. 5a. It is evident that as U R increases, the momentum produced from flapping is smaller compared to the incoming flow, so in U R = 0.29 we appreciate the highest thrust (or the minimum drag). Among these values, a minimum is attained at a forcing frequency around f + = 1.47. We note that this vaue does not match to the solid elastic resonance f + = 1 nor to the maximum elastic displacement observed in Fig. 2a at f + = 1.30.

Trefftz volumes, defined by the associated lengths Trefftz are determined from the experimental parameters and plotted on Fig. 5b. For some cases, U T becomes higher than the PIV domain behind the foil, so we cannot use this scheme. Nonetheless, the relationship between drag and forcing frequency presents strong changes for low relative velocity values, U R = [0.29; 0.37], where Trefftz volumes are well defined. For the other cases, the induced momentum, as already depicted in Fig. 3b, is just high enough to compensate the losses in the wake due to transversal fluctuations.

In this context, previous works have hypothesized on the causes of an optimal propulsion. Triantafyllou et al. [START_REF] Triantafyllou | Optimal thrust development in oscillating foils with application to fish propulsion[END_REF][START_REF] Triantafyllou | Wake mechanics for thrust generation in oscillating foils[END_REF], for the case of wake dynamics of a pitching rigid foil postulated that the frequency of maximum spatial amplification in the wake provides optimal thrust production per input power. In the same sense, later, Moored et al. [START_REF] Moored | Hydrodynamic wake resonance as an underlying principle of efficient unsteady propulsion[END_REF][START_REF] Moored | Linear instability mechanisms leading to optimally efficient locomotion with flexible propulsors[END_REF] experimentally showed that resonant peaks in thrust occurred for discrete values of the effective flexibility, a non-dimensional parameter measuring the ratio of added mass forces to internal bending forces. In their study, not only friction but added-mass contributions play a role on more complex dynamics. Their linear stability analysis used the Orr-Sommerfeld equation as viscous terms cannot be neglected and the flow finds more combinations of resonant frequencies for the wakes. They found on the one hand that for flexible propulsors each peak in efficiency occurs when the driving frequency of motion is tuned to a wake resonant frequency, not a structural resonant frequency. On the other hand, panel flexibilities that attain global optimally efficient locomotion are those for which structural and wake resonant frequencies are tuned.

Although, the analysis is focused on a self-propelled condition, we found that some of these concepts can be applied to our present results, as we discuss in what follows.

U θ ω f L δ x y Γ U ∼ U cos(θ)δ Γ u θ ∼ ω f Lδ (b)
FIG. 6: a) Non-dimensional circulation as a function of forcing frequency. Filled symbols come from PIV measurements and non-filled symbols from the kinematic model represented schematically in b).

The model predicts circulation from the parameters of the foil kinematics and the free flow (see text).

D. Drag forces and Circulation

Considering a potential flow theoretical framework, the forces induced on the body are proportional to each element of the vortex street circulation Γ i and to their distance y i respect the symmetry axis y = 0, as described by the impulse equation :

F = d dt Γ i y i (10)
The circulation Γ i = ∫ s ωds can be determined for each case as we have clear quantitative results of vorticity and vortex definition from the PIV cluster analysis. The results are displayed on Fig 6a . An alternative estimate for the circulation can be obtained using the foil kinematics. Since there is a competition between vorticity production from the free flow and from flapping, as depicted by Fig. 6b, Γ i = ∫ s ū dl, we express:

Γ i = Γ uθ -Γ U = ω f Lδ -U cos(θ)δ (11)
The values computed from this kinematic approximation are also presented in Fig. 6a and a reasonable agreement with the calculation from the PIV fields is found for most cases. Also, the kinematic computation brings an estimate of the circulation where no PIV data is available, thereby extending the range of frequencies explored. For U r = 0.29 the circulation prediction from the kinematic model exceeds the calculation from the PIV measurements for frequencies over the maximum, f + > 1.30. Maxima are delayed for increasing U R so differences in other cases are not so clearly displayed. We hypothesize that the wake is less sensitive to higher frequencies so forcing at f + > f + max ceases to create effective vorticity.

Examining U R = 0.29, for f + = 1.14 and f + = 1.47 we observe that while the value of vortex circulation from PIV measurements does not change significantly, the drag coefficient on Fig. 5a is about the double. Recalling Eq. [START_REF] Blake | Fish functional design and swimming performance[END_REF], this fact confirms that the change in the resulting force must be due to a modification in vortex arrangements.

The spatial development of the vorticity field is depicted in Fig. 7, for a similar forcing phase period under three different forcing frequencies. We observe that the intensity and position of vortices in the near wake does no differ from one case to the other but for larger distances (x L > 2.5) f + = 1.47 shows a more ordered structure. Traces for every phase of the period are also plotted for both cases, describing crossing (f + = 1.14; 1.30) and parallel (f + = 1.47) vortex streets. In an air flow, for Re ∼ 2000, where we can neglect friction and added mass contributions on the wake dynamics, instability of these trajectories can be attributed to an inviscid mechanism. We study this particular aspect of forced wakes by means of a hydrodynamical stability analysis based on previous approaches [START_REF] Triantafyllou | On the formation of vortex streets behind stationary cylinders[END_REF][START_REF] Thiria | Stability properties of forced wakes[END_REF] and hypotheses [START_REF] Moored | Hydrodynamic wake resonance as an underlying principle of efficient unsteady propulsion[END_REF][START_REF] Moored | Linear instability mechanisms leading to optimally efficient locomotion with flexible propulsors[END_REF].

E. Hydrodynamical stability

In order to analyze the stability properties of the wake, we explore its receptivity to disturbances in the region x L > 1. Neglecting viscosity effects, time mean average velocity profiles u x (y) are analyzed with the inviscid Rayleigh equation. Global frequencies for the wake are obtained for two forcing cases U R = 0.29 and 0.37 where clear propulsive force peaks were observed. Details on the method and implementation are included in Appendix A 1.

Figure 8 summarizes the results of the analysis. Both panels in the figure share the same abscissa f + . The top plot reproduces the drag coefficient curves from Fig. 5a whilst the bottom panel presents the non-dimensional frequencies selected by the wake (written as Strouhal numbers), as obtained by the stability analysis. The horizontal bars represent Strouhal numbers calculated as f n 2δ R U while vertical bars mark values of drag minima for each reduced velocity. These lines coincide in the vicinity of non-dimensional frequencies issued from linear stability analysis.

It is thus resonance (e.g. U R = 0.29, f + = 1.47) between the stability properties from the wake and the elastic vibrations from the foil the criteria that leads to maximum propulsion (or drag minima). Such condition is in agreement with parallel trajectories of vortices as already depicted on Fig. 7c. Again, y i in eq. (10 remains constant and assures the maximum impulse. On the contrary, when foil elastic and wake frequencies are not synchronized (i.e. (U R = 0.29, f + = 1.14)), even for similar values of circulation produced in the near wake as measured and predicted on Fig. 6a, vortex paths result unstable producing higher transversal fluctuations and less added net momentum. From these facts, we found that the interaction between the foil induced motion and the free flow (through U R ) is not only in terms of non-linear damping but also in setting the wake stability conditions. The shaded horizontal lines correspond to the elastic non-dimensional frequencies of the foil for each case. Resonance between the wake and the elastic foil occurs when they intersect the Strouhal curves. The top plot reproduces the C D curves of Fig. 5a for the same reduced velocities. Shaded vertical lines mark the frequencies for the drag minima to highlight that these coincide with the aforementioned wake-foil resonance.

IV. CONCLUSIONS

We characterized the dynamics of a flexible foil performing pitch oscillations immerse in an uniform airflow. For small to moderate elastic deformations, scale laws can describe the relationship between the maximum of amplitude, and the corresponding forcing frequency of the foil and the reduced velocity. In the same sense, inertial drag resistive contribution importance is depicted through a Cauchy number.

Study of the flow main structures was performed through PIV measurements and applying clustering methods. Propulsion forces were determined from velocities by means of an integral momentum scheme. Maximum propulsion was not obtained for the forcing and elastic resonance condition f + = 1 but rather when the wake resonant frequency is tuned with the foil dynamics. A convective instability associated with a jet-like structure that develops in the wake was studied by hydrodynamic linear stability analysis. Reduced velocity, forcing oscillations and the corresponding elastic deformations condition the shape of such jet velocity profiles, defining the resulting wake global resonant frequency. In agreement with these findings, we observe that vortex trajectories corresponding to propulsion maxima are more stable and remain parallel in the reverse Bénard-von Kármán configuration. FIG. 9: Schematic representation of the necessary steps involved in the linear stability analysis: (a) Time-averaged net mean velocity flow contours for u x with three profiles (circle symbols) fitted by eq. (A2) . (b) A mean velocity profile at ⟨U y (y, x = x c )⟩, for a fixed x-coordinate x c . Solving the Rayleigh equation for this profile, for a range of wavenumbers k = k r + ik i produces (d) a map for k i = constant in the ω-plane. At a critical point, a cusp on the curve, the point (ω 0r , ω 0i ) is determined for x c . In (c) the results for each x-position (symbols) are plotted in a ω-plane. Criteria for the prediction of a global frequency in a spatially developing flow is given by the condition ∂ω ∂x x=xs = 0 [START_REF] Chomaz | A frequency selection criterion in spatially developing flows[END_REF]. The saddle point ω s = ω(x = x s ) (big black circle) is obtained through an analytic continuation of the complex function ω 0 (x) considering x = x r + ix i (lines). ω s corresponds to the appearance of a cusp point. (dashed curve) [START_REF] Hammond | Global dynamics of symmetric and asymmetric wakes[END_REF]. Considering all the cases, the points ω s = (ω, σ) provide the data for Fig. 8 A generalized eigenvalue problem solves (A1)

ωB(k)ψ = A(k)ψ
We use NumPy linear algebra libraries in order to calculatie the pseudoinverse of B, for a k = k r + ik i and then computing the most unstables eigenvalues ω 0 = ω 0r + iω 0i . The global frequency for the spatially developing flow is determined using the criterion proposed by [START_REF] Chomaz | A frequency selection criterion in spatially developing flows[END_REF] which consists in finding a saddle point, ∂ω 0 ∂x x=xs = 0 of the complex function ω 0 (x) through use of the Cauchy-Riemann equations and analytic continuation to complex values of x = x r + ix i . Figure 9 summarizes the procedure, Hence, we could select the onset frequency f = ω 0 (x s ) 2π for this spatially developing flow to characterize the wake resonance. ωL U ∞ -7.00 -5.00 -3.00 -1.00 1.00 3.00 5.00 7.00 (2) Each ωi is assigned to the k-centroid closest according to D k i .

( (2) Each ωi is assigned to the k-centroid closest according to D k i .

( FIG. 10: Left: Steps for the construction of Clusters. Right: application example from our PIC data. Three diferent snapshots of voricity fields corresponding to the same cluster contributes to its determination through ensamble average.

Clustering

A clustering technique, named k-means algorithm [START_REF] Steinhaus | Sur la division des corps matériels en parties[END_REF][START_REF] Lloyd | Least squares quantization in pcm[END_REF], was applied in order to obtain clear and representative states of the flow. Burkardt et al. [START_REF] Burkardt | Pod and cvt-based reduced-order modeling of navier-stokes flows[END_REF] used the method for the first time in a fluid mechanics context. Later, Kaiser et al. [START_REF] Kaiser | Cluster-based reduced-order modelling of a mixing layer[END_REF] developed Cluster-based Reduced Order Modeling (CROM) which constitutes and alternative to classical POD-Galerkin reduced model schemes.

Strictly, k-means is a classification algorithm which minimizes the average distance between a specified number N c of points (called centroids) and the data in a phase space. For example, a point ω i in the phase space representing one PIV snapshot can be defined by its vorticity values in the spatial domain of interest, in our case, x ∈ [1 3.5], y ∈ [-1 1D]. Other components of the flow field may be chosen or even projections on some suitable base functions. Taking into considerations our PIV measurements, given their spatial resolution, ω i corresponds to 10836 values. We can use i.e. the euclidean norm d(ω i , ω j ) = ∑ 10836 s=1 (ω i (s) -ω j (s)) 2 as a proper distance between two points of the phase space.

The N c centroids c k , k ∈ [1, N c ] can be initialized by picking N c random PIV fields. Figure 10 presents the algorithm calculations and illustrate how coherent structures are clearer and more identifiable.

When the method is stopped, N c centroids are obtained, each one of them representing average states issued from the data and as separated as possible one from the other. A representation of a characteristic state of the flow can then be obtained by averaging all snapshots belonging to the same cluster. This method can be considered as a generalized phase averaging technique without the necessity of a strict periodicity of the phenomenon. An object oriented, open source k-mean algorithm is available in the Scikit-Learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] Python library .

FIG. 1 :

 1 FIG. 1: (a) General view of the experimental setup of the pitching flapping foil in the wind tunnel. A laser sheet defines the field of vision. (b) 2D view at midspan describing the flow configuration and geometrical description. (c) PIV instantaneous vorticity field corresponding to U R = 0.29, f + = 1.47. (d) Contours of the time-averaged net velocity flow u(x) = u x -U showing a jet-like flow structure.

FIG. 2 :

 2 FIG.2: a) Amplitude of non-dimensional displacements δ + as a function of the non-dimensional forcing frequency f + under various reduced velocities U R . The maxima of (f + max , δ + max ) are determined for each curve and plotted on b), showing the amplitude attenuation and frequency shift. In the inset, a quasi-linear dependence with respect toU R = U (ω n L) is shown. c) Renormalized version of a) δ + δ + max vs. f + f + max . d)Strouhal numbers defined from rigid and from total elastic displacements. The inset shows a scaling through the Cauchy-like number C * Y .

FIG. 5 :

 5 FIG. 5: a) Drag Coefficients as a function of the non-dimensional forcing frequency f + for different values of the reduced velocity U R . b) Corresponding Trefftz lengths. In cases marked by black circles, U T > PIV domain and the force computation with eq. (8) cannot be performed.

47 FIG. 7 :

 477 FIG. 7: Contours of vorticity for a given phase of the flapping period at U R = 0.29 at three different forcing frequencies. White contours shows the path of vortices during the whole period of oscillation. Full animated trajectories are available in Supplementary Material.

( 1 )

 1 Compute the distance matrix D k i : D k i = d(ωi, c

)

  For an ensemble W of Ns asigned snaphsots to a cluster, calculate the barycenter as b k = ∑ i ωi Ns, i ∈ W .

( 4 )( 4 )( 1 )

 441 Centroids are displaced to new barycenters c (1) k = b k . If c (n) Compute the distance matrix D k i : D k i = d(ωi, c

)

  For an ensemble W of Ns asigned snaphsots to a cluster, calculate the barycenter as b k = ∑ i ωi Ns, i ∈ W .

( 4 )( 4 )<

 44 Centroids are displaced to new barycenters c (1) k = b k . If c (n) tol no centroid is displaced, the method is stopped, otherwise it is iterated from (1).

  Profiles of the net mean velocity ⟨ux⟩ -U at y = 0 for U R = 0.29 and different forcing frequencies f + . Profiles of the intensity of velocity fluctuations ⟨u ′2 y ⟩ 1 2 at y = 0 for U R = 0.29 and different forcing frequencies f + .

		f +	0.98	1.14	1.30		1.47	1.63	U R	0.29	0.37 0.46 0.54 0.59
		1							1
		0.8							0.8
	ux -U	0.4 0.6						ux -U	0.4 0.6
		0.2							0.2
		0							0
	1 (a) 1 1.5 2 2.5 3 3.5 4 x L	1.1	1.2	1.3 f +	1.4	1.5	1.6
		f +	0.98	1.14	1.30		1.47	1.63
									UR	0.29	0.37 0.46 0.54 0.59
		1
									1
		0.8
									0.8
	⟨u ′2 y ⟩ 1 2	0.4 0.6						⟨u ′2 y ⟩ 1 2	0.4 0.6
		0.2
									0.2
		0
		1	1.5	2	2.5	3	3.5	4	0
	x L (c) 1	1.1	1.2	1.3 f +	1.4	1.5	1.6
									(d) Maxima of velocity fluctuations profiles ⟨u ′2 y

(b) Maxima of velocity profiles ⟨ux(y = 0)⟩ -U for each (f + , U R ) case.

  The bottom plot presents Strouhal numbers (St) calculated in the wakes through linear stability as a function of f + for U R =0.29; 0.37.

	CD	-0.6 -0.4 -0.2
		0.25	UR 0.29
			0.37
		0.2
	St
		0.15
		0.1
		1	1.2	f +	1.4	1.6
	FIG. 8:

We acknowledge support of the French-Argentinian International Research Project IVMF, CNRS-INSIS (France), CONICET (Argentina).

Wake flow stability properties are surveyed through a linear inviscid 1D approach, using the Rayleigh equation:

being U y a time-averaged velocity profile at a fixed x-coordinate, k the wavenumber, ω the frequency and ψ the stream function.

An analytical fit for U y to avoid numerical noise is used, as in the study by [START_REF] Triantafyllou | On the formation of vortex streets behind stationary cylinders[END_REF] for cylinder wakes:

The boundary conditions are on the one hand, symmetry about y = 0, corresponding to BvK vortex streets. This also is the case of reverse BvK, which applies directly to our measurements. On the other hand, for y ≫ 1, U y ≃ U the free flow velocity, then ψ → 0.