Rubber-like elasticity in laser-driven free surface flow of a Newtonian fluid

Significance Upon deformation, Newtonian fluids are expected to exhibit viscous behavior, and only when deformed on very short timescales, below the molecular diffusion time of a single molecule, is a solid-like elastic response expected. We have revealed a strong, rubber-like elasticity in the Newtonian fluid glycerol by analyzing the dynamics of a laser-driven free surface bubble. Not only do we find an elasticity persistent for four orders of magnitude longer than the diffusion time but also observe tolerance to large deformations only found in rubber-like materials. Our observations are independent of surface tension and require the existence of a transient state with solid-like long-range correlations different from the bulk state. This invites us to revisit our understanding of the liquid state.


A. Light absorption in water and glycerol
Water and glycerol have strong light absorption in the infrared spectral range around 3 µm wavelength (3400 cm -1 wavenumber), see Fig. S11. The absorption coefficient at 2.94 µm wavelength and with our experimental ablation fluences is approximately 8300 cm -1 for water (1). For glycerol, fluence-dependent absorption coefficient information is lacking in the existing literature and the best available value is 6500 cm -1 (2). The strong absorption is based on resonant excitation of O-H bond stretching modes prevalent in both compounds. Due to the strong interaction between the molecules by means of hydrogen bonding, the excitation energy redistributes quickly over the entire ensemble. For a comprehensive absorption spectrum of water from 40 nm to 1 mm wavelength, see Fig. 1 of ref. (3).
Details on the interaction between a sub-nanosecond laser and water at high fluences in ambient conditions can be found in ref. (4). The fluences used for such ablations are five times higher than those necessary for bubble formation. A higher fluence results predominantly in phase explosion instead of photomechanical spallation (see Sec. B). Similarly, numerous studies using mid-IR nanosecond lasers have tackled the laser ablation of liquids under ambient conditions, particularly for water or water-rich biological targets such as tissues (5)(6)(7). Cheng et al. (5,6) employed a 1-D fluid dynamics model to simulate plume dynamics after ablation, building on the experimental work of Aptiz and Vogel (8), who demonstrated ablation with a 70 ns, 2.94 μm laser. However, these studies also primarily concentrated on the phase explosion process, which dominates at higher fluences, and not the photomechanical spallation, where we observe laser-driven bubbles.

B. Photomechanical ablation efficiency and local temperature increase
Photomechanical effects can alter the efficiency of material ejection (9, 10) during general laser ablation events. In a pure photomechanical ablation process, a thin layer near the liquid-vacuum interface is heated by a short laser pulse to temperatures insufficient for vaporization but at such a high rate that it cannot undergo thermal expansion. This creates a traveling compressive stress wave. Due to a notable difference in acoustic impedance between the liquid and vacuum, the wave reflects at the free surface. As a consequence, a tensile component can form and grow in amplitude with increasing depth. If this tensile wave amplitude surpasses the material's tensile strength, a thin layer can be detached (11,12). A schematic description of such an event is available in Fig. S12. The efficiency of photomechanical ablation depends on the material's properties, especially its tensile strength, absorption length, and the Grüneisen parameter. The Grüneisen parameter determines the initial shock pressure being induced in the system (11,12).
The occurrence of photomechanical effects can increase the ejection efficiency of ablation, which means that of the available pulse energy, less is converted into heat resulting in a reduced temperature in the target. We are unaware of studies where the local temperature increase was probed in detail in the case of pulsed laser ablation of liquids. The presumably transient nature of the associated material parameters during laser interaction and afterward will make the problem particularly challenging to tackle. There is little consensus regarding the Grüneisen parameter and the dynamic tensile strength, especially for glycerol, at these extreme conditions. Thus, we will perform a best-effort estimation of the maximum shell temperatures in the system. We can deduce a qualitative understanding of the temperature increase during the stress-confined ablation utilizing the isochoric heat capacity . While this parameter is not at our disposal for glycerol at the relevant pressure ranges (>100 MPa), we rely on recent work by Ahmadi et al. (13) at a pressure of 55 Mpa, where ≈ 2142 J/(Kg·K) is reported for a 298 K isotherm. With knowledge of the volume of material interacting with the laser pulse, the temperature increase can be calculated. The interaction volume is given by the area over which a bubble is formed times the depth where 90 percent of the laser energy is absorbed. For an absorption coefficient of μ = 6500 cm −1 , 90% of the laser energy will be absorbed within 3.5 µm. For this cylindrical volume of 3.5 µm height and radius of 80 µm, 43 µJ will be deposited by the laser pulse with a peak fluence of 240±11 mJ/cm 2 . Assuming all 90% of the laser energy is converted to heat, we get an upper estimate of the average local temperature increase of about 200 K resulting in an absolute temperature of about 500 K. This is well below the boiling point of glycerol (≈ 560 K) (14) and suggests a liquid shell, as supported by our observations. However, due to the exponential absorption of laser energy, layers up to a depth of less than 1.5 µm could be heated above the boiling point in this limiting scenario.

C. One-Dimensional solution of the thermoelastic wave equation in a liquid and minimum spall fracture depth for glycerol calculation
A detailed discussion regarding the photomechanical spallation process can be found in the article by Dingus and Scammon (11). The important takeaway for us from this work is the qualitative prediction they have concerning the spallation depth. They estimate the spallation to occur at a depth one or a few times the absorption length of the ablation laser. In our case, the absorption length is around 1.5 µm (for μ = 6500 cm −1 ), giving a rough spallation depth of a few micrometers.
To better understand the process, we used a simple model describing the photomechanical responses in liquids after laser pulse absorption, which follows ref. (15). Using this model, a qualitative understanding of the transient stresses induced during our ablations can be gained. To this extent, we follow a slightly modified process described by Paltauf et al. (12) for top-hat pulses (spatial domain) to estimate the thickness range of the spall layer in our pico-second laser-induced ablations.
In the case of stress-confined ablations, when for the pulse width tp and the acoustic dissipation time tac holds tp << tac, the absorption of the top-hat laser pulse produces an initial pressure distribution 0 ( ): Where Γ is the Grüneisen parameter of the liquid, relating the initial pressure to the absorbed volumetric energy, 0 is the fluence of the laser, μ is the liquid's absorption coefficient for the ablation laser's specific wavelength, and d is the depth in the sample. In the case of liquid glycerol at room temperature, we find Γ = 0.9 (16) and μ = 6500 cm −1 (2).
Following the absorption, compressive stress waves are generated. These generated waves consist of two parts with the same spatial profile as 0 , but half the amplitude, traveling in opposite directions. The part propagating towards the free surface is reflected at the surface with a reflection coefficient . Z1 and Z2 are the acoustic impedances of the liquid (for glycerol 2.4*10 6 kg/m 2 s (17)) and the surrounding medium, respectively. In our case, the surrounding medium is a vacuum. Thus, irrespective of the acoustic impedance of the sample, we have the reflection coefficient = −1. Then the amplitude δ ( , ) of the stress wave at a depth d > 0 can be written as a sum of 3 pressure components 1 , 2 , 3 , Where is the speed of sound in the medium for glycerol ≈ 1890 m/s, and 0,max = μ Γ 0 .
The influence of a finite laser pulse duration can be incorporated by convolution with the temporal profile of the laser for a revised pressure amplitude as a function of time and depth: Our modeling uses ( ) as a Gaussian temporal distribution of sigma ≈ 170 pico-seconds (corresponding to a FWHM of 400 ps) centered at around 600 pico-seconds. We do not know when the spallation is happening as this is dependent on transient parameters of the liquid outside of these equations, like a temperature-dependent spallation strength. Hence, the modeling is carried out for a broad temporal window of 8ns (Fig. S13A).
We calculated the time-dependent hypothetical spall depth (HSD) (Fig. S13B), which is simply the depth at which the pressure amplitude reaches the dynamic tensile strength (18) of glycerol (85 MPa at 293 K). At these given conditions, HSD can reach a depth of up to 12 µm within 8 ns. The lowest depth at which a potential spall fracture could happen is around 300 nm very early after absorption. Given the temperature-dependent variations of Γ and fluence-dependent variations in these values may vary significantly. However, a clear behavior concerning the spatial fluence distribution and its effect on the shell thickness can be garnered from the calculations. In this model, lower pulse energies lead to lower fluences for the top hat pulse and an increased HSD (Fig. S13B). In our experiment, we observed higher shell thickness for the parts closer to the bulk and thinner shells for the top features of the plume. This observation can be attributed to the Gaussian spatial intensity distribution of our ablation pulse, where the 0, is not constant in the spatial domain. The top-hat pulse used in modeling fails to account for this experiment's boundary conditions. With the HSD calculations available for different pulse energies (Fig. S13B), we can infer that the wings of the Gaussian spatial intensity distribution with lower fluence likely result in higher HSD and thicker shells.

D. Relation of the observed dynamics to laser-induced forward transfer (LIFT)
The LIFT process (19) is a laser-driven material ejection process where also bubble formation is observed (20,21). However, comparing the LIFT dynamic with our process more rigorously, the similarities seem superficial. The mechanism of material ejection due to photomechanical stress waves differs from the phase explosion induced material ejection in LIFT. In the latter, the material is driven by an expanding vaporized volume, while photomechanical spallation happens ideally under the absence of vaporization, and only a minimal amount of energy contributes to phase change, apart from some evaporation within voids and cracks (12).
Between LIFT and photomechanical spallation exist considerable differences in typical values of parameters and observables regarding laser wavelength, pulse length, and laser energy. Furthermore, there are differences in the geometry of the ablation process and the induced physical processes. Typically, the LIFT process is performed in a back-illumination geometry, where a highly focused laser enters the backside of a target and creates an ablation of the front side of the target, see Fig. 5 of (19). The ablation of the front is initially driven by a pressurized volume of vaporized target created in some depth below the front surface.
Our experimental setup (SI Appendix, Fig. S6) is realized in front illumination, where the target is hit from the front, and the front section is ablated. The ablated material propagates in the direction from which the laser came. This geometrical difference implies that the ablated material has interacted with the highest intensities of the laser. At LIFT, the ejected material is from a region that did not experience the high fluences present in the vaporization region.
The lasers employed in LIFT, particularly in experiments (20) with liquid targets such as glycerol, typically have nanosecond pulse lengths and wavelengths in the shorter part of the visible spectrum, e.g., around 350 nm and below. If we consider the absorption spectrum of water, which we find sufficiently similar also to the one of glycerol, we find the region around 350 nm to be the one with the longest absorption length and hence weakest absorption. In fact, the difference in absorption at 350 nm and 3000 nm, where we operate in our setup, is eight orders of magnitude (see Fig. 1 of ref. (3)). In our experiment, the matter with the most intense interaction with the laser is the most relevant for the observable dynamics. In LIFT, the matter with the least interaction with the laser light forms the matter of the early plume dynamic.
In LIFT, the laser is highly focused, allowing for fluences that are 100-1000 times larger than in our case (20,21). As a result, the pressure generated during this phase explosion driven bubble expansion is generally one order of magnitude higher than that due to the mere thermoelastic effect (22). These differences in laser parameters, energy absorption, and geometry lead to considerable differences in the observable dynamics. Our dynamics are pretty much covered by a few microseconds, while in LIFT, the dynamics last several dozens of microseconds and hence ten times, or more, longer.
The LIFT process is typically carried out under atmospheric conditions. Our process is not realizable if not carried out in a vacuum environment. The shell is only formed if not subjected to any backpressure. At atmospheric pressure, the environmental air molecules immediately destroy the shell due to the creation of many nucleation sites from air-liquid interaction.

E. Kelvin-Voigt and Maxwell rheological models
Rheological models are frequently visualized by a network of springs, resembling the elastic contribution, where the force is linearly dependent on the elongation following Hooke's law, and dashpots, resembling the dissipative viscous contribution, where a Newtonian fluid is assumed and the force is proportional to the temporal change of elongation.
The simplest manifestation of the Maxwell model (23) contains two hypothetical elements, a single dashpot in series with a single spring. In such a series network, the stress in the dashpot is always identical to the stress in the spring. The total strain , however, is always the sum of the individual strain in each element. If this combination is suddenly deformed, then the stress decays on a characteristic timescale , the relaxation time given by the viscosity of the dashpot divided by the elastic constant of the spring. It shows the viscous flow for a timescale much longer than the relaxation time and elastic resistance as an initial response to fast strains. The constitutive equation is given by: The generalized Maxwell model is a more complex network of parallel models. In particular, a single spring is parallel to a number of n Maxwell models (each with spring and dashpot in series) like Eq. S7 and generally different viscosity, elasticity, and hence relaxation times. Due to the single spring in parallel, it can have elasticity for infinite times.
The Kelvin-Voigt model (23) is as simple as the Maxwell model and is given by a parallel network of spring and dashpot, respectively. In the parallel network, the strains, not the stresses, are identical in each branch and identical to the total strain. The total stress is the sum of the individual stresses in each branch. We find a constitutive equation of the form: Suppose we would add a mass point to the otherwise mass-less network and solve for its motion Eq. S8 can be transferred to the standard form of the harmonic oscillator: with new parameters and 0 representing a dimensionless damping ratio and the frequency of the periodic motion, respectively. These parameters are, of course, functions of parameters of Eq. S8 and . The well-known solutions of Eq. S9 with amplitude and phase are: The velocity (and ) dependent term 2 0̇ of Eq. S9 leads to a sinusoidal motion with damped amplitude and altered frequency.
The Kelvin-Voigt model shows instantaneous viscous behavior since the straining of the spring and the dashpot, which must work for every deformation, co-occurs. The Kelvin-Voigt model is a simple model for solids with viscosity and less meant for fluids. For a constant force, the deformation is damped by viscosity and the system will finally reach a state of pure elasticity. The Kelvin-Voigt model is an entry point for describing an elastic solid rather than a liquid subject to viscous energy dissipation.
While the Maxwell element constitutes a high-pass with low attenuation for fast dynamics, the Kelvin-Voigt element resembles a low-pass with high attenuation for fast dynamics.
We simulated the behavior of the membrane under the assumption of a Maxwell rheological model. For that, we kept the spring's elastic constant, here in the form of the shear modulus ′ , constant and used different parameters for the viscosity and relaxation time. We used ′ = 3 * 10 6 Pa, which is about ten times larger than our experimentally found elastic constant and reduced the membrane's kinetic energy by going to a tenth of the velocity. This ensured we could reach the quarter period point within the simulation, at which we have maximum elongation. The challenge lies in the significant shape changes of the simulation domain.  Fig. 2B and C. We find that for different viscosities and hence relaxation times according to = * ′ the motion of the mass point follows a sine-like behavior at least for the time till the maximum position is reached. This is true for the fully elastic membrane but also allows Maxwell models with different viscosities to produce a motion of , that is potentially identifiable with a sine-like appearance, given the present experimental resolution. What these Maxwell models separate from an actual harmonic behavior becomes evident if we inspect the distribution of the total energy to the branches of kinetic and elastic energy. Fig. S14B shows exemplarily the distribution of energies for the simulation with = 1.4 Pa*s. While in an actual harmonic motion, the energy at maximum elongation is entirely converted from kinetic energy to elastic potential energy, for the Maxwell model with = 1.4 Pa*s, the energy is not stored but got dissipated and, as such, is not available for any dynamic, like a hole opening dynamic at high velocity or a jetting of matter. Such a situation can be experienced in the slow blowing of chewing gum bubbles (24). Once the driving pressure is removed, the shell is dynamic-less and just subject to gravity.

F. Dimensionless flow numbers
Surface flows can be categorized by dimensionless numbers (25) expressing the relative influence of different forces such as capillarity (surface tension), viscosity, inertia, or gravity. The capillary number (Ca) expresses the relative importance of viscous drag forces over capillary forces. High Ca indicates the dominance of viscous forces over capillary forces. The capillary number is defined by the viscosity , a characteristic velocity and the surface tension : An elasto-capillary number suitable for the present experimental conditions can be defined by following ref. (24) and utilizing the shear modulus : relates elastic forces with capillary forces expressed by the surface tension. expresses again typical dimensions of the flow and is used in analogy to the radius in a spherical shell.

G. Sine-like tip motion for shells under the Maxwell model
Not only fully elastic membranes can exhibit a sine-like motion for the tip, but also shells assumed to be of a material following the Maxwell model might. Under the Maxwell model (SI Appendix, Sec. E), the material has an elastic constant G' and a viscosity η. From these two parameters follows directly a relaxation time = / ′. Fig. S14A shows simulation results of the tip position for different sets of parameters G', η. The results of the different simulations can be reasonably fitted by simple sine functions. The difference to the sine motion in a pure elastic shell case is in the temporal development of the kinetic and elastic energy in the shell. Fig. S14B shows that for the Maxwell model, the kinetic and the elastic branch of energy in the shell can be (entirely) depleted in the sine's quarter period, leaving no energy in the shell to perform any dynamic past the quarter period point. For the fully elastic shell all initial kinetic energy would be readily available in the form of elastic energy at the quarter period point, e.g., for a fast hole opening dynamic.

H. Lift-off
The lift-off is a defect of the shell in the lower third part, close to the bulk, where the prevalence to multiple nucleation and rupture reduces the area of the shell that is connected to the bulk. Finally, this can largely detach the shell from the bulk. Our work shows that higher straining rates coincide with higher elasticity, stabilizing the shell. The region where the lift-off occurs has higher shell thickness and lower strain rates than other areas.
Both findings are adverse for a strong elasticity and a stabilization of the shell. Additionally, the droplets traveling in the volume encompassed by the shell are potential sources of rupture once they get in contact with the shell. While this contact seems unlikely to happen in the upper region of the shell due to the lower velocity of the droplets in comparison to the shell, it seems likely in the region where the rupture for the lift-off takes place. The multiple rupture sites in the lift-off region, where holes also dynamically combine to form larger holes prevent us from performing a statistical analysis to determine the hole growth velocity and infer elastic parameters. We cannot make a comparison to the holes in the upper section and get information for different sections of the shell.

I. Storage modulus for liquid glycerol at different frequencies
The idea that there might be a measurable storage modulus ( ') in Newtonian liquids well above the molecular relaxation time scales is nothing new (26). In Table S3, we present the studies which have looked at the storage modulus ' of glycerol at different frequencies.

J. Traveling rim
The rim's origin, shape, and velocity under varying influence of viscous, inertia, or surface forces draws steady attention (28)(29)(30)(31)(32)(33)(34)(35)(36). For circular hole growth, we can follow ref. (32), where the authors differentiate three regimes depending on the Ohnesorge number ℎ = /√2 relating viscosity to inertia (density , film thickness ) and surface tension . For high ℎ > 10 it is found that the toroidal rim is not forming; instead, the whole film thickness is increasing over a certain region. Under our experimental conditions, the Ohnesorge number is ℎ exp ≈ 3. If we assume a 10-fold increased viscosity, the Ohnesorge number increases by the same factor and we would enter a parameter region where we can expect the influence of viscosity on the shape of the rim and the velocity with which the rim travels. An exponential velocity growth over time towards the Taylor-Culick velocity should develop, which is not observed. However, caution is advised in this conclusion. We have shown that the shear elasticity can disguise as apparent surface tension. This apparent surface tension can be calculated by solving the expression of the Taylor-Culick velocity for the surface tension: (S14) With our parameters, we find an apparent surface tension about 320 times larger than the bulk surface tension. The Ohnesorge number would decrease by almost a factor of 20; hence, the dynamic viscosity could increase by another factor of 20 while still having Ohnesorge numbers of 10 and below. Hence, a much more significant dynamic viscosity increase during the plume dynamic could be possible, although the rim shape suggests conditions of an inviscid liquid.

K. Holes fracture
The hole fracture, which we analyzed to infer the rim speed, happens in the upper half of the shell and generally at times when the expansion of the shell has already slowed down. Nevertheless, an underlying expansion is potentially present during the hole growth dynamic. In the analysis of the rim speed, we neglect any influence of this underlying motion on the rim speed. Fig. S15 depicts a small section of the expanding shell and the velocity components near the onset of a hole fracture. The cyan disc indicates the hole opening. The mass points in close vicinity to the rupture have different velocity components. A velocity component indicated by the red arrows reflects the outward motion of inflating volume and is parallel to the surface normal and perpendicular to the rim motion. Then there is a velocity component indicated in green that reflects that the surface is increasing. For an observer in the point of rupture, all other points on the surface move away and increase distance. This effect is entirely independent of any traveling rim. The red velocity field, perpendicular to the motion of the traveling rim, overlays with the rim motion without interference and effect on the magnitude of the rim motion in the surface plane. The green velocity components are small in magnitude compared to the rim speed, which is 60 µm/µs, since the green motion is the relative motion of neighboring points. The maximal velocity difference between any two points on the entire shell, at any time of the dynamic, is only a fraction of the tip's maximal 300 µm/µs velocity since the shell is a flat disc in the beginning that bulges. Hence, the rim's traveling speed is very high compared to local differences in the shell's velocity, even for well-separated points on the shell. That is why the holes finally "eat up" the entire shell.
Furthermore, the green velocity components can be understood in terms of how they increase the stress in the shell because the stress is finally responsible for the speed of the rim. We have shown in the report that the elasticity in the shell can be expressed as an apparent surface tension, as expressed in Eq. 11 in the Materials & Methods section. At some point of inflation, this surface tension becomes independent from further extension due to power scaling in the formula. Consequently, the forces do not increase, although the extension increases. So, even if the velocity components do increase the general extension of the shell, the influence on the forces is small and the velocity might only be influenced by the fact that a stretched and hence thinner shell has a rim with less mass at a given hole diameter. This change in thickness is again a slow process compared to the dynamic of a traveling rim.
In total, the effect of the velocity components of the shell is mainly that they dynamically affect the topology the rim is traveling on, but not relevantly the velocity with which the rim travels on the surface of the shell. The underlying expansion changes the time it takes for the rim to reach a point on the shell at a certain distance to the fracture point but not the velocity with which it travels the potentially larger distance toward this point.

L. Elasticity in charged liquids
Charging and potentially related strong electric fields have been identified in static experiments to promote elastic properties (37). Short pulse laser-matter interaction can potentially charge the liquid via plasma creation. However, the use of 3 µm wavelength radiation greatly obfuscates the required multi-photon interaction. We do not detect light emission accompanying any plasma creation at the employed laser intensities, although detectable in our setup at considerably higher energies.               Movie S1 (separate file). Temporal evolution of the laser-induced bubble in glycerol. The bubble's temporal evolution for the first 3 µs under vacuum conditions of 10 -3 mbar, initiated by an ablation fluence of 240±11 mJ/cm 2 . The movie-like series of snapshots is not the temporal evolution of a single bubble but the temporally ordered sequence of individual snapshots of different ablation events recorded at different time points spanning over 3 µs.