Delayed elastic contributions to the viscoelastic response of foams

We show that the slow viscoelastic response of a foam is that of a power-law fluid with a terminal relaxation. Investigations of the foam mechanics in creep and recovery tests reveal that the power-law contribution is fully reversible, indicative of a delayed elastic response. We demonstrate how this contribution fully accounts for the non-Maxwellian features observed in all tests, probing the linear mechanical response function. The associated power-law spectrum is consistent with soft glassy rheology of systems with mechanical noise temperatures just above the glass transition [Fielding et al., J.Rheol. 44, 323(2000)] and originates from a combination of superdiffusive bubble dynamics and stress diffusion, as recently evidenced in simulations of coarsening foam [Hwang et al., Nat.Mater. 15, 1031 (2016)].


I. INTRODUCTION
Liquid foams are dense packings of bubbles dispersed in a surfactant-rich continuous phase. 1 At large enough packing fractions, the foam elasticity is determined by the resistance of bubbles to be deformed, which is proportional to the surface tension and inversely proportional to the bubble radius. 1,2 However, unlike other over-jammed systems, foams are not only elastic. Indeed, due to differences in Laplace pressure between small and large bubbles, foams coarsen in time. 1,3 This process continuously perturbs the stress balance between contacting bubbles, eventually resulting in configurations that exceed the local yield conditions, which in turn triggers local bubble rearrangements. [4][5][6] This structural relaxation mechanism provides the means of relaxing macroscopically imposed stresses 4,7 such that foams typically display a very slow fluid-like response. 2,7,8 Assuming that bubble rearrangement events occur stochastically, 6 one would a priori expect stress relaxation to be exponential, a characteristic of "Maxwell" fluids.
However, the slow viscoelastic response of foams strongly differs from that of simple "Maxwell" fluids that are solely characterized by an instantaneous elastic response and a single relaxation time. 9 Indeed, while the foam responds to the application of a small constant stress by an instantaneous elastic deformation, an intermediate creep phase precedes the terminal relaxation. 2,8 Similarly, the frequency-dependent response to an oscillatory strain is characterized by a storage modulus that gradually decays with decreasing frequencies, before exhibiting the hallmarks of viscous relaxation. 2,7,10,11 In this work, we show that these non-Maxwellian features relate to a time-dependent creep strain that is reversible. This delayed elastic contribution displays a power-law behavior, which is directly probed in strain recovery experiments following the release of a constant applied stress. We find that accounting for delayed elastic contributions in addition to Maxwellian behavior fully describes the linear mechanical response functions of foam at long times and equivalently at low frequencies. Our work thus identifies foam as a power fluid with a final relaxation, reminiscent of systems near the jamming transition. Within the class of over-jammed systems, this feature denotes foams as unique: they are over-jammed, yet do relax due to coarsening-induced structural relaxation events, which govern both delayed elastic storage and dissipation.

II. EXPERIMENTAL
Our experimental system is a Gillette shaving cream with coarsening characteristics similar to those of Gillette foams used in numerous previous studies. [4][5][6][7][8]10 It comprises propellant gases isobutane and propane, dispersed at a volume fraction of ϕ = 91.5% ± 0.3% in an aqueous mixture of long-chain fatty acids and sodium lauryl sulfate. All rheological tests are performed at 25 ○ C starting at tw = 5000 s, with tw defining the time elapsed since the foam production. By that time, the average bubble radius is 52 ± 5 μm, and the foam has evolved to the dynamical scaling regime of coarsening, 1,3,11 as shown in Fig. S1 of the supplementary material.
Rheological experiments are performed by using a stresscontrolled rheometer from Anton Paar (MCR series), equipped with a wide-gap Couette geometry (gap size 5.02 mm), which ensures that the mechanics probed is that of a bulk foam. Both cup and cylinder of the Couette cell are sanded to avoid wall slip, and the cell is sealed using a solvent trap to minimize drying and excessive gas exchange with the environment. The foam is loaded right after production through a hole at the bottom of the cup with the cylinder locked in measuring position, as described in Ref. 7. To minimize drift effects due to the air-bearing connection of the cylinder to the motor, we set the measuring position to the angle at which the residual stress observed after releasing a constant stress applied on a viscous fluid (glycerol) is minimal (<10 −4 Pa). 12 Such precaution becomes particularly important in recovery experiments, where residual torques may affect the signals when the shear rates approach the resolution limit. 12 The following rheological tests are performed within the linear range of viscoelasticity as demonstrated in Fig. S4 of the supplementary material: (i) stress relaxation tests measuring the temporal evolution of the stress σ(t) after application of a step strain of γ 0 = 0.1%, (ii) oscillatory strain experiments measuring the frequency dependence of the storage G ′ (ω) and loss modulus G ′′ (ω) with γ 0 = 0.1%, and (iii) creep and recovery tests measuring the temporal evolution of the strain γ(t) after applying a step stress of σ 0 = 0.2 Pa for a time of tc = 1000 s, and subsequently resetting the stress to zero. This procedure guarantees quasi-stationary creep conditions as the bubble size does not change by more than 8% during stress application, as shown in Fig. S1(b) of the supplementary material. Since the oscillatory stain experiments typically last 5000 s, we correct the values of the moduli by taking into account the drop of the foam elasticity during the test, as detailed in Figs. S2 and S3 of the supplementary material. All experiments are repeated at least six times following the loading procedure described above, where we report the averaged data in the figures hereafter and use the standard deviation as a measure of error.

III. RESULTS AND DISCUSSION
The results of all three tests clearly indicate that the foam is not simply characterized by an elastic modulus and a single relaxation process, as this is the case for Maxwell fluids. 9 Instead, stresses relax following a more complex pathway. This can be appreciated in Fig. 1(a), where we report the time dependence of the relaxation modulus G(t) = σ(t)/γ 0 obtained in stress relaxation tests. Estimating the foam viscosity as ηR = ∫ +∞ 0 G(t)dt ≈ 5 × 10 4 Pa. s and the short-time or equivalently high-frequency modulus as G∞ ≈ 190 Pa, the expected Maxwellian response would be G(t) = G∞e −t/τ R , with a relaxation time of τR = η R /G∞ ≈ 260 s. 9 This description, depicted as a dashed line in Fig. 1(a), is, however, clearly inconsistent with our experimental findings. Deviations from Maxwellian characteristics are also apparent in the frequency dependence of G ′ (ω) and G ′′ (ω), as shown in Fig. 1(b). Indeed, while the mechanical response of a Maxwell fluid is characterized by a sharp crossover between the high-frequency regime, where G ′ (ω) is independent of ω and G ′′ (ω) decreases with increasing frequency, and the low-frequency range where both moduli increases, the foam displays a markedly different behavior. In the high frequency range, the loss modulus increases as a square root of frequency, a feature that has been attributed to nonaffine motion along random slip planes. 7,13 Moreover, deviations from the Maxwell behavior are also observed in the lower frequency range. Indeed, the storage modulus gradually decays upon decreasing the frequency, instead of being almost frequency independent before exhibiting the hallmarks of relaxation. While these observations are consistent with results of previous investigations, 2,7,8,10,11 the description of the low-frequency characteristics of G ′ (ω) and G ′′ (ω), or equivalently the long-time behavior of G(t), remains unclear.
New insight is gained by examining the results of the creep and recovery tests. Figure 2(a) displays the time dependence of the strain during creep and recovery, normalized by the stress σ 0 applied during the creep test, J(t) = γ(t)/σ 0 . Formally, J(t) is the compliance for data acquired in creep tests, while for recovery, the normalization with σ 0 is a convenient way to compare creep and recovery responses within the same framework. For the creep experiment performed by applying a constant stress, the foam elasticity appears as an instantaneous rise of J = 1/G∞, with oscillations due to the tool inertia. 14 This is followed by an increase of J(t) that reflects the foam relaxation processes. As expected, a compliance describing the behavior of Maxwell fluids J(t) = 1/G∞ + t/η R fails to describe the data. 9 Reporting the derivative of J(t) as a function of time reveals that this discrepancy corresponds to a gradual decrease of the compliance rate, which reaches the constant value 1/η R only at long times, as shown in Fig. 2 The nature of this decrease is uncovered upon releasing the stress in recovery tests. Indeed, beyond the elastic recoil of 1/G∞ expected for a Maxwell fluid, we find that J(t) further decreases in time, as shown in Fig. 2(a). This delayed recovery is characterized by a power-law recovery rate, which matches the decay of the creep rate, as shown in Fig. 2(b). This suggests that the initial creep is also best described by a power law and more importantly that it is due to a time-dependent, yet recoverable contribution, thereby indicating a delayed elastic contribution to creep. [15][16][17] To describe our data, we thus write the creep compliance as a sum of three terms where the delayed elastic contribution is described by the Jeffreys-Lomnitz law of amplitude m and exponent 0 < α < 1. 18,19 Note that for t ≫ τ 0 , the Jeffreys-Lomnitz law reduces to a simple power law of J(t) ∼ t α . As shown in Figs. 2(a) and 2(b), Eq. (1) fits our data remarkably well, yielding η R = 4.5(±0.6) × 10 4 Pa. s and G∞ = 192 ± 7 Pa, in good agreement with the values inferred from the analysis of the stress relaxation test. For the parameters of the central term, we find τ 0 = 0.3 ± 0.1 s, m = 0.09 ± 0.01, and a rather small creep exponent α = 0.17 ± 0.03. For recovery times small compared to the creep test duration tc = 1000 s, we expect that recovery is well approximated by a simple reversal of the elastic contributions, with J(tc) being the compliance reached at the end of the creep experiment and t being the experimental timescale of the recovery test. 9 The good agreement between this description and the experimental data shown in Figs. 2(a) and 2(b) corroborates that short-time creep is caused by a delayed elastic storage process, which is reversible. The power-law behavior suggests that this contribution is due to a process with a broad distribution of timescales, which appears naturally when describing the central term of Eq. (1) as a series of Kelvin-Voigt elements, where ρ(τ) denotes the retardation spectrum. 9,19,20 This spectrum can be obtained as  Note that the peak in the relaxation spectrum is at ≈500 s, considerably above τ R . The red-dotted line is the relaxation spectrum expected for τ R → ∞. (b) Dynamic susceptibility χ ′′ (ω) = G ′′ (ω)/G∞ as a function of frequency. Symbols denote the experimental data corresponding to those shown in Fig. 1(b). Red-solid line shows the behavior expected from the description of the creep data. Maxwell relaxation (dashed-black line) leads to the appearance of a low-frequency relaxation peak on top of the delayedelastic "wing" corresponding to τ R → ∞ (red-dotted line).
To fully validate our modeling, we examine whether the stress relaxation of our foam can be described within the same framework. Since ρ ≠ 0, we need to assume that stress relaxation is governed by a distribution of relaxation times, 9,17,20,21 where the relaxation spectrum μ relates to ρ as 9,20,21 taking the integral over u as a principal value. 21 With ρ, m, and τR = η R /G∞ known from the fit of the creep data, we compute μ reported in Fig. 3(a) and then calculate G(t) using Eq. (5). The resulting G(t) describes the experimental data remarkably well, as shown in Fig. 1(a). Equivalently, we can determine G ′ (ω) and G ′′ (ω) as a function of frequency by again summing up the Maxwell responses over the relaxation spectrum μ. 9 As shown in Fig. 1(b), the result accounts for the characteristics of the experimental data in the low frequency range of interest. Note that the result is also more accurate than that obtained by using the usual Fourierbased method to transform G(t) into G * (ω), 9 as shown in Fig. S4(d) of the supplementary material. The fact that we can use the sole description of the creep data by Eq. (1) to accurately recast the results obtained in stress relaxation and oscillatory shear tests unambiguously shows that the process leading to time-dependent reversible strain in recovery is at the origin of the slow non-Maxwellian relaxation of foam.
Compared to other conversion methods, 2,7,22 the approach used here has the advantage of being analytically solvable, providing details about the distribution of relaxation times μ. As denoted in Eq. (6), this distribution solely depends on ρ, the parameter m, and the terminal relaxation time τR = η R /G∞. In the limit τR → ∞, μ decays as μ(τ) ∼ 1/τ 1+α , as shown in Fig. 3(a). This decay is consistent with the predictions of soft glassy rheology (SGR), assuming a mechanical noise temperature of x = 1 + α ≈ 1.17 just above the glass transition. [23][24][25] The corresponding relaxation modulus and low-frequency response are shown in Figs. 1(a) and 1(b), respectively. For our foam, however, τR is finite, resulting in a cutoff of the power law and the appearance of a peak in the spectrum at ≈500 s. This timescale is considerably above the Maxwellian relaxation time obtained from the creep data (τR = 230 s) and is consistent with the inverse crossover frequency inferred from the extrapolation to G ′ (ω) = G ′′ (ω) of the data shown in Fig. 2(b). Such a shift of the effective relaxation to lower frequencies is caused by the ρ-dependent integral in Eq. (6) and reflects that the delayed elastic contributions effectively delay terminal relaxation.
Note that the use of the Jeffreys-Lomnitz law for the central term in Eq. (1) is also motivated by the possibility of interpolating between the behavior of Maxwell fluids with α = 1 and the logarithmic creep ∼ ln(1 + t/τ 0 ) with α → 0, 18,19 which is generally observed in glassy systems. 16,19,26 Conveniently, these two cases are also consistent with the predictions of SGR for the liquid phase (x > 2) and the glass phase (x < 1), respectively. 25 The Jeffreys-Lomnitz law, therefore, provides an empirical formula that allows us to obtain retardation spectra analytically for a range of soft materials where power-law or logarithmic creep is commonly observed. 16,[26][27][28][29]

IV. CONCLUSION
In conclusion, our experiments show that the rheology of foam under quasi-stationary conditions retains features of power-law fluids but eventually exhibits a final relaxation. The experimental approach adopted highlights the advantage of using creep and recovery tests where the power-law contribution is assessed unambiguously, in contrast to other rheological tests where the different elastic and dissipative contributions are not simply additives. Most importantly, the recovery tests show that power-law creep is related of Chemical Physics ARTICLE scitation.org/journal/jcp to a process that allows for delayed elastic recovery. [15][16][17] To account for this behavior, let us recall that the coarsening process in foam leads to intermittent restructuring events. [4][5][6] In creep experiments, a restructuring event leading to the relaxation of local stresses entails that the macroscopic stress will be carried by the rest of the system. 16 This results in an increase in strain that can only be recovered in time; this is because upon releasing the stress, a full recovery of the stress-bearing part of the system will be hindered by the zones that were previously relaxed and now resist the deformation imposed by the instantaneous elastic recoil. Further restructuring events will then be needed to allow for full recovery. 16 Since the restructuring by intermittent events is a time-dependent process, full elastic recovery is delayed. How intermittent dynamics in coarsening foams lead to powerlaw rheology has been recently established in simulations, 30 where both bubble dynamics and stress fluctuations were monitored. In agreement with recent experiments, 31 the mean-square displacement of the bubbles was found to be superdiffusive, while the stress fluctuations were diffusive. Using microrheology arguments, 32,33 the authors show that these two features combined lead to G * (ω) ∼ ω α with α = 0.17, 30 in good agreement with our experimental findings. However, the simulation results do not account for the contribution of the final relaxation observed experimentally and instead resemble the response for τR → ∞ shown in Fig. 1(b). A similar response is also predicted in SGR for fully equilibrated systems with a mechanical noise temperature just above the glass transition. [23][24][25] However, for partially equilibrated systems, SGR predicts a creep behavior similar to that observed experimentally for timescales below the age of the system. 25 Beyond this limit, the transition between equilibrated and nonequilibrated creep modes leads to an age-dependent transient relaxation, followed by a long time power-law creep. 25 For our foam, we find such a transition when performing long experiments during which the bubble size significantly increases such that the quasi-stationary conditions realized in the experiments presented here no longer apply; this work will be presented in a future publication.
Finally, let us note that the frequency dependence of the loss modulus of our foams bears resemblances with the "excess wings" in the dielectric loss of deeply supercooled molecular systems. 34 Indeed, the mechanical analog χ ′′ (ω) = G ′′ (ω)/G∞ shown in Fig. 3(b) displays a low-frequency relaxation peak with a "wing" at higher frequencies. For supercooled liquids, recent simulations attribute the wing to rare and localized rearrangement events, while the peak relates to the growth by "dynamical facilitation" of regions in the material that have already rearranged. 35,36 Within the framework of foam mechanics, the wing can be considered corresponding to delayed elastic storage, consistent with the idea that the few intermittent events occurring within a short lag time will lead to an increase of strain in a matrix that has not yet been restructured-a contribution that can be fully recovered. Whether the final relaxation of foam can be attributed to a growth of already rearranged zones or to a percolation of the rearranged zones has yet to be established.
Let us close by highlighting the somewhat unique features of foams. They are over-jammed systems that are to be considered athermal due to the mesoscopic size of the bubbles. However, unlike other over-jammed systems, such as compressed emulsions, 37 they structurally relax due to dynamics driven by coarsening. These features confer foams the properties of systems that are near the jamming transition, where mechanical noise effectively leads to local structural relaxation. However, while the properties of jammed materials sensitively depend on shear history, 38,39 the coarsening process in foam naturally erases shear history in time by constantly re-injecting mechanical noise into the system. This makes foams ideal systems to study the mechanics of almost jammed systems in the linear range. 40

SUPPLEMENTARY MATERIAL
See the supplementary material for details on foam evolution by coarsening and the linearity of all rheological tests performed.

ACKNOWLEDGMENTS
We thank Emanuela Del Gado, Jan Vermant, and John Crocker for very helpful discussions. Financial support from the Swiss National Science Foundation (Grant No. 200021_172514) is gratefully acknowledged.

Conflict of Interest
The authors have no conflicts to disclose.

DATA AVAILABILITY
The data that support the findings of this study are openly available in Zenodo at http://doi.org/10.5281/zenodo.6381567.

Evolution of bubble characteristics and its impact on foam elasticity
To be able to assess the impact of foam evolution on our rheological tests, we characterize both the bubble size distribution and the elastic modulus as a function of the foam age .
The age-dependent bubble size distribution ( ) is characterised using the method described by Gaillard et al. [1], with results shown in Fig. S1 (a). For coarsening foam, we generally expect the system to reach a stage where such distribution is self-similar, while the mean bubble size � keeps increasing [2,3]. For our foam, we find that for ≥ 3000 s, ( ) indeed collapses onto a unique master-curve when renormalizing by the mean bubble size � , as shown in the inset of Fig. S1(a).
As shown in Fig. S1(b), we find that the age dependence of the mean bubble size is well described by with � (0) ≈ 20 μm, 0 ≈ 820 s and ≈ 0.49. This result is consistent with �~1 /2 expected for bubble coarsening of dry foams in the scaling regime [2,3], and is in good agreement with findings of previous investigations on Gillette foams [4][5][6][7][8].  For our experiments starting at = 5000 s, these findings show (i) that the foam is well within the scaling regime, and (ii) that the mean bubble size increases from � ≈ 52 μm to � ≈ 56 μm over the duration of a creep experiment ( = 1000 s), as denoted in Fig. S1(b). This relatively small increase of ≈ 8% ensures quasi-stationary conditions as mentioned in the main text.
As the elastic modulus of a foam typically scales as the inverse of the bubble radius [2,9], we expect the foam elasticity to decrease as a function of the foam age. To assess this dependence, we measure the elastic modulus of the foam at different ages in two different ways, using creep experiments on the one hand, and oscillatory strain experiments on the other hand.
For the creep experiments, we exploit the creep ringing effect observed due to tool inertia. As an example, we show the data obtained for = 5000 s in Fig. S2(a). Following reference [10], we fit the short-time oscillations of the compliance to where , , and are ringing parameters, and ∞ is a measure of the high frequency modulus. We repeat this operation at different foam ages and report ∞ ( ) in Fig. S2(c). As expected from the evolution of the bubble size, ∞ ( ) decreases with and is well described by where is an apparent plateau modulus and a characteristic frequency. We then retrieve the values of ′ and ′′, as shown for = 5000 s in Fig. S2(b). Since the acquisition of high frequency data (100 -10 rad/s) is fast, the foam elasticity barely changes during the acquisition performed by decreasing the frequency from 100 rad/s downwards. The foam age is thus unambiguously defined by at the start of the experiment. As shown in Fig. S2 (c), ( ) is somewhat smaller than ∞ ( ), but exhibits the same age-dependence. Indeed, a fit to the form of Eq. S3 yields 0 ≈ 1000 , and ≈ 0.49 in close agreement with the fit parameters of ∞ ( ), while (0) ≈ 408 Pa is somewhat smaller than ∞ (0).
As already mentioned for the evolution of the bubble size, these findings show that the modulus does not change by more than ≈ 8% during our creep experiments ( = 1000 s), such that we can consider the foam to be in a quasi-stationary state during the test.
However, the acquisition of the frequency dependent moduli down to = 5 × 10 −3 rad/s shown in Fig. S2(b) requires a measuring time of 5000 s. During this time laps the foam elasticity decays by ≈ 25%. As this is non-negligible, the resulting complex modulus * ( , + ) depends on the time passed since the start of the experiment at . To correct for this, we exploit the known dependence of ( ) on the foam age to define a stationary modulus as * ( , ) = * ( , + ) ( )/ ( + ), from which we obtain ′ and ′′ reported in Fig. 1(b) of the paper. Note that performing the correction with ∞ instead of leads to the same result, as both moduli exhibit the same age dependence. To test the robustness of this correction we perform an experiment in which the frequency is increased from low to high frequencies, which is the inverse of our standard protocol. The raw data obtained for both type of experiments, each started at = 5000 s, display significant differences, as shown in Fig. S3(a). After correction, both data sets are in good agreement, as shown in Fig. S3(b). This validates our correction method, which permits to retrieve the frequency behaviour of the foam in quasistationary conditions.
Let us finally remark that we accounted for creep ringing in our analysis of the creep compliance described in the paper. We however do not include the ringing term in Eq. 1 in order to focus on the characteristics of the foam mechanics, rather than on features set by the tool inertia.  Fourier transform of the stress relaxation data (purple symbols). The good agreement at lower frequency shows that linearity is fulfilled. The transform fails at high frequency due to the limited time-resolution of the stress relaxation test. The description obtained from creep and discussed in the paper is shown as red lines. Grey lines correspond to the description of the high frequency data according Eq. (S4).

Ensuring linearity
To warrant experiments in the linear range, we explore creep at different stresses and ensure that the stress chosen for our creep experiments ( 0 = 0.2 Pa) is within the range of stresses where the compliance remains unchanged, see Fig. S4(a). For the choice of the strain used in stress relaxation and oscillatory strain experiments, we measure ′ and ′′ at a fixed frequency of 5 rad/s and vary the amplitude with results shown Fig. S4(b). We chose to work with 0 = 0.1%, which is a strain within the linear range, while ensuring sufficient torque to obtain reliable data at all frequencies investigated. A very powerful cross-check of linearity consists in verifying that the following relation holds [12] � ( − ′) ( ′) ′ ∞ 0 = . (S6) As shown in Fig. S4(c), the convolution integral on the left-hand side computed from the experimental data of and closely matches the line = . This shows unambiguously that both and are probed within the linear response regime at all times during the tests. denotes the onesided Fourier transform [12,13]. As shown in Fig. S4(d), the data obtained are rather noisy at high frequency due to the limited time resolution of our stress relaxation test. At lower frequencies however, the data obtained by Fourier transform are in good agreement with the frequency sweep data.
Let us finally remark that our approach based on fitting the creep compliance to recast the lowfrequency relaxation has many advantages, as shown in Fig. S4(d) and Fig. 1(b) of the main text: -the results extend to a mid-frequency range (0.1 − 1 rad/s) where the simple transform ℱ[ ( )] fails due to limited time resolution in relaxation tests. -the development of ′ nicely connects to the high frequency description by Eq. (S4).
-the possibility to analytically treat the description of the creep data enables us to compute the results expected for → ∞, which are otherwise not accessible. -most importantly, while the transform ℱ[ ( )] works satisfactorily at low frequency, it does not convey any physical information about the processes involved. By contrast, our approach based on creep and recovery highlights the role of delayed elastic storage at low frequencies.