On the linear viscoelastic behavior of semidilute polydisperse bubble suspensions in Newtonian media

In this work, we investigated the linear viscoelasticity of semidilute polydisperse bubble suspensions via small amplitude oscillatory shear (SAOS) tests performed in a rheo-optical setup. For all tested suspensions, the measured viscoelastic moduli ( G 0 , G 00 ) aligned with the theoretical predictions of the Jeffreys model for average dynamic capillary numbers ( h Cd i ) greater than unity. But at lower h Cd i values, experimental G 0 values exceeded theoretical predictions. To investigate this, we considered the effects of suspension polydispersity and various SAOS measurement artifacts, including bubble rise, coalescence, and changes in suspension microstructure over time. Polydispersity could not cause the observed deviation, because the viscoelastic trends deviate from those of classic single-mode relaxation only for bimodal bubble size distributions with equal volume fractions of very small and very large bubbles; in any other case, the polydisperse suspension behaves as monodis-perse with a bubble radius equal to the volume-weighted mean radius. Furthermore, bubble rise proved to play a minor role, while SAOS rheo-optical experiments revealed that bubble size and organization varied negligibly during our measurements. The G 0 deviation at low h Cd i values was linked to bubble fluid dynamic interactions induced by the bubble spatial distribution. Image analysis showed that at low bubble volume fractions, stronger and prolonged preshearing reduces these interactions by increasing the average interbubble distance. But this effect is negligible in denser suspensions, which show similar G 0 trends for any applied preshearing. Finally, a multimode Jeffreys model fitted to the experimental data showed that bubble interactions complicate the relaxation process, introducing multiple relaxation modes. © 2024 Author(s). All article content, except


I. INTRODUCTION
Bubble suspensions are ubiquitous in many industrial applications, from the manufacturing of cosmetics and pharmaceuticals [1] to the production of foamed cement [2].In all these contexts, bubbles can alter considerably the texture and processability of the products [3], a fact that highlights the need for fundamental investigations on the rheology of these systems.
Bubbles have been shown to induce shear-thinning and other viscoelastic phenomena, even when dispersed in Newtonian ambient fluids [4][5][6][7][8][9][10][11].The relevance of these phenomena rises with the bubble volume fraction, f.In addition to this parameter, two dimensionless numbers characterize the rheology of bubble suspensions: the capillary number and the dynamic capillary number.The former, denoted as Ca, characterizes the deformation of bubbles in equilibrium, relating it to their relaxation time, λ, and the applied shear rate, _ γ , where η s is the viscosity of the ambient fluid, R is the radius of the relaxed nondeformed bubbles (assumed to be all identical), and σ α,β is the surface tension of the ambient fluid in air [4,6,7].For vanishingly small values of Ca, the suspended bubbles essentially do not deform, retaining their original spherical shape; as the applied shear rate and, in turn, Ca increases, the bubbles start deforming in the flow direction, leading to the shear-thinning behavior observed in bubble suspensions.For polymeric liquids, the equivalent dimensionless number is the Weissenberg number (Wi) [12,13].
The dynamic capillary number, denoted as Cd, describes the steadiness of the flow (from a Lagrangian point of view) and the viscoelastic character of the bubble suspension, comparing the relaxation time of the bubbles, λ, with the time scale over which the shear rate changes significantly (or, more generally, over which "a typical fluid element experiences a significant sequence of kinematic events" [14]).For polymeric liquids, the equivalent dimensionless number is the Deborah number (De) [12,13].For a simple oscillatory shear flow, the time scale of the flow coincides with the inverse of the oscillation frequency, ω, so that hand, for Cd ) 1, the bubbles do not relax fast enough to reach their equilibrium configuration and, therefore, behave as elastic solids.In this case, the flow is unsteady, and the bubble suspension reveals viscoelasticity, the elasticity arising from the bubble interfaces [15].For monodisperse bubble suspensions, the onset of viscoelasticity is observed at Cd 1 [4].
The steady-shear rheology of bubble suspensions, in the limit of Cd !0 for varying Ca, has been thoroughly investigated for monodisperse [4,11,16] and polydisperse [5,6,8,10] systems, with studies showing that bubble deformation [4,7,9] and clustering [8] are responsible for the shearthinning behavior.Time-dependent flows (from a Lagrangian point of view) with infinitesimal bubble deformation, that is, in the limit of Ca ! 0 for varying Cd, have also been investigated to characterize the linear viscoelastic behavior of bubble suspensions [4,10,17,18].These flow conditions can be achieved experimentally via small amplitude oscillatory shear (SAOS) rheological tests, where a sinusoidal deformation γ ij is applied to the system and the resulting shear stress τ ij is measured, where γ 0 ij is the amplitude of the applied deformation, ω is the oscillation frequency, G 0 is the elastic modulus, which describes the elastic character of the suspension, and G 00 is the loss modulus, which describes the viscous character of the suspension.Below, we discuss the main models that describe the linear viscoelastic behavior of bubble suspensions as a function of the bubble volume fraction and the dynamic capillary number.
Llewellin et al. [4] suggested a constitutive equation for dilute monodisperse bubble suspensions, expressed in the form of the linear Jeffreys model [19].The equation holds for Ca ( 1 and varying Cd and has its theoretical foundation in the analysis of Frankel and Acrivos [20] on the timedependent flow of dilute monodisperse emulsions with infinitesimal droplet deformation.The equation reads where τ ij is the deviatoric stress tensor, e ij is the rate-of-strain tensor, η s is the viscosity of the Newtonian ambient fluid, and the overdot (in _ τ ij and _ e ij ) denotes the partial time derivative.The parameters α 1 , β 1 , and β 2 are functions of the bubble volume fraction f and the relaxation time of a single bubble λ.The expressions of these parameters are presented in Table I.Solving Eq. ( 4) for a linear oscillatory shear flow, Llewellin et al. [4] obtained the following expressions for G 0 and G 00 : where G 00 red denotes the loss modulus without the viscous contribution of the solvent (i.e., ωη s ).Note that in the limit of large Cd, Eq. (5b) yields a negative value of G 00 for f .0:6; however, this is not a problem, because the model is valid only for dilute suspensions.Mitrias et al. [17] simulated the oscillatory shear flow of dilute monodisperse bubble suspensions to determine their viscoelastic moduli.They then compared their results with those from Eq. ( 5), finding good agreement for bubble volume fractions lower than 0.5%.
As mentioned above, Eq. ( 4) rigorously holds for dilute monodisperse bubble suspensions.To extend its validity to polydisperse bubble suspensions above the dilute regime, Llewellin et al. [4] modified the parameter β 1 as follows: where b is a fitting parameter.For b ¼ 9, the authors found that Eq. ( 4) described well their experimental data for polydisperse systems with bubble volume fractions up to 46%.However, their method was purely empirical and only applicable to their experimental system; thus, it cannot be considered conclusive in terms of determining the effect of polydispersity on the viscoelastic behavior of bubble suspensions.This was also stated by Mader et al. [6], who employed a more systematic, yet simple, approach to elucidate the role of polydispersity.In brief, they suggested that every dilute polydisperse suspension can be treated as the sum of N noninteracting monodisperse components, each with a characteristic radius R i and bubble volume fraction f i .Thus, to obtain the total material functions of the polydisperse suspension, one can calculate the respective quantities for each bubble size class and then sum the individual contributions.The authors provided a worked example for the steady-shear viscosity of a dilute polydisperse bubble suspension but did not cover linear viscoelastic material functions.Seo and Youn [18] suggested a phenomenological rheological model for dilute monodisperse bubble suspensions.Their constitutive equation is also expressed in the form of the linear Jeffreys model [i.e., Eq. ( 4)] and was based on the analysis of Maffettone and Minale [21] on the deformation of an ellipsoidal droplet in a simple shear flow.Similar to the Seo and Youn [18] 3 model of Doi and Ohta [22], the model of Seo and Youn [18] describes the macroscopic stress of the suspension, accounting for the evolution of the bubble size and shape.The parameters α 1 , β 1 , and β 2 for their model are shown in Table I.Note that the models of Llewellin et al. [4] and Seo and Youn [18] represent modified versions of the Jeffreys model [19], the difference in the coefficients α 1 , β 1 , and β 2 stemming from the analysis in the limit of a single droplet of Frankel and Acrivos [20] and Maffetone and Minale [21], respectively; therefore, these models do not account for fluid dynamic interactions or crowding effects among bubbles.Joh et al. [10] solved the constitutive equation proposed by Seo and Youn [18] for a linear oscillatory shear flow, obtaining the following expressions for G 0 and G 00 : To describe their polydisperse experimental data, they modified Eq. ( 7), based on the method of Mader et al. [6].But instead of considering a discrete bubble radii distribution, they generalized the linear superposition approach by treating the bubble radius R as a continuous variable.Therefore, in terms of a volume fraction density function, f (R), related to the bubble radii distribution, one can write Eq. ( 7) as follows: Above, the volume fraction density function is defined so that f (R)dR yields the volume fraction of bubbles with radius in the range dR around R. Comparing their polydisperse experimental data with the predictions of Eq. (8b), Joh et al.
[10] noted a deviation of G 0 at low values of ω.They attributed this to bubble fluid dynamic interactions, without however validating this assumption.As said, several studies have focused on the timedependent flow of suspensions with infinitesimally deformed bubbles (i.e., Ca ( 1), advancing constitutive equations to describe their viscoelastic behavior.However, there are sparse experimental data to validate these equations.In this work, we aim to provide a systematic experimental characterization of the linear viscoelastic behavior of semidilute bubble suspensions by employing a rheo-optical setup to visualize the behavior of the suspensions during SAOS rheological tests.This coupled approach offers a higher confidence in the experimental measurements, because it allows investigating the effects of bubble size distribution and various SAOS measurement artifacts, including bubble rise, coalescence, and changes in bubble spatial organization over time, which can influence the rheological measurements of bubble suspensions.

A. Chemicals
Bubble suspensions were generated using a Newtonian mixture of mineral oil (RTM37 Mineral Oil Rotational Viscometer Standard-produced by Paragon Scientific, Birkenhead, UK) and 5% w/w of the liquid, nonionic surfactant Span 80 (Sigma Aldrich, St. Louis, USA).The chosen ambient fluid is viscous enough to ensure time-stable suspensions and minimal bubble rise during the rheological measurements.The properties of the individual materials and the final mixture were measured at 12 C, in compliance with the rheological tests, and are presented in Table II.Additionally, in Sec. 1 of the supplementary material [23], we provide a detailed rheological characterization of the used ambient fluid.

B. Generation and rheological characterization of the bubble suspensions
To produce the bubble suspensions, we used an in-house manufactured aeration device consisting of a rotating propeller with attached aeration plates covered with 2 μm porous filters (Fig. 1).This system offered aeration and mixing simultaneously, ensuring uniform distribution of bubbles in the ambient fluid and generation of homogeneous suspensions.The aeration time depended on the desired bubble volume fraction, with longer aeration times resulting in denser suspensions.After producing the suspensions, a Silverson (L5 Series) high-shear mixer was used to reduce the size of the bubbles.Then, we determined the bubble  volume fraction gravimetrically, weighing 100 ml of the generated suspensions and using the following equation: To measure the density accurately, the suspensions were allowed to rest until their temperature reached that of the ambient fluid before aeration.Before the rheological measurements, images of the bubble suspensions were taken using a bright field microscope (Zeiss Axio Observer 5) and their bubble size distributions were obtained by image analysis using a MATLAB in-house code.
The rheological characterization of the generated bubble suspensions was carried out in an Anton Paar MCR302 stress-controlled rotational rheometer.The instrument was equipped with a Peltier plate for controlling the operating temperature and a parallel-plate geometry (d ¼ 40 mm) with sand-blasted surfaces (surface roughness of 2 μm) to prevent wall-slip effects.All tests were carried out with a gap of 1.3 mm at a constant temperature of 12 C.For the tested samples, the chosen gap was always ten times larger than the average bubble diameter, hence avoiding wall effects originating from the confinement between the plates.Moreover, the operating temperature allowed for minimal bubble rise during the rheological measurements.
To retrieve the viscoelastic moduli, we carried out SAOS rheological tests, applying a linear sinusoidal shear strain to all tested samples.To select an appropriate value for the strain amplitude, we performed strain sweep experiments on a bubble suspension with f ¼ 9:5%.We tested three oscillation frequencies, ω 1 ¼ 0:56 rad/s, ω 2 ¼ 1 rad/s, and ω 3 ¼ 10 rad/s, and we scanned a range of shear strains between 0.05% and 10%.We observed that both viscoelastic moduli remained constant with the applied strain, a result that suggests that any strain value within this range is appropriate for testing the linear viscoelastic properties of our system.Therefore, in the SAOS rheological measurements, we opted for a shear strain of 5%.The G 0 and G 00 curves obtained from the strain sweep experiments are reported in Sec. 2 of the supplementary material [23].
The SAOS measurements were performed by ramping down the oscillation frequency from 50 to 0.1 rad/s.We chose to operate from higher to lower frequencies, i.e., from shorter to longer characteristic flow times, to minimize the effect of bubble coalescence and rise on the obtained viscoelastic curves.To validate the accuracy of our experimental data, in terms of low torque and phase angle resolution issues, we performed a dedicated study, which is detailed in Sec. 3 of the supplementary material.Prior to the SAOS measurements, we presheared the bubble suspensions for 3 min at 0.1 s À1 to remove potential loading effects, without altering the bubble size distribution (as confirmed by the experiments discussed in Sec.III A).

C. Rheo-optical oscillatory shear experiments
The bubble suspensions were visualized under oscillatory shear, using the rheo-optical setup depicted in Fig. 1.To enable optical access to the tested samples, we replaced the bottom plate of the rheometer with a quartz glass plate (Anton Paar Peltier Universal Optical Device, P-PTD 200/ GL).We conducted the SAOS tests as previously described and captured images throughout the process using a FLIR GS3-U3-32S4M-C 1=1:8 00 camera (acquisition frequency: 5 Hz, image resolution: 1536 Â 2048 pixels), a Nikon mono zoom lens, and a white LED light as illumination.To prevent reflections from the metal, we applied black spray-paint to the sandblasted parallel plate geometry.The depth of field for the current visualization setup was 600 μm.To calibrate the focus plane, we fabricated a disk of 0.6 mm thickness by sticking together transparent laminated sheets of 0.3 mm thickness.Then, we attached a 1 mm calibration tape on top of the disk and focused on it.The pixel-to-mm conversion scale was found to be equal to 306 pixels/mm.All visualization experiments were conducted at a focus plane of 600 μm, with minor adjustments for different samples.Note that, although we focused on a precalibrated plane, owing to the large depth of field our visualization encompassed multiple bubble layers.This enabled us to visualize bubbles even when they had risen beyond the calibration plane, covering a visualization gap of up to 300 μm above the calibration plane.

A. Bubble size distributions
All the tested samples were polydisperse, with bubble radii following the gamma distribution between 20 and 120 μm (Fig. 2).For each bubble suspension, a volumeweighted average radius, hRi, can be determined using the De Brouckere average radius R 43 , which is the ratio of the fourth to the third moments of the bubble number density function, n(R); in terms of volume fraction density function, the same average radius is given by the ratio of the moments of order one and zero, hRi ; We chose to operate with a volume-weighted mean radius because it accounts for the bubble volume fraction, a property that is linked closely to the viscoelasticity of bubble suspensions.For the reported bubble volume fractions, f 1 ¼ 4:2%, f 2 ¼ 13:5%, and f 3 ¼ 19:2%, the respective volume-weighted mean radii were equal to 61, 64, and 65 μm.To ensure that the initial bubble size distribution did not change significantly during the preshear stage, we compared the bubble sizes computed from the microscope images to those derived from the images of the visualization experiments.For each sample, we used the image immediately after the end of preshearing, and we identified the bubbles and their radii (in pixels) using MATLAB.We then converted the pixels to mm, using the scale we determined through the focus plane calibration, i.e., 306 pixels/mm.Our image analysis showed that after 3 min of preshearing at 0.1 s À1 , the bubble radii distributions remained essentially the same.

B. SAOS results
To characterize the viscoelastic behavior of the bubble suspensions, we must compare the experimental values of G 0 and G 00 red to the predictions of an appropriate theoretical model.Since the aim is to use the simplest model that can describe the behavior correctly, we started with modeling the polydisperse suspension as a monodisperse one, using Eq. ( 4) with the parameters α 1 and β 2 expressed in terms of the volume-weighted average bubble radius hRi and of an average dynamic capillary number, hCdi ; hλiω, (11) where hλi denotes the average relaxation time of the bubbles, defined using hRi.
To proceed, it is necessary to determine which version of Eq. ( 4) is more appropriate: (i) the original Jeffreys version, (ii) the Llewellin et al. [4] modification, or (iii) the Seo and Youn [18] modification.By observing Eqs.(5b) and (7b), we see that for large values of hCdi, the second term in the expressions for G 00 red becomes dominant, resulting in negative G 00 red values.This is not physically reasonable.In particular, for our experimental system, substituting its physical parameters (i.e., solvent viscosity η s ¼ 53:063 Pa s, bubble volume fractions of the tested suspensions f 1 ¼ 4:2%, f 2 ¼ 13:5%, and f 3 ¼ 19:2%, and respective average relaxation times of the bubbles hλi 1 ¼ 0:88 s, hλi 2 ¼ 0:80 s, and hλi 3 ¼ 0:75 s) into Eqs.(5b) and (7b), we obtained negative G 00 red values for Cd larger than 0.8.This issue was not addressed in [4] and [10], despite their experimental parameters being similar to ours.We must note that both studies accounted for the viscous contribution of the solvent, which for semidilute suspensions is dominant.Thus, even if G 00 red was negative, it was not readily noticeable from the η 0 ; G 00 =ω and G 00 graphs presented in [4] and [10].Therefore, to describe the linear viscoelastic moduli of the suspensions investigated in this work, we opted for the original Jeffreys model [19] [i.e., Eq. ( 4) with α 1 ¼ β 2 ¼ hλi], because this model does not pose this problem.For an oscillatory shear flow, the original Jeffreys constitutive equation yields the following expressions for G 0 and G 00 red in a suspension of identical bubbles: For each of the reported bubble volume fractions, we determined the theoretical G 0 and G 00 red values using Eq. ( 12) expressed in terms of hCdi.In Fig. 3, we report the experimental G 0 , G 00 red curves (black and hollow points), along with the theoretical predictions of the Jeffreys model for a monodisperse suspension (red and black dashed lines), as functions of the average dynamic capillary number.As we see, there is good agreement between the experimental and theoretical values for hCdi larger than unity.But at lower hCdi values, the experimental values of G 0 are significantly larger than the theoretical ones.Moreover, the shape of the experimental G 0 curves is more complex, resembling that of a suspension with multiple relaxation modes.As shown in Fig. 4, the onset of the G 0 deviation is related to the bubble volume fraction, with denser suspensions deviating earlier, that is, at higher hCdi values.Furthermore, within the same low range of hCdi, the experimental values of G 0 deviate more from the theoretical predictions at higher bubble volume fraction, with the densest suspension deviating almost twice as much compared to the most dilute one.For a dedicated discussion on this, see Sec. 4 of the supplementary material [23].
This complex G 0 trend was unexpected, because it is not predicted by the model.This led us to believe that the model we considered may not be appropriate for describing the experimental data.A possible reason for the poor predictions at low hCdi values may be the polydisperse character of our suspensions.To account for this issue more thoroughly, we used the linear superposition approach suggested by Mader et al. [6].Considering a polydisperse bubble suspension with N discrete bubble radii classes, Eq. ( 12) can be generalized as follows: where Cd i ; λ i ω, with λ i and f i representing the bubble relaxation time and volume fraction of the ith bubble class, respectively.Using Eq. ( 13), we calculated the theoretical G 0 and G 00 red curves corresponding to the tested suspensions.These curves are also reported in Fig. 4 (black and blue solid lines).As seen, these new curves are almost identical to those derived by modeling the polydisperse suspension as monodisperse with a bubble radius equal to the volume-weighted average radius of the bubbles population.This is not unexpected, as one can prove that the theoretical predictions of the simpler Jeffreys model (i.e., Eq. ( 12) expressed in terms of hCdi ) differ from those of the more refined one (i.e., that based on the summation of the G 0 , G 00 red contributions of each bubble size class) only if the suspension bubble size distribution is bimodal, with equal volume fractions of very small and very large bubbles.In any other case, the two versions of the Jeffreys model yield nearly the same predictions.This is discussed in more detail in Sec. 5 of the supplementary material [23].
Because the generalized Jeffreys model, which accounts for polydispersity more accurately, yields almost the same results as its simpler version (that for monodisperse suspensions), we can conclude that polydispersity is not the root cause for the observed G 0 deviation at low hCdi values.To investigate this behavior further, we considered the influence of artifacts that may be induced by the rheological measurements.Specifically, due to their extended duration, SAOS tests may be susceptible to phenomena such as bubble rise, coalescence, and, more broadly, changes in the suspension microstructure over time.These phenomena are examined in detail in Secs.III C, IV A, and IV B.

C. Effect of bubble rise
To investigate the effect of bubble rise on the obtained G 0 curves, we first calculated the bubble rise velocity, using the Hadamard-Rybczynski equation [24,25], where hRi is the volume-weighted average bubble radius, ρ s and η s are the density and viscosity of the ambient fluid, respectively, and g is the gravitational acceleration.This equation is accurate in the limit of a single bubble rising in a clear liquid.The presence of other bubbles and surfactants in a suspension has been shown to retard the bubble rise [26,27], causing Eq. ( 14) to overestimate the rise velocity.We recognize that, while this equation provides the most conservative estimation of bubble rise during our experiments, its use also means that the computed rise velocities may be exaggerated, so that in our experiments the actual effect of bubble rise may not be as pronounced.In Table III, we present the calculated bubble rise velocity for all tested suspensions.
Based on the computed bubble rise velocity, as we see from Table III, at the end of the SAOS measurements that last about 35 min, bubbles should have risen appreciably, covering a distance equal to nearly 45% of the rheometer gap.Even if these values are overestimated, they still indicate that the effect of bubble rise may become important toward the end of the SAOS tests.To investigate whether bubble rise is pronounced from the beginning of the SAOS measurements, we calculated the distance covered by the rising bubbles during the initial 5 min of the experiment.This timeframe corresponds to the frequency range 50-10 rad/s, after which we observed the onset of the G 0 deviation.For all the tested suspensions, this distance was found to be negligible, indicating that bubble rise cannot be responsible for causing the observed G 0 deviation at low Cd h i values.To validate these calculations, we placed a sample of a suspension with f ¼ 10:4% on the rheometer plate and waited for 35 min, allowing the bubbles to rise freely before initiating the SAOS test.For the same suspension, we then performed a normal SAOS measurement without waiting time and we compared the obtained G 0 curves.As we see in Fig. 5, the two curves almost overlap, the complex G 0 trends being essentially the same regardless of bubble rise.This indicates that in our experiments, bubble rise is not the cause for the characteristic G 0 deviation observed at low Cd h i values.Nevertheless, it could potentially amplify this deviation toward the end of our experiments by affecting the mean interbubble distance.This observation was confirmed via time sweep experiments and SAOS measurements performed by increasing and decreasing the oscillation frequency.The results of these experiments are presented in Sec.6 of the supplementary material [23].
Specifically, Fig. XIIIa of the supplementary material shows that the G 0 deviation at low Cd h i values is present even in SAOS measurements performed by ramping up the oscillation frequency, where the influence from bubble rise is minimal.When the measurements are performed inversely, more time elapses before recording the low Cd h i viscoelastic data.During this time, bubbles tend to rise, leading to a decrease in the average interbubble distance, which amplifies the G 0 deviation.This observation suggests that there is a correlation between the deviation of G 0 and the average interbubble distance or, more broadly, the overall bubble spatial organization in our samples.
It is worth noting that in dense colloidal suspensions, the presence of clusters has been shown to increase the bulk elasticity of the suspension as the particles come into close proximity and create a network structure [28][29][30].This network effectively stores elasticity, leading to increased resistance to deformation and larger G 0 values.Even though the presence of bubble clusters or threads may not result in a network as organized and robust as that in colloidal suspensions, it certainly alters the average interbubble distance, influencing the microstructure and, in turn, the rheology of the suspension.From Fig. XIIIb of the supplementary material, we see that before the initiation of the SAOS measurements, bubbles do not appear to be uniformly distributed within the tested samples.This nonuniform bubble spatial distribution is most likely created upon loading our samples the rheometer plate and affects the average interbubble distance, leading to the characteristic G 0 deviation at low Cd h i values.We further investigate this aspect in Sec.V.

IV. VISUALIZATION OF BUBBLES UNDER OSCILLATORY SHEAR
Having examined the effect of bubble rise during our SAOS measurements, we then investigated whether the observed G 0 trends could be attributed to changes in the suspension microstructure over time.These changes can be related to the general bubble organization and the bubble size.To this end, we visualized the bubble suspensions under oscillatory shear and conducted a statistical image analysis.The results are discussed below.

A. Changes in bubble organization over time
To examine whether the pre-existing spatial organization of the bubbles right after the preshear changed significantly as a result of the SAOS measurements, we used the following dimensionless parameter, called coverage: C ; Total bubble surface (in pixels) Image surface (in pixels) : This parameter indicates the percentage of image covered with bubbles.When the bubbles align or come close and form clusters, they tend to overlap more; this leads to a decrease in the area of the image occupied by bubbles, with a consequential decrease in coverage.Hence, a reduction in coverage can be used as a qualitative metric for bubble clustering and alignment, or more broadly, for changes in the microstructure of the dispersed phase.To minimize the amount of data, we conducted the image analysis using 15 representative images for each oscillation frequency.We converted each image into a black and white matrix, with the bubbles portrayed as white pixels on a black background.
Then, we summed all the white pixels and divided them by the total number of pixels in the image, obtaining a characteristic value of coverage.For each frequency, we determined a mean coverage value C by averaging the coverage values of the corresponding 15 images.In Fig. 6(a), we present the mean coverage as a function of the average dynamic capillary number for the three reported bubble volume fractions.As we see, C remains constant for all tested suspensions, suggesting that the general bubble organization does not change significantly throughout the rheological measurements.

B. Bubble coalescence
To investigate potential bubble coalescence phenomena that would affect the suspension microstructure, we evaluated how the average bubble radius changed during the rheological tests.To this end, we analyzed the same images used to calculate the coverage to measure the bubble size distribution and the corresponding mean radius hRi during the oscillatory tests.For each oscillation frequency, we calculated an average value for hRi, denoted as hRi, obtained by averaging the hRi values of the 15 representative images.In Fig. 6(b), we show how hRi changed for the three reported bubble suspensions throughout the SAOS tests.As we see, the average bubble size slightly increases throughout the rheological experiments (these tests start at large hCdi values and then, in time, progress toward lower hCdi values), reaching a maximum increase of 7%-9% of the initial size.Therefore, based on the observed trend of coverage and the minimal increase in the average bubble size, we conclude that the suspension microstructure remains mostly stable during the rheological tests.
The failure of the Jeffreys model, and of other models (e.g., [18]) previously used for dilute and semidilute bubble suspensions, to predict accurately the G 0 trends at low hCdi values led us to investigate the impact of other nonlinear phenomena, specifically bubble fluid dynamic interactions, on the experimental data.As suggested by Joh et al. [10], such interactions may induce phenomena with longer relaxation times, resulting in the characteristic G 0 deviation at low hCdi values.Specifically, when bubbles are in close proximity or interact with each other, the fluid flow induced by the motion of one bubble can influence the motion and the behavior of the neighboring bubbles, leading to groups of bubbles behaving collectively.In this context, interactions are influenced by two key factors: (i) the bubble volume fraction and (ii) the mean interbubble distance.In an idealized system, where bubble can be considered perfectly distributed, these two parameters are directly related.But in our experimental system, technical aspects, such as rheometer loading and subsequent preshear protocol, can change the average interbubble distance for a fixed bubble volume fraction.To investigate this, we performed additional SAOS tests on the generated suspensions, using different preshearing conditions.Specifically, we presheared our samples for 33 min at 0.9 s À1 , then comparing the corresponding SAOS results with those obtained after 3 min of preshearing at 0.1 s À1 .The aim was to determine whether a change in the bubble spatial distribution, induced by longer preshear protocols, could mitigate the deviation from the theoretical predictions of G 0 .At this point, we must note that the second preshearing protocol was found to be the maximum preshear that could be applied without causing any significant changes to the bubble size distribution and the total bubble volume fraction of the suspensions.Our findings are discussed below.

V. BUBBLE FLUID DYNAMIC INTERACTIONS AND EFFECT OF PRESHEAR A. Rheological results
Following the method described in Sec.III A, we checked whether the initial bubble size distribution changed significantly during the stronger and more prolonged preshear stage.Our image analysis showed that when the samples were subjected to 33 min of preshearing at 0.9 s À1 , the bubble radii shifted to slightly higher values, while still following a gamma type distribution (Fig. XIV in the supplementary material).In Table IV, we present the updated volume-weighted average bubble radii for the reported suspensions.For each of them, we calculated the corresponding value of hCdi.
In Fig. 7, we compare the G 0 and G 00 red curves of the reported bubble suspensions, obtained after 3 min of preshearing at 0.1 s À1 and after 33 min of preshearing at 0.9 s À1 .
As we see, preshearing affects the resulting G 0 curves only for low bubble volume fraction.For f ¼ 4:2%, stronger and more prolonged preshearing leads to G 0 values closer to the theoretical ones, suggesting a more uniform redistribution of the bubbles, and therefore weaker fluid dynamic interactions among them.There is still a deviation between the experimental and theoretical values for average dynamic capillary numbers lower than 0.1, which we believe is due to bubble rise during the SAOS tests, which directly impacts the average interbubble distance and, in turn, the bubble interactions.As the bubble volume fraction increases, the experimental G 0 curves are almost insensitive to the applied preshearing conditions, suggesting that preshearing affects the mean interbubble distance negligibly when the suspensions are denser.The reported results indicate that there is no significant bubble rise during the longer preshearing stage.This is supported by both the rheological and optical findings, which do not show any significant differences in the high frequency plateau of G 0 and in the bubble size distributions obtained just after the preshearing (see Fig. XIV in the supplementary material), respectively.

B. Local spatial distribution of bubbles
To quantify the fluid dynamic interactions among bubbles due to their local spatial distribution, we performed a statistical image analysis, using the method described by Kudrolli et al. [31].In brief, for each of the reported bubble suspensions, we analyzed an image obtained after the end of the preshearing stage, considering both preshearing conditions.We divided each image into 64 equal squares, and we identified the number of bubbles in each square, using MATLAB.We then determined the distribution function for the local bubble number (n), i.e., the probability of having a certain number of bubbles in a cell.This method allows evaluating bubble interactions based on the assumption that a smaller number of bubbles in a given area indicates a larger average interbubble distance and thus weaker fluid dynamic bubble interactions.In Fig. 8, we present the probability distribution of the local bubble number for the tested bubble volume fractions and preshearing conditions.If the bubble suspensions were monodisperse and uniformly distributed, we would expect to have always the same number of bubbles in each cell.In this case, the probability function would be a narrow peak centered over this specific value of n.However, in our case, where the tested suspensions are polydisperse and the bubbles are not uniformly distributed, we obtain a wider distribution of local bubble numbers.Comparing between the two different preshearing conditions, we see that for the lowest bubble volume fraction, stronger and more prolonged preshearing results in narrower local bubble number distributions, which are also shifted toward lower values, meaning that the average interbubble distance decreases, resulting in weaker fluid dynamic interactions among the bubbles.However, as the bubble volume fraction increases, the two distributions almost overlap, confirming our previous rheological measurements and suggesting that in denser suspensions preshearing does not reduce the bubble interactions significantly.We highlight that for all volume fractions tested, the bubble size distribution, as opposed to the local bubble spatial distribution, is not altered by preshearing (as reported in Figure XIV of the supplementary material [23]), thus indicating that the shift observed in Fig. 8(a) can be solely associated with the more effective bubble redistribution.

VI. FITTING A MULTIMODE JEFFREYS MODEL
As the rheological measurements indicate, there is an unexpected increase in the suspension elastic modulus at low hCdi values.Based on the statistical image analysis, during the SAOS experiments the microstructure of the suspension remains mostly stable.Therefore, we believe that the G 0 increase is mainly due to bubble fluid dynamic interactions caused by the initial spatial distribution of the bubbles after the loading of the samples on the rheometer plate.As shown, these interactions are closely related to the bubble volume fraction and the applied preshearing conditions.When the suspension is subjected to high oscillation frequencies, the characteristic flow time is very small; thus, the measured viscoelastic moduli capture the response of the smallest length scale of the suspension microstructure, namely that of an individual bubble.This explains the good agreement between the real and the theoretical values at higher hCdi, where the generalized Jeffreys model predicts the viscoelasticity arising solely from the bubble interfaces accurately.As the oscillation frequency decreases, the characteristic flow time increases, enabling more complex relaxation phenomena associated with bubble fluid dynamic interactions.These interactions manifest as an increase in suspension elasticity, with the measured G 0 reflecting the response of bubbles that behave collectively owing to their proximity and mutual influence.
As previously discussed, the Jeffreys model, as well as the models of Llewellyn et al. [4] and Seo and Youn [18], does not account for bubble fluid dynamic interactions, thus failing to predict the G 0 deviation at low hCdi values.To address this limitation, we considered a model developed by Palierne [32], aimed at modeling the rheology of dilute and semidilute emulsions consisting of viscoelastic and Newtonian fluids.This model accounts for droplet fluid dynamic interactions by assuming that the local strain exerted on a single droplet is modified by the deformation of the surrounding droplets and that the interactions between droplets are of a dipole kind.For emulsions formed by two viscoelastic fluids, the model predicts a G 0 profile with a double shoulder.This profile typically exhibits a relaxation mode associated with the relaxation of the viscoelastic matrix at high oscillation frequencies and a secondary relaxation mode associated with droplet interface relaxation at low oscillation frequencies.While the G 0 profiles observed in our experiments also present multiple relaxation modes, we must highlight a fundamental difference: our suspensions consist of two Newtonian fluids, i.e., air as dispersed phase and a mixture of mineral oil and Span 80 as ambient fluid.These fluids relax instantly, unlike for emulsions of two viscoelastic fluids.Therefore, in our systems, the relaxation mode at high oscillation frequencies cannot be attributed to the relaxation of the ambient fluid.To confirm this, in Sec. 8 of the supplementary material [23], we compared the predictions of the Palierne and Jeffreys models with our experimental data.As we can see, the G 0 trend predicted by the Palierne model qualitatively resembles that of the Jeffreys model, both featuring only one characteristic relaxation time.This is not surprising given that, for emulsions of Newtonian fluids, the literature [33,34] widely acknowledges the similarity between the two models.Similar to the Jeffreys model, the Palierne model fails to predict the second G 0 shoulder at low hCdi values.In addition, it does not capture the high hCdi plateau as accurately as the Jeffreys model.These shortcomings likely stem from the assumptions of the Palierne model about the type of droplet interactions.Specifically, the model assumes that the interactions are of a fluid dynamic dipole nature and that, within the interaction range, droplets are uniformly dispersed.Palierne [32], along with subsequent relevant studies, [35,36] clearly acknowledges that the model fails to accurately predict experimental results for systems where the aforementioned assumptions are not met.
The interactions that we believe cause the observed G 0 deviation in our systems differ from those considered by the Palierne model.In our case, bubbles are locally distributed nonuniformly, as evidenced by our rheo-optical experiments.Consequently, in some regions of the sample, bubbles are closer than they would be if they were uniformly distributed, making their local interactions stronger and causing them to behave and relax collectively.Similar to our case, Bousmina and Muller [35] as well as Carreau et al. [36] report the failure of the Palierne model to predict the low frequency G 0 shoulder, which they also attribute to interactions caused by the presence of aggregates in their systems.Despite the absence of a model that well describes our experimental data at low hCdi values, our rheo-optical experiments indicate a nonuniform spatial distribution of bubbles within the ambient fluid.Drawing from existing literature and having excluded experimental artifacts related to the performed rheological measurements, we infer that these higher-order interaction phenomena lead to the longer relaxation times observed in our experiments.
To validate this argument and quantify the complex relaxation phenomena arising from the bubble interactions, we fitted a multimode Jeffreys model to the experimentally determined values of G 0 and G 00 red .Then, we examined how the bubble volume fraction and the applied preshearing conditions affected the suspension relaxation modes.
We performed the fitting employing the Curve Fitting App of MATLAB.Initially, we fitted Eq. (13a) to the experimental G 0 values, using the relaxation times λ i and the volume fractions f i as fitting parameters.In all cases, we started from a single-mode Jeffreys model and progressively added modes, until we achieved the best fit.The maximum number of relaxation modes was determined through a simple convergence analysis, assessing whether the addition of an extra relaxation mode resulted in a noticeable reduction is the sum of squared errors (SSE).We then substituted the determined λ i and f i values in Eq. (13b), to obtain the corresponding G 00 red curve.In Table V, we present the values of the fitting parameters for the various bubble volume fractions and preshearing conditions, along with the corresponding R 2 value for each fit.
As can be observed, the smallest relaxation time was very similar in all cases, closely matching the relaxation time of an individual bubble calculated using the Jeffreys model and the volume-weighted mean bubble radius of each suspension.This finding supports our argument that at higher oscillation frequencies and, in turn, shorter characteristic flow times, we effectively capture the relaxation of a single bubble.The number of relaxation modes increases with the bubble volume fraction, suggesting that bubble interactions, and consequently complex relaxation phenomena, are more pronounced in denser suspensions.This complements our rheological measurements, which revealed that the onset of the G 0 deviation starts at higher hCdi values for larger bubble volume fractions.
We should mention that the relaxation time derived from the Jeffreys model, i.e., λ as defined in Eq. ( 1), refers to the case of a solitary bubble relaxing independently.In such a scenario, the bubble relaxes to its original undeformed state without the influence of neighboring bubbles.In the case of emulsions, it has been shown that the presence of neighboring droplets can influence the shape relaxation process of a single droplet.To further investigate the effect of bubble shape relaxation on the obtained G 0 trends, we calculated the bubble shape relaxation time in our suspensions using the Palierne expression [34], which considers the effect of f on the shape relaxation of the single droplet.This expression is a modification of the relaxation time derived from the Jeffreys model, incorporating a function of the bubble volume fraction.The equation and the calculated bubble shape relaxation times are reported in Sec. 9 of the supplementary material [23].The bubble shape relaxation times were found to be close to the relaxation times derived from the Jeffreys model.In fact, comparing the shape relaxation times with the fitting results reported in Table V, we observe that they fall between the first and second relaxation times computed through fitting for all tested suspensions.However, our fitting results revealed the presence of additional relaxation times which are at least an order of magnitude larger compared to the computed bubble shape relaxation times.This suggests that even though the relaxation mode associated with the shape relaxation of the bubbles can be influenced by crowding effects, it cannot be considered responsible for the characteristic G 0 deviation at lower hCdi values, which is clearly associated with longer relaxation times.
It is important to note that the sum of the f i values obtained from fitting does not always match the measured total bubble volume fraction of the suspension.This is reasonable considering that, to study the intricate behavior of the bubble suspensions, we resorted to an idealized bubble suspension.In Fig. 9, we report a schematic of this idealization.As shown in the schematic, the tested suspensions are not uniformly distributed when the SAOS measurements begin.Due to their spatial distribution, some bubbles are in close proximity, behaving collectively as one group.When fitting the predictions of the model to the experimental results, we assumed that each group of bubbles behaves like a single bubble with a radius equivalent to an effective (interaction) radius.In this context, we do expect that the sum of the f i values computed via fitting may exceed the actual bubble volume fraction of the tested suspension-especially in denser suspensions, where interactions become more pronounced.
Nevertheless, in the case of the most dilute suspension (f 1 ¼4.2%) subjected to stronger and more prolonged preshearing, we observed a close alignment between P f i and the bubble volume fraction obtained experimentally.This suggests that the applied preshear effectively redistributed the bubbles, causing the range of bubble fluid dynamic interactions to coincide with the volume-weighted radius of a single bubble.In this case, the interactions among bubbles were negligible, causing most bubbles to behave individually, without being influenced by the neighboring bubbles.This is confirmed by the decrease in the number of relaxation modes and the increase in the relative bubble volume fraction associated with the shortest relaxation time.However, as the overall bubble volume fraction increased, the applied preshearing conditions did not alter the number and the characteristics of the relaxation modes significantly.This aligns with our previous rheological measurements and image analysis results, suggesting that in denser suspensions preshearing impacts negligibly the bubble fluid dynamic interactions and the associated complex relaxation phenomena.
It is worth mentioning that, in the current study, the effect of bubble fluid dynamic interactions is noticeable due to the applied flow field.In the steady-shear experiments presented in our previous article [8], we successfully recovered the 3 min preshear at 0.1 s −1 33 min preshear at 0.9 s −1  9. Schematic depicting the range of bubble fluid dynamic interactions, as emulated when fitting our experimental G 0 , G 00 red data using a multimode Jeffreys model.
Taylor equation [37] for the zero-shear viscosity, a result that indicates that bubble interactions were not appreciable.However, in SAOS experiments, designed to probe the suspension microstructure, interactions, even if mild, become relevant, because they affect the relaxation process of the suspension and, consequently, the observed viscoelastic trends.Additionally, in our previous work [8], we explored bubble volume fractions up to 10.4%, while in this study, we extended the range to 19.2%.Although both studies involve semidilute suspensions, this work involves a significantly higher bubble volume fraction, which is expected to cause more pronounced bubble fluid dynamic interactions.

VII. CONCLUSIONS
In this work, we characterized the linear viscoelastic behavior of semidilute polydisperse bubble suspensions with a Newtonian ambient fluid.To determine the suspensions viscoelastic moduli, G 0 and G 00 red , we performed SAOS rheological tests with a preshear stage of 3 min at 0.1 s À1 .Comparing our experimental G 0 curves with the theoretical predictions of the original Jeffreys model, we found a good agreement for hCdi values larger than unity.But for lower hCdi values, the measured G 0 was larger than expected, this deviation occurring earlier in more concentrated suspensions.To elucidate this behavior, we systematically explored various aspects that could cause this behavior, including polydispersity, bubble rise, coalescence, and changes in suspension microstructure over time.
We showed that for the suspensions investigated in this work, the predictions of the generalized Jeffreys model accounting for polydispersity and of the simple Jeffreys model for monodisperse suspensions (used for bubbles with radius equal to the volume-weighted mean bubble radius of the suspension) yield essentially the same results.This indicates that polydispersity in itself is not the reason for the observed deviation between the experimental results and the model predictions.
As the observed trends could not be attributed to polydispersity, we investigated the effects of microstructure, examining whether it changed significantly over time due to artifacts related to the performed rheological measurements.To this end, we visualized the generated bubble suspensions under linear oscillatory shear and performed a statistical image analysis, examining the effects of bubble rise, coalescence, and spatial organization.Our findings indicated that, in general, the suspension microstructure was preserved during the SAOS measurements, suggesting that the observed viscoelastic trends were not the result of experimental artifacts.The failure of the Jeffreys model to accurately predict the G 0 trends at low hCdi values led us to investigate the impact of bubble fluid dynamic interactions on the experimental data.We believe that these interactions are induced by the spatial arrangement of the bubbles on the plate of the rheometer, present before the initiation of the SAOS measurements, and lead to complex relaxation phenomena, which become evident at longer characteristic flow times.To see whether a variation in the bubble spatial distribution, induced by longer preshearing protocols, could mitigate the deviation from the theoretical predictions of G 0 , we performed additional SAOS experiments with a preshear stage of 33 min at 0.9 s À1 .Our results showed that, for dilute bubble suspensions, stronger and more prolonged preshearing led to G 0 values closer to the theoretical predictions.But as the bubble volume fraction increased, the applied preshearing conditions had no significant impact on the experimental G 0 values.
To validate this, we fitted a multimode Jeffreys model to the viscoelastic moduli obtained experimentally.In line with the rheological measurements and the image analysis results, we showed that bubble interactions cause a complex relaxation process, consisting of multiple relaxation modes.The number of relaxation modes increased with the bubble volume fraction, indicating that the effect of bubble interactions amplifies in denser suspensions.Finally, all our results suggest that stronger and prolonged preshearing can effectively reduce the fluid dynamic interactions among bubbles when the bubble volume fraction is low.But as the bubble volume fraction increases and the bubble fluid dynamic interactions start dominating, preshearing does not impact the local bubble spatial distribution significantly, leading to the same G 0 trends independently of the applied preshearing conditions.

TABLE III .
Bubble rise velocity for different bubble volume fractions.
FIG.5.Effect of bubble rise on G 0 and G 00 red of a polydisperse bubble suspension with f ¼ 10:4%.

TABLE IV .
Volume-weighted average bubble radius ⟨R⟩ for different preshearing conditions.

TABLE V .
Fitting parameters for different bubble volume fractions and preshearing conditions.