Experimental Resonances in Viscoelastic Microfluidics

Pulsatile flows of viscoelastic fluids are very important for lab-on-a-chip devices, because most biofluids have viscoelastic character and respond distinctively to different periodic forcing. They are also very important for organ-on-a-chip devices, where the natural mechanical conditions of cells are emulated. The resonance frequency of a fluid refers to a particular pulsatile periodicity of the pressure gradient that maximizes the amplitude of flow velocity. For viscoelastic fluids, this one has been measured experimentally only at macroscales, since fine tuning of rheological properties and system size is needed to observe it at microscales. We study the dynamics of a pulsatile (zero-mean flow) fluid slug formed by a viscoelastic fluid bounded by two air-fluid interfaces, in a microchannel of polymethyl methacrylate. We drive the fluid slug by a single-mode periodic pressure drop, imposed by a piezoactuator. We use three biocompatible polymer solutions of polyethylene oxide as model viscoelastic fluids, and find resonances. We propose a model accounting for surface tension and fluid viscoelasticity that has an excellent agreement with our experimental findings. It also provides an alternative way of measuring relaxation times. We validate the method with parameters reported in the literature for two of the solutions, and estimate the relaxation time for the third one.


INTRODUCTION
The study of oscillatory fluid flow at microscales has become relevant due to the increasing number of applications that use this type of motion. For example: chemical synthesis inside microfluidic channels [1], liquid-liquid extraction [2], mixing by oscillatory cross flow [3][4][5][6][7], cooling of microelectronic circuits by micro oscillating heat pipes [8], inertial focusing of particles of a few microns [9,10], DNA elongation studies [11] and studies of oscillatory movement of liquid plugs displaced by air in microchannels as model pulmonary flows [12,13].
Pulsatile flows of viscoelastic fluids are very important for most organ-on-a-chip devices, where the natural mechanical conditions of cells are emulated [14][15][16][17], since most natural processes occur at certain characteristic frequencies. The characterization of viscoelastic fluids under non-steady pressure forcing is also important for lab-on-a-chip clinical analysis of biofluids such as blood, mucus, or synovial fluid. The dynamics of polymeric viscoelastic solutions under pulsatile forcing in microchannels is an area of recent development [18]. Flow of these solutions is strongly influenced by chemical properties of the polymer, its molecular weight and ramifications, concentration, the nature of the solvent, temperature and pressure [19].
The fluid response to an oscillatory pressure gradient has often been described by the dynamic permeability, a theoretical linear response function that has been obtained for numerous confined fluids: Newtonian, Maxwellian and general linear viscoelastic fluids, in a wide range of confining geometries [20][21][22][23][24][25][26][27]. It has also been obtained theoretically for Newtonian and viscoelastic fluids confined in elastomeric materials at microscales [28,29] and for compressible binary fluids [30]. A distinctive feature of the dynamic permeability, when elastic elements are present in the system, is that it presents resonances, which refer to particular pulsatile periodicities of the pressure gradient that maximize the amplitude of fluid velocity. Experimental observation of resonances consists of an increase of flow velocity amplitude at a specific frequency range of the driving pressure gradient, that maximizes the momentum transfer to the fluid. For single fluids, resonances have only been reported experimentally at macroscales [31,32], since fine tuning of rheological properties and system size is needed to observe them at microscales, in a desired frequency range.
Recently, a model to study the dynamics of a pulsatile (zeromean flow) fluid slug, consisting of a Newtonian fluid and two air-fluid interfaces, driven by a periodic pressure gradient in a rectangular microchannel, has been proposed [33]. In that model, a stress tensor for a Newtonian fluid, together with Laplace condition for the pressure jump at both sides of the curved air-fluid interfaces, has been considered. Analytical solution of the model showed, for relatively low frequencies, a monotonic increase with frequency of the magnitude of the dynamic permeability as well as the emergence of a resonant behavior, due to the presence of surface tension. Microfluidic experiments were designed and implemented to observe both the low-frequency dynamics and the resonance. The model was then validated against the experimental results and used as a proposed strategy to measure surface tension in dynamic situations.
There are different ways to impose oscillatory frequencies to a fluid inside a microchannel. By using syringe pumps, low frequencies of oscillation (below 10 Hz) are achieved [1]; in contrast, the use of high-speed valves and gas-pressurized fluids [3], mechanical motors [34], heating [35] or mechanical displacement of an air bubble [36] can increase the forcing frequency range to 10-1,000 Hz. Alternative possibilities to impose pulsatile forcing in this frequency interval are the use of a moving train of droplets [37] and the coupling of a loudspeaker diaphragm to a microfluidic chamber [4]. Finally, coupling the displacement of a piezoelectric to a fluid encompasses a wide range of forcing frequencies of the methods described previously [38].
There are several sophisticated theoretical models to study the rheological behavior of PEO solutions, in different ranges of concentration and molecular weights. Of particular importance are the Phan-Thien-Tanner (PTT) model [39,40] and the Cross model [19]. They have been adequate to study several experimental conditions and driving forces where a complex rheological response, involving elongational and shear thinning effects, has been experimentally observed and theoretically reproduced. However, there is also experimental evidence that a Maxwellian model predicts correctly and accurately the behavior of small ejected, low molecular weight PEO (1x10 6 g/mol) droplet jets [41]. Moreover, within microchannels of constant sectional area, several works suggest that for spatially-uniform pressure gradients elongational and shear thinning effects, like the ones considered by the PTT and Cross models, are irrelevant [40,42,43]. Furthermore, despite the fact that viscoelastic fluids generally involve several relaxation times, many studies of fluids with complex rheological behavior often report a single dominant Maxwellian-like relaxation time, fitted from their experimental data, since the Maxwell model is used as an archetype in the field.
In this work, we perform experimental and theoretical studies of the dynamics of a pulsatile (zero-mean flow) microfluidic slug, formed by a viscoelastic fluid bounded by two air-fluid interfaces in a rectangular microchannel, and find resonances in the dynamic permeability. We have driven the fluid slug by a single-mode periodic pressure drop, imposed by a piezoactuator in the range from 0.5 to 200 Hz, managing to keep the amplitude of the dynamic pressure drop practically constant at all frequencies. We have determined the displacement of the viscoelastic slug by visualization of the oscillatory movement of air-fluid interfaces. We have used three biocompatible polymer solutions of polyethylene oxide (PEO), as model viscoelastic fluids, because the rheological behavior of PEO has been widely assessed [19,44,45]. We propose a linear model accounting for surface tension and fluid viscoelasticity, that has a good qualitative agreement with all of our experimental findings and a quantitative agreement for low pressure drops, where the linear theory is expected to describe the system. Such agreement provides an alternative way of measuring relaxation times. We validate the method against parameters reported in the literature for PEO of two different molecular weights: 1x10 6 g/mol (PEO1) and 5x10 6 g/mol (PEO5); and estimate the relaxation time for PEO of 8x10 6 g/mol (PEO8). This is of great relevance because relaxation times are sometimes difficult to measure for low polymer concentration in conventional rheometers [45][46][47][48].
The paper is organized as follows: section 2 describes the experimental procedure and the data analysis; section 3 describes the experimental results including resonances of the dynamic permeability; section 4 introduces a theoretical model for pulsatile viscoelastic slugs; section 5 compares experimental results for the dynamic permeability with predictions obtained from the theoretical model, it also introduces a proposal to measure relaxation times; section 6 summarizes the most important conclusions and perspectives.

Microfluidic Device
We machined a straight microchannel (37.48 ± 0.11 mm long, 1.00 ± 0.04 mm wide and 0.31 ± 0.05 mm deep) on a 2 mm-thick polymethyl methacrylate (PMMA) plate using a CNC machine (CNC3018). The channel was sealed with a second PMMA plate with four inlets (Figure 1) exposing both parts to volatilized chloroform for 4 min and pressing them by a pair of slides and clamps. The bonding was completed by sonication of the device in ethanol at 50 • C for 15 min [50, 51].

Experimental Setup
A piezoelectric actuator equipped with flexural hinges as an amplifier device (APF705, Thorlabs) was attached on one of its sides to an elastic membrane that covered a rigid polyethylene cylinder (15 mm long, 4.7 mm diameter). The opposite side of the cylinder has a seal with a tubing (0.51 mm ID, 1.19 mm OD and 1 cm long; Microbore PTFE Tubing, Cole-Parmer) inserted in the middle. The other end of the tubing was introduced into the first microdevice inlet. The movement of the actuator displaces the air in the cylinder and transduces an oscillatory movement to a slug of PEO solution (1.0 cm in length; 3.1 µL) situated in the middle of the microchannel (Figure 1). The oscillation frequency and amplitude of the piezoelectric motion was controlled by a multifunction data input/output device (USB-6351, National Instruments) and magnified by a Trek PZD350A High-Voltage amplifier (75-150 V). The pressure drop was measured by a differential pressure sensor (Honeywell 142PC01G) attached by PTFE tubing to the second and third inlets of the microfluidic channel. The fourth inlet was open to the atmosphere. The displacement of the liquid slug was visualized with the aid of an inverted microscope (DM IL LED, Leica) and the movement of the interface closest to the atmosphere outlet (IF2 in Figure 1) was recorded with a high-speed camera (Phantom Miro M110, Vision Research). Depending on the driving frequency, videos from 30 fps up to 3,000 fps were acquired after a 10 s   stabilization period of cycling movement to ensure recording after transient states. The size of the fluid slug was verified after each measurement to confirm that no evaporation had occurred.

Data Analysis
The piezoactuator movement was adjusted by regulating the input voltage to keep the same reference pressure drop for all the range of frequencies studied. For each PEO solution, oscillatory pressure drops of four different amplitudes were used: 225, 450, 700, and 900 Pa, each in the frequency range from 0.5 to 200 Hz (Figure 2). The sinusoidal shape of pressure drop allowed us to fit a sinusoidal wave to obtain the amplitude of each signal.
The videos of the interface movement were analyzed using MATLAB utilities, that track the position of all interface points through time, then velocity was obtained by numerically differentiating position data.
To prove that the frequency imposed by the piezoactuator was consistent with the interface movement, Fourier transform of the pressure drop, interface displacement and interface velocity were performed. A dominant peak for the spectrum of all these signals was observed at the same frequency of oscillation of the piezoelectric transducer, indicating that the fluid slug follows the dynamics imposed by the piezoactuator (Figure 3).
The dynamic contact angle of the interface was determined from video image analysis of the advancing and receding time lapses. The interface profile at every time, was fitted to a fourthdegree polynomial function. The fit reproduces very well the interface profile and was extrapolated to compute the contact angle at the wall.

RESULTS
The air-fluid interfaces of the polymeric solutions display a characteristic curvature that flattens and bends in every oscillation cycle. This is illustrated in Figure 4 middle. In Figure 5, we show the air-fluid interface of PEO5 0.1% through an entire oscillation cycle at four different frequencies.
A maximum amplitude of the displacement, A max , for this experiment occurs at 40 Hz. Peak-to-peak amplitude of the displacement is highlighted in red so the frequency dependent behavior can be assessed visually. It is clear that the neighboring smaller or larger frequencies display a smaller displacement.
By changing the frequency but keeping constant the amplitude of pressure drop driving the movement, a clear non-monotonic behavior of the amplitude of the interface displacement is observed for different driving frequencies, even for the PEO of lowest molecular weight. Top panels in Figure 6 show that the highest displacement for each fluid increases in magnitude as the imposed pressure rises, which is an expected result for an increasing driving force. The maximum interface movement is observed for PEO1, the fluid with the smallest molecular weight and viscosity. The interface velocity, at the center of the microchannel, as a function of the oscillating frequency shows an asymmetric bell-shape curve for each imposed pressure (Figure 6 middle panels). As expected, the maximum velocity amplitude rises as the pressure increases for all PEO solutions. The peak of each curve is the resonance frequency, meaning that at this frequency the amplitude of flow velocity is maximum in the frequency range studied. A non-trivial effect is observed in which the resonance frequency decreases with an increasing pressure drop. This is part of the non-linear behavior of the system response. Since we are driving the system with a onemode pressure drop, we expect that, at least in the linear regime, the amplitude of the interface velocity would be given by the maximum amplitude of the interface displacement multiplied by its frequency. For this reason, the resonances observed in velocity have higher frequencies than those obtained for the amplitude of interface position.
We also analyzed the dynamic permeability of each polymer solution as a function of frequency and amplitude of pressure drop (Figure 6 bottom panels). We present an operational definition of the amplitude of the local dynamic permeability at the center of the channel, as the ratio between the maximum amplitude of the velocity at the center of the channel, divided by the pressure gradient-given by the quotient of pressure drop, p, and slug length, L,-that is, Derivation of Equation (1) will be given in section 4 Model for small pressure drop values; however, we will use this operational definition for all pressure drops studied, since it is a convenient way to cancel out the small differences in pressure drop amplitude for the different frequencies tested (see Figure 2). In this sense, K can be interpreted as a velocity rescaled by a pressure gradient. Accordingly, in the bottom panels of Figure 6, we observe that the resonance frequencies of the local permeability are roughly the same as those obtained for velocity data. Also, as pressure drop increases, the resonance frequency of the permeability decreases. This effect is more pronounced for low molecular weights. For the dynamic permeability value at resonance, there is no clear trend when the amplitude of the pressure drop changes. For details and graphs describing this behavior, see Supplementary Material. Regarding the dynamics of the contact angle, we found that advancing angles are larger than receding ones. As an example, Figure 7 shows the dynamic angles obtained when the amplitude of pressure drop is 225 Pa. It has been reported in the literature that, when velocity increases, the advancing angle augments and the receding angle decreases, so a larger difference between them should be observed [52][53][54][55]. The difference of advancing and receding angles, θ = θ adv − θ rec (sometimes called contact angle hysteresis [52]), is shown as a function of frequency in Figure 8. For a single-mode oscillatory flow, the fluid velocity increases with frequency up to the resonance and then decreases. We can observe that the same trend exists for θ as a function of frequency. The fact that the dynamic contact angle difference is affected by the interface velocity, has previously been reported for capillary numbers close to the ones of our experiments (Ca = 10 −5 to 10 −3 ) [53,55]. This phenomenon is attributed to surface roughness and chemical heterogeneity [56,57], but there is an influence of the fluid rheological properties, like shear thinning or elasticity [53]. In our slug, an important component of elasticity is given by the presence of two interfaces.

MODEL
In order to explain theoretically our experimental results, we build up a model containing two basic features: the viscoelastic character of the fluid and the elastic character of interfaces. The  interplay between these two elasticities will lead to a complex behavior, the simpler the rheology of the viscoelastic fluid, the easier the understanding of the physical interaction. Because of this, we present a model for the uniaxial dynamics of a viscoelastic slug, taking into account the presence of interfaces and the viscoelastic character of a single-relaxation-time fluid, that could be extended to models containing more characteristic times. We build up a model for the dynamics of viscoelastic slugs over a model presented for the dynamics of a Newtonian slug [33].
A study of a pulsatile fluid slug consisting of a Newtonian fluid and two air-fluid interfaces driven by a periodic pressure gradient, has been recently proposed and validated [33]. In that model, a stress tensor for a Newtonian fluid of the form τ = −η∇v, together with Laplace condition for the pressure jump at both sides of the air-fluid curved interfaces, p = p driving + σ κ 1 + σ κ 2 , has been considered. In these expressions, η is the fluid viscosity, v is the axial fluid velocity, σ is the surface tension of the air-fluid interfaces, κ 1 and κ 2 are the left and right hand side curvatures, respectively, and p driving is the pressure drop external to the fluid slug (on the air side). The dynamics for such Newtonian slugs is described by an integro-differential equation in space and time, which, in frequency domain, can be written as a simple equation, differential in space and algebraic in frequency, that reads: where ω denotes angular frequency, ρ is the fluid density, L is the length of the fluid slug, andv andp denote Fourier transforms of velocity and pressure, respectively. This equation incorporates momentum conservation, the stress tensor for a Newtonian fluid, Laplace equation for the pressure jump at the interfaces, an approximation of interface curvatures as concavities, and continuity of velocities at both interfaces. This model has given a correct description of the experimental dynamics of a water slug and of a 70% glycerol solution in water slug, when interfacial curvatures are considered to be a response to a dynamic external pressure gradient p driving /L. Details of the derivation of Equation (2) can be seen in [33]. A Newtonian slug stress tensor, that integrates stresses of the Newtonian fluid and the interfaces, of the form substituted in the linearized momentum conservation equation for uniaxial flow in the x direction, gives exactly Equation (2) for the dynamics of a Newtonian slug. A stress tensor of the form (3) for a material consisting of a volume of fluid and two air-fluid interfaces, has not been introduced in the literature, to the best of our knowledge, since classical treatments describe both fluid phases and apply boundary conditions at air-fluid interfaces, rather than describing the system fluid-interfaces as a composite material. Vazquez-Vergara et al. [33] together with discussion in the previous paragraph, show that introduction of a slug stress tensor, as the one in Equation (3), is a consistent approach to describe the zero-mean flow, linear pulsatile dynamics of Newtonian slugs.
Viscoelastic fluids might involve, in general, several relaxation times. The coupling of such times with the characteristic time given by the presence of interfaces, is expected to lead to a complex dynamics of viscoelastic slugs. To understand such coupling, we start by introducing the simplest model of a viscoelastic slug, consisting of a volume of linear Maxwellian fluid, that has a single relaxation time, and two air-fluid interfaces, with surface tension σ . 1 We propose the following expression for the viscoelastic slug stress tensor: where the parameter t r is the Maxwell relaxation time. Equation (5) reduces to the constitutive equation of a Newtonian slug (Equation 3) in the limit t r → 0, and reduces to the Maxwell model in the absence of interfaces [58], that is, in the limit σ → 0. Along with the previous consistency proofs, experimental validity of the model given by Equation (5) Solution of Equation (6) subject to no-slip boundary conditions, for flow in a rectangular microchannel whose plates are separated by a distance 2l, gives a linear relation between velocity and pressure drop in frequency domain, that is, v(z, ω) = −K(z, ω) p air L , where z is the coordinate perpendicular to the plates. The complex local dynamic permeability, K(z, ω), is given by It is worth noticing that Equations (7) and (8), which are the solution for the dynamics of a Maxwellian slug given by our model, are consistent with the pulsatile solution of the linear Maxwell model for a single fluid in a rectangular cell in the limit of zero surface tension given in [26] (equivalent to the solution in [23,24,59] in the cylindrical case).
Expressions (7) and (8) are valid for general periodic timedependent pressure drops, consisting of an arbitrary number of sinusoidal modes. In particular, for a one-mode driving pressure tensor proposed for the fluid slug would beτ OB ∇v.
drop of frequency ω 0 , it can be shown that, in time domain, the amplitude of the velocity at the center of the cell, v max , is related to the amplitude of the pressure drop, p max , as v max = K(z = 0, ω 0 ) p max L , (9) which is in agreement with the operational definition of K, used to compute the permeability from experimental data in Equation (1). Our model has been deliberately developed for zero-mean pulsatile flows, due to the oscillatory nature of our experimental driving force. It therefore cannot be used to model fronts in imbibition-like systems, where the pressure drop has always the same sign and the interface curvature is due to wetting. In our case, curvature effects due to wetting cancel out since they have opposite signs on the left and right side interfaces [33], and the dynamics of the slug will be governed by the instantaneous interfacial curvatures caused by pulsatile forcing.

COMPARISON WITH EXPERIMENTAL RESULTS
We compare the experimental results for K at low-amplitude pressure drops (225 Pa) from Figure 6 with K(0, ω 0 ) derived from the linear model developed in the previous section. For the figures, we simply use K to denote the local dynamic permeability at the center of the microchannel. Figure 9 shows experimental and theoretical predictions of K as a function of the driving frequency for three PEOs. A log-log scale has been used to highlight the tendency of low frequency data. Figure 9 left shows, with green dots, the permeability K obtained from experimental data for PEO1. With a red continuous line, it shows the theoretical permeability, as predicted by our model for a Maxwellian slug, obtained from Equation (8) with z = 0 and a relaxation time, t r = 1.78 ms, reported as Maxwellian in the literature [41]. As reference, we have also plotted, in a blue continuous line, the theoretical permeability for a Newtonian slug. As Figure 9 left shows, the agreement between experimental data and theoretical prediction for a slug of a viscoelastic fluid obeying Equation (8) is excellent, both, at low frequencies and around resonance; Figure 9 middle shows equivalent curves for PEO5. Since the relaxation time, at the concentration used in our experiments, is not reported in the literature, we took one reported for PEO4, as surrogate [19]. The agreement between the green dots, obtained from experimental data, and the red line predicted by our model, is very good, both, at small frequencies and around resonance, despite the fact that the relaxation time used was obtained from a fit to a Cross model [19,60]. Before discussing Figure 9 right for PEO8, for which there is no relaxation time reported in the literature, at the desired concentration, we will discuss the theoretical behavior of the resonance frequency, in terms of characteristic frequencies of the system.
We can define three characteristic frequencies of the system that depend on viscosity, η, surface tension, σ , relaxation time, t r , and system's geometry, l (half the microchannel thickness), as In terms of these frequencies, the argument of the cosine term, Al, in Equation (8), can be written as We find two different regimes for the resonance frequency. For ω relax ≪ ω η , the resonance frequency, ω res , is given by while for ω relax ≫ ω η , it is given by The resonance frequency, with the proper scaling to make the second of these regimes collapse, is given in Figure 10.
We have found numerically the maximum of the local dynamic permeability in Equation (8) using characteristic values of ω η and ω σ for each of the three PEO solutions employed in the experiments. We have plotted the resonance frequency as a function of the independent variable ω relax /ω η . The red curve corresponds to PEO1, the blue curve to PEO5 and the gray curve to PEO8. Red and blue dots over the curves (for PEO1 and for PEO5, respectively) represent pairs (t lit r , ω res (t lit r )), in the corresponding rescaled units of the graph, where ω res (t lit r ) is the theoretical resonance obtained using a relaxation time from the literature.
Our resonance curves in Figure 10 could serve "as a rheometer" to validate or estimate viscoelastic relaxation times. Such a procedure is schematized with horizontal dashed lines in Figure 10, which relate a reasonable resonance frequency range-obtained by visual inspection of Figure 9-, with a range of possible relaxation times, obtained from the vertical dashed lines, through the theoretical curve in Figure 10. To validate our method, we have estimated that the resonance FIGURE 10 | Illustration for the way of estimating the value of relaxation time (through ω relax ) by visual estimation of the resonance frequency, and its corresponding error bars. The theoretical model predicts a relation between resonance frequency and relaxation frequency that can be exploited to validate or estimate relaxation times. We have used PEO1 solution (red dashed lines), and PEO5 solution (blue dashed lines) to validate the method. Then we have used the method to estimate the relaxation of the PEO8 solution, from its experimental resonance (gray dashed lines). Red and blue intervals, marked with dashed lines, are in agreement with theoretical data of resonance frequency and relaxation times (placed as a circle red dot) for PEO1 and (placed as a circle blue dot) for PEO5. Black dashed line represents an approximated relation for ω res in Equation (13). Continuous lines are exact results for resonance, obtained by numerical means.
frequency for PEO1, is in the range [120:136] Hz, and in the range [100:120] Hz for PEO5. As Table 1 shows, the range of relaxation times estimated by the theory-from measurement of the experimental resonance frequency-are in agreement with the values for the relaxation time reported in the literature. We consider that, as a proof of concept, this validates our method for estimating relaxation times of single-relaxation-time viscoelastic fluids. Accordingly, we estimate from Figure 9 that, for PEO8, the resonance frequency is in the range [80:100] Hz, and estimate that the relaxation time would be in the range [6.8:14.8] ms. Even though this range is wide, when considering a value in the middle of this range, we obtain the red curve in Figure 9 for the dynamic permeability,which, despite being less accurate (a) Taken from [41]. (b) Taken from [19].
than the ones for the lower molecular weight polymers, still gives a reasonable agreement for both, the tendency of the experimental permeability at low frequencies, and the region around resonance. It is important to note that the relaxation times obtained with our method follow the trend that, the larger the molecular weight of PEO, the larger the relaxation time. This is in qualitative agreement with global trends observed in literature for a vast range of concentrations and molecular weights in PEO aqueous solutions [19,41,61]. Our method promises to be valuable for low molecular weight polymers, for which relaxation times are difficult to obtain experimentally by conventional means. For high molecular weight polymers, more sophisticated models, including shear thinning, might result necessary to describe the dynamics and to obtain rheological parameters from it, as it happens for high concentrations of low molecular weight polymer solutions [19,45]. It is worth mentioning that the size of the fluid slug can be experimentally adjusted to fall in the range where the approximated black dashed line in Figure 10 is valid, that is, whenever it agrees with the exact value of resonance frequency, given by the continuous lines; at such regime, there is a simple analytical relation between resonance frequency and relaxation time, and the indirect determination of relaxation times could be carried out easily using Equation (13). In contrast, the plateau observed for each continuous line on the right side of Figure 10, corresponds to the regime where the resonance frequency is given by Equation (12), when the relaxation time is not relevant. The curves differ from each other due to their viscous frequency, ω η . Our theory is also able to predict the contact angle, given an initial shape of the interface. Several models for front dynamics in the literature establish a difference between a static contact angle at equilibrium situations, and the angle observed due to an imbibition-like front dynamics, where the pressure gradient has always the same sign, this is typically called a dynamic contact angle. In addition to such descriptions, recent studies have dealt with a different time-varying dynamic contact angle, which is affected not only by imbibition phenomena but also by oscillatory driving forces [62]. This is the case of our dynamic contact angle. We explain in a nutshell how do we compute dynamic contact angles from our theoretical model: our differential equation gives the slug velocity v(z, t) as a function of pressure drop. Since interface shape, u(z, t), and velocity are related through ∂u/∂t = v, once the velocity v(z, t) is known, we can integrate this equation to obtain u(z, t) = v(z, t)dt + u 0 (z), where u 0 (z) is an integration constant which gives the interface profile at rest (or, equivalently, at very high frequencies). Once the interfacial profile is known in time, the arctangent of its slope close to the wall gives the dynamic contact angle.
It is worth mentioning that our experiments measure an angle along the channel width, not along the channel height, which is the dimension modeled in Equation (6). So, it is necessary to find out a relation between the angles measured in both planes. We follow Tabeling results [63] on steady flow where the relation between the slope of flow velocity at the channel walls in both planes obeys a linear relation of the form ∂v ∂y y=W/2 = m ∂v ∂z z=l (14) where m is a factor that only depends on the aspect ratio W/2l (see Figure 11 left). Since for a single-mode oscillatory flow, interface shapes and velocity are linearly related through the driving frequency, that is, ω 0 , as v(y, z, t) = ω 0 u(y, z, t), an equivalent relation can be proposed for interfacial profiles, u, as ∂u ∂y y=W/2 = m ∂u ∂z z=l (15) This geometrical relation between the slope at both confining dimensions is illustrated in Figure 11 left. The dynamics of both angles is illustrated in Figure 11 right. With the correction explained above, the experimentally obtained contact angles are compared to the ones predicted by our model. This is shown in Figure 12. We find that the contact angle predicted by the theory is larger than the one obtained from experiments; however, our theory correctly predicts the phase difference between the angle and the pressure gradient.
A model for single Maxwellian fluids accounting for channel width and height, has found that the effect of the second confining dimension is relevant only at high frequencies [64]. We therefore consider that if such analysis were extended to Maxwellian slugs, it would not affect the conclusions regarding resonances and dynamic permeabilities presented in this work.

CONCLUSIONS AND PERSPECTIVES
We have made a thorough experimental study of the dynamic behavior of pulsatile fluid slugs made by three biocompatible viscoelastic fluids [65]. We have studied their responses in a wide frequency range, from [0. 5:200] Hz, at four different amplitudes of the pressure drop, which have been maintained practically constant in all the frequency range. We have chosen the dynamic permeability as a parameter to characterize the fluid dynamics since, in the regime where flow and pressure drop are linearly related, it can be analytically demonstrated that it is a response function of the system to a pulsatile pressure drop. At higher pressure drops, it is a convenient way to represent a rescaled velocity, which cancels, to linear order, the small differences in pressure drop amplitude applied at different frequencies. We have found that the permeability of the three viscoelastic slugs present resonances, that is, a special frequency range of the FIGURE 11 | Left: Illustration of an instantaneous velocity profile in a cross section of a rectangular microchannel, computed theoretically following [63]. Right: Illustration of the time behavior of the angle (given by the arctan of the slope) along the channel width (red), and along the channel height (blue), for a channel whose width is three times its height. FIGURE 12 | Comparison between predicted and measured contact angles for a slug of PEO1 solution pulsated at frequency of 120 Hz, and a pressure drop amplitude of 225 Pa. Blue and red lines correspond to the contact angle, obtained from experimental measurements, at the left and right walls of the microchannel, respectively; black line corresponds to the theoretical prediction. Agreement between theoretical and experimental results is noticeable, particularly because their oscillation is in phase. oscillatory pressure gradient that maximizes the amplitude of flow velocity.
We have also developed a continuum-mechanics linear model of viscoelastic slugs, containing both, the elasticity of the fluid and the elasticity given by air-fluid interfaces, through surface tension. The slug model gives excellent agreement with our experimental results at low-amplitude pressure drops. Such an agreement was a necessary condition for model validation.
The dynamic permeability at all frequencies coincides very well for PEO1 and PEO5 solutions, both, at small frequencies and at resonance. Coincidence of experimental and theoretical resonances provides an alternative strategy for measuring relaxation times. We validated this strategy with relaxation times reported in the literature for PEO1 and PEO5 solutions, and estimated the relaxation time for PEO8. With such estimation, the dynamic permeability for PEO8 gives the correct slope as a function of frequency (in log-log scale), at small frequencies, and gives a very good agreement around resonance. It is worth mentioning that viscoelastic fluids, in general, have several relaxation times. If one uses a single-relaxationtime model, this one should correspond approximately to the larger characteristic time experimentally observed. Our model for viscoelastic slugs could be extended, in a more or less straightforward manner, to models with several relaxation times. The problem for validation in that case, would be the lack of experimental information in the literature of the several relaxation times reported for a specific molecular weight and a specific concentration of a polymer solution obeying a specific model. Finally, we have compared the dynamics of the contact angle, and found that theory and experiments predict similar amplitudes and exactly the same phase shift with respect to the oscillatory pressure drop. In our experiments, the interface motion is negligible close to the wall, when compared to the motion of the interface at the center of the channel. For this reason, we have not included slip in our model and the contact line displacement is not a concern for our relatively smallamplitude pressure drops. This is the first time that experimental resonances of the dynamic permeability of viscoelastic fluids, confined at microscales, are reported in the literature. Our results are relevant for flow and shear rate control, with potential applications to many systems, like organ-on-a-chip devices where the natural mechanical conditions of cells are emulated [14][15][16][17]; and biofluids, which are typically viscoelastic, are present. Our results could also be useful to study how cells would respond to different imposed, non-physiological, external stresses [66][67][68][69][70][71].

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