Causal Discovery for Materials Degradation

Most engineered materials degrade over time, and the measurements that best predict how fast they fail are usually slow and expensive to run. This project finds which cheap, fast measurements actually drive degradation, rather than which ones merely correlate with it. We built and verified a physics simulator, trained a neural-operator surrogate that reproduces it to within about 1.5% error, and designed a five-part audit that catches models that predict well for the wrong reasons. On a synthetic testbed where the correct answer was known in advance, the audit identified the true driver and flagged the decoy in all five independent tests. Current work extends the method to experimental degradation data.

• ~100× fasterDataset generation cut from about 14 hours to about 25 minutes by a float32 CUDA solver rewrite (GPU versus CPU).
• 1.3 to 1.5%Relative L2 error from the FNO surrogate, in and out of distribution.
• 5 of 5Causal diagnostics agreed with a per-sample ground-truth key across 8,000 simulations.
• 3 schemesIndependent solver schemes in agreement before any training data was generated.

The problem: causes, not correlations

Almost all engineered materials degrade, and identifying which measurable properties govern the rate of degradation is central to designing durable systems. At present these properties are usually established through slow, costly characterization. This work identifies which cheaply and quickly measured features are the strongest causal drivers of degradation, and verifies that a predictive model depends on genuine causes rather than on correlations that only track the outcome in the training data. The distinction matters because a model that predicts well for the wrong reasons fails precisely when it is extrapolated to a new material.

Before touching experimental data, we built a fully controlled synthetic testbed in which the correct answer is known in advance. Every claim below was checked against that ground truth.

A solver verified before any training

The synthetic system is the Allen–Cahn equation with a spatially varying mobility $M(x)$, which models how a two-phase microstructure evolves:

Here $u$ is the order parameter, $\varepsilon$ the interface width, and $g$ a constant tilt of the double-well potential. We used a second-order semi-implicit spectral scheme, cross-checked against an independent fourth-order Runge–Kutta solver and a py-pde reference. A first-order stepping scheme missed the cross-scheme tolerance, so we moved to a second-order SBDF2 scheme, which brought the schemes into agreement at $7\times10^{-5}$.

Verification came before any training data. The Ginzburg–Landau energy decays monotonically, and the shrinking-circle test recovers the curvature law

to within 0.8% of theory. The coarsening length obeys $\ell^{2}\propto \varepsilon^{2}t$ with $R^{2}=0.979$. The textbook exponent $\ell\sim t^{1/2}$ is unmeasurable in a unit box, so we restated that check in its reachable differential form rather than pretending to measure it.

Data generation was the first real bottleneck. A single-threaded NumPy implementation projected roughly 14 hours per dataset. Rewriting the solvers in float32 CUDA gave approximately a 100× speedup, to about 25 minutes, with the float64 NumPy path retained as the verified reference.

A neural-operator surrogate

Even at CUDA speeds, generating data at scale with the solver is expensive, so we trained a Fourier neural operator (FNO) that maps the input fields directly to the final state. Two decisions shaped it, both learned the hard way.

The first was horizon length. The initial surrogate used a horizon of $T=4$ and memorized its training set while generalizing at only 17%, because a one-shot map across an entire coarsening cascade is poorly conditioned. Shortening the horizon to $T=1$ resolved this; longer horizons are left to autoregressive composition.

The second was a hidden input. At $T=1$ the model still overfit at 11% error because $\varepsilon$ was not supplied as an input, so the target was not a well-defined function of the inputs. Adding $\varepsilon$ as a constant channel reduced the error roughly sixtyfold, to 1.5%. That same hidden-variable structure later became a deliberate design element of the causal audit.

The final surrogate reached 1.52% relative $L_2$ error in distribution and 1.26% on a held-out parameter band. Because the two independent solvers agree to within $\sim\!9\times10^{-4}$, essentially all of the residual error is attributable to the model rather than the solver.

Auditing what the model actually relies on

The central question is whether a model that predicts well has learned the causes of the dynamics or a spurious correlate. To answer it with certainty, we constructed a dataset with a known ground-truth key: a nuisance factor $S$ that is provably non-causal and a driver $C$ that is provably causal, with training environments in which $S$ correlates with the outcome and evaluation environments in which that correlation weakens, breaks, or reverses.

The construction hides the causal driver, because the mobility field $M(x)$ is unobserved, so the nuisance carries genuine predictive information in the biased environments. In the nucleation-growth regime with a tilted well, traveling-wave analysis gives a front velocity

and a critical radius that is independent of the mobility amplitude. Seeding every nucleus above this radius makes severity monotonic in the hidden amplitude. The full pool contains 8,000 simulations across seven environments, each with a per-sample answer key.

Five independent diagnostics assess causal against spurious reliance, each ruling something in or out, and each checked against the key.

Across all five diagnostics the verdicts agreed with the key. The clearest signal is the Sobol sensitivity: the nuisance factor is negligible in the physics yet dominant in the model, so the model's sensitivity structure inverts the true physics. The latent-ablation step doubles as a mitigation: removing the identified nuisance direction restored predictive validity on adversarial data while leaving the decoded physics roughly 96% intact.

Two diagnostics were initially mis-specified and corrected without changing any threshold: a multivariate stability test was replaced with ICP-correct univariate probes, and the IRM check was reformulated around out-of-environment risk. One secondary check returned a negative result: an IRM-penalized probe generalized worse than plain empirical risk minimization at every penalty weight. We report this outcome as found.

From synthetic to experimental data

The synthetic testbed established that the audit reliably distinguishes genuine causes from misleading correlations and can repair a model that has adopted a shortcut. The current phase applies the same reasoning to experimental degradation measurements: training predictors on real data, running the diagnostics to identify which features are load-bearing, and ranking candidate features by their causal effect on the degradation rate. The intended outcome is a clearer account of which properties most influence degradation resistance, and which high-impact features low-cost sensors should be designed to monitor.

Spectral Signatures of Lahore's Air Pollution

Lahore has some of the worst air quality in the world, and pollution that thick leaves a fingerprint on the light that passes through it. This project reads the sky's visible and near-infrared spectrum for that fingerprint. Using a low-cost spectrometer pointed at open sky, we measured downwelling spectra under two very different pollution levels and decomposed each into its spectral components, then proposed reading those components as signatures of the pollutants in Lahore's air.

• AQI 211 vs 420Sky spectra measured under two contrasting air-quality conditions in Lahore.
• 340 to 1180 nmVisible to near-infrared range, captured with a StellarNet BLUE-Wave spectrometer.
• 47 and 43 peaksMulti-Gaussian deconvolution of each spectrum into its spectral components.
• FrameworkA method proposed for reading spectral features as signatures of specific pollutants.

The idea: reading pollution off the sky

When sunlight travels through the atmosphere, gases and particles absorb and scatter it wavelength by wavelength, so the spectrum that reaches a spectrometer on the ground is shaped by everything it passed through. To first order, the intensity at each wavelength is attenuated according to the Beer–Lambert law,

where the optical depth $\tau(\lambda)$ collects the contribution of every absorbing and scattering species along the path,

Each species has its own absorption cross-section $\sigma_i(\lambda)$, so a change in the mix of pollutants should change the shape of the measured spectrum. That is the premise: if Lahore's pollution imprints measurable features, then comparing sky spectra taken at different pollution levels is a way to look for them. The measurement adapts PhysLab's AstroLAB blue-sky spectroscopy lab, which normally treats the blue sky as a stand-in for the solar spectrum, and repurposes it toward an air-quality question. We proposed this framing and the analysis plan.

Measurement in the field

Spectra were collected with a StellarNet BLUE-Wave spectrometer coupled by a fibre-optic cable to the SpectraWiz software. Following the AstroLAB method, the fibre fed a simple diffuser tube whose open end was covered with foil and pricked with a pin, giving a small, uniform aperture. Mounted on a tripod and pointed at an open patch of sky away from the sun, this records the downwelling sky spectrum. Integration time was tuned to sit between undersampling and saturation. We ran the acquisitions under two sharply different reported air-quality conditions, at AQI 211 and AQI 420.

Decomposing each spectrum

Each spectrum was baseline-corrected with an end-points-weighted baseline and then decomposed into a sum of Gaussian components through peak-analysis fitting, so the broad envelope and the finer features could be examined separately. The AQI 211 spectrum was fit with 47 Gaussians over 1684 points; the AQI 420 spectrum with 43 Gaussians. In both, the envelope of scattered skylight rises through the blue and green and falls away toward the near-infrared.

Several fitted components near 700 to 760 nm are negative, consistent with a sharp absorption band in that region; the molecular-oxygen A-band near 760 nm is the natural candidate. Features like this are the ones the method is built to pick out.

Matching peaks to functional groups

Once a spectrum is decomposed, each Gaussian gives a peak centre $\lambda_k$, and the question is which chemical species each centre belongs to. The anchor is again Beer–Lambert: absorption at any wavelength is a sum over species, each with its own cross-section $\sigma_g(\lambda)$ and column density $N_g$,

Two readings fall out of this. The direct one is assignment: given a library of reference band positions $\{\lambda_g\}$, each peak is matched to the nearest band within a tolerance $\Delta$,

with a match score $s_k = \exp\!\big(-(\lambda_k - \lambda_{g^{*}})^2 / 2\delta^2\big)$ that weights closer matches more heavily. The stronger reading is unmixing: treat the measured absorption as a non-negative combination of reference cross-sections and solve for the column densities directly,

In the visible and near-infrared the candidate bands include the molecular-oxygen A-band near 760 nm, the ozone Chappuis band across 500 to 700 nm, nitrogen-dioxide absorption below about 500 nm, water-vapour bands through the near-infrared, and O–H and C–H overtone bands past 900 nm. Paired with real-time pollutant readings and a regression or principal-component step, these assignments become a quantitative link between a spectral feature and the pollutant behind it.

The same plotting pipeline used on the raw spectra extends directly into the assignment step:

import numpy as np
import matplotlib.pyplot as plt

# SpectraWiz export: two columns, wavelength (nm) and count
wl, count = np.loadtxt("sky_211aqi.txt", unpack=True)

plt.figure(figsize=(7, 5))
plt.plot(wl, count, color="0.2", lw=0.8)
plt.xlim(340, 1180)
plt.xlabel("wavelength (nm)")
plt.ylabel("count")
plt.title("Sky spectrum (AQI 211)")
plt.tight_layout()
plt.savefig("spectrum_211.png", dpi=150)

% SpectraWiz export: two columns, wavelength (nm) and count
data  = load('sky_211aqi.txt');
wl    = data(:,1);
count = data(:,2);

figure;
plot(wl, count, 'k', 'LineWidth', 0.8);
xlim([340 1180]);
xlabel('wavelength (nm)');
ylabel('count');
title('Sky spectrum (AQI 211)');

import numpy as np

# Gaussian peak centres (nm), from the deconvolution "Center Grvty" column
peaks = np.array([369.95, 417.53, 448.64, 461.13, 468.03, 512.71,
                  583.65, 649.03, 724.90, 735.25, 760.00])

# Reference bands: name, centre (nm), tolerance (nm)
reference = [
    ("NO2 (electronic)",  430.0, 40.0),
    ("O3 Chappuis",       600.0, 80.0),
    ("O2 B-band",         688.0, 10.0),
    ("H2O (NIR)",         720.0, 20.0),
    ("O2 A-band",         760.0, 10.0),
    ("H2O (NIR)",         820.0, 25.0),
    ("O-H overtone",      970.0, 30.0),
    ("C-H overtone",     1150.0, 40.0),
]

def assign(peak, ref):
    """Nearest reference band within tolerance, with a Gaussian match score."""
    best, best_d = None, np.inf
    for name, lam, tol in ref:
        d = abs(peak - lam)
        if d <= tol and d < best_d:
            best, best_d = (name, lam, tol), d
    if best is None:
        return None
    name, lam, tol = best
    score = np.exp(-(peak - lam) ** 2 / (2 * (tol / 2) ** 2))
    return name, lam, best_d, score

for p in peaks:
    m = assign(p, reference)
    if m:
        name, lam, d, score = m
        print(f"{p:7.1f} nm  ->  {name:16s}  "
              f"(ref {lam:6.1f}, off {d:4.1f} nm, score {score:.2f})")
    else:
        print(f"{p:7.1f} nm  ->  unassigned")

Long-Range Solar-Powered UAV

An aircraft that harvests more energy in daylight than it spends flying can, in principle, stay up indefinitely. This project builds toward that: a four-metre glider carrying solar cells across its wing, sized to fly ten hours continuously with four of those hours after the sun has gone. The engineering question is whether the numbers close, whether a 2.1 kg airframe flies slowly enough that its cruise power sits below what the cells can gather and the batteries can hold. They do, with margin.

• 4 m span, 2.1 kgA wing loading of 2.1 kg per square metre at an aspect ratio of 16.
• ~6.5 m/s cruiseAbout 23 km/h, from level-flight equilibrium at the design mass.
• ~230 kmDistance covered over a ten-hour flight at that cruise speed.
• 4 h after dark118 Wh of stored energy against an estimated 12 W cruise draw.

Airframe

The aircraft is a built-up glider with a 4000 mm wingspan on a 250 mm chord, giving a wing area of about 1.0 m$^2$ and an aspect ratio of

A slender wing like this is what makes long-endurance flight possible, because induced drag falls as aspect ratio rises. The fuselage is 1790 mm long, the tailplane spans 800 mm, and the wing carries 5 degrees of dihedral outboard of a 2.5 degree inner panel for lateral stability without active correction. At an all-up mass of 2.1 kg the wing loading is only

about a fifth of a typical radio-controlled sport model. Light loading means slow flight, and slow flight is what keeps the power draw low enough for solar operation to work.

Building the airframe

The geometry above had to become parts. Wing ribs and fuselage formers were nested and cut from sheet plywood on a CNC router, which keeps every rib in a set identical and wastes little material. Identical ribs matter more on a four-metre wing than on a short one: a twist introduced by inconsistent ribs reads as a standing rolling tendency, and the flight controller would then spend the whole flight trimming it out.

The wing is a built-up structure rather than a foam core. Ribs are threaded onto the spar, aligned, and taped while the joints cure, then the leading edge and upper surface are sheeted so the skin carries torsion. Sheeting the top surface is also what makes the cell layout practical, because a cell laminated over open bays would flex with the skin and crack.

The fuselage is a composite shell laid up wet in a two-part female mould, so the outer surface takes its finish from the mould rather than needing to be faired afterwards. A moulded shell also keeps the nose stiff around the motor mount and the battery bay, which is where the concentrated loads are.

How fast it has to fly

Cruise speed is not chosen freely; it follows from the mass and the wing. In steady level flight lift balances weight,

which rearranges to give the speed the aircraft must hold at a given lift coefficient:

Substituting $m = 2.1$ kg, $S = 1.0$ m$^2$, and sea-level density $\rho = 1.225$ kg/m$^3$ gives a cruise speed that depends only on how hard the wing is working:

Taking $C_L = 0.8$ as the nominal cruise point puts the aircraft at roughly 6.5 m/s, or 23 km/h, and a ten-hour flight therefore covers about 230 km. The speed is low by design: because $V$ scales with the square root of wing loading, the light airframe flies at a walking pace, and the power required scales with $V^{3}$, so every reduction in speed pays back sharply in energy.

Whether the energy budget closes

Cruise speed sets the power draw. With an aspect ratio of 16 and a span efficiency $e \approx 0.9$, the induced drag coefficient at $C_L = 0.8$ is

and adding a profile and parasite contribution of roughly $C_{D_0} \approx 0.025$ for a built-up airframe of this kind gives a total near 0.039, so the lift-to-drag ratio is about 20. Drag at cruise is then $D = W/(L/D) \approx 1.0$ N, and the aerodynamic power required is

At a combined motor and propeller efficiency near 55%, that is about 12 W drawn electrically, with the flight controller, GPS, and telemetry radio adding a few watts more. The two 4S 4000 mAh lithium-polymer packs hold

which over a four-hour leg after sunset allows an average of about 30 W. The night leg therefore closes with roughly a factor of two in hand against the cruise estimate. In daylight the wing is tiled with 5 by 5 inch monocrystalline cells rated at 3.6 W each; a cell covers about 0.016 m$^2$, so even a partial tiling of the 1 m$^2$ wing generates tens of watts at peak, above the cruise draw and leaving a surplus to recharge the packs and pay back the energy spent climbing.

Power path and avionics

The cells feed a Genasun 12 V (4S) maximum-power-point charge controller, which conditions the array output and charges the batteries continuously in flight while supplying the motor and avionics. Propulsion is an Emax GT-2826 1090 kv motor on a 60 A electronic speed controller, sized for the climb rather than the cruise, since climbing is where the aircraft spends energy that daylight then pays back.

Long flights rule out manual piloting, so navigation runs on a Pixhawk PX4 flight controller with a three-axis gyroscope and accelerometer, a NEO-M8N GPS module locking up to ten satellites, and a 915 MHz telemetry link that returns flight data to a ground station. Four servos drive the control surfaces. The controller flies preset waypoint paths and can return to the landing site on stored energy once the light goes.

Flight planning

Because the aircraft flies a preset route rather than being hand-flown, the route has to be built and checked on the ground before launch. Missions are planned in the PX4 ground station: we lay a waypoint path over satellite imagery of the site, mark the launch and return point, and bound the area the aircraft is to cover, then upload the plan to the flight controller.

Flight testing

Calculations of this kind are only worth as much as the airframe that flies them, so the glider was flown to check that it trims and holds a stable, slow cruise at the predicted attitude.

Rupture Propagation in Interconnected Soap Films

Dip a slinky in detergent and pull it upright and it fills with soap films, one spanning each gap between successive turns of the wire. Puncture any one of them and the collapse does not stay local: it runs through the whole stack in under three seconds. Isolated soap films have been studied closely for over a century, but a connected stack behaves differently, because gravity drains liquid downward and leaves the films thicker at the base than at the top. That gradient makes the system sensitive to something a single film cannot be sensitive to, namely where the rupture starts. We swept four variables across 48 conditions and found that rupture position changes the propagation speed more than anything else we could control.

• 48 conditionsFour sweeps over concentration, extension, temperature, and rupture position, three runs each.
• 2.2× fasterBottom-initiated ruptures propagated at 35.1 cm/s against 16.2 cm/s for top-initiated ones.
• Near-constant 2.8 sTotal collapse time held flat from 30 to 70 cm, so front speed tracked the length of the stack.
• 23% slowerDoubling the soap-to-water ratio cut front speed from 14.1 to 10.8 cm/s at 30 cm.

Background: films between the turns of a helix

A soap film spanning two coaxial rings takes the shape that minimises its area, which is the catenoid, a surface of revolution generated by a catenary. In a slinky each film is bounded by two adjacent turns of the same wire, so the stack is a chain of catenoids sharing a single helical support. The helix supplies the connectivity; the individual films are catenoids, and it is worth being precise about that, since a helicoid is a different surface entirely and would require a boundary along the central axis.

What happens after a puncture was settled for a single film by Taylor and Culick, who showed that the hole opens at a steady rate set by the balance between surface tension pulling the rim outward and the inertia of the liquid it sweeps up. The resulting speed is

with $\gamma$ the surface tension, $\rho$ the liquid density, and $h$ the film thickness. Two features of this expression drive everything below. Speed rises as the film gets thinner, as $h^{-1/2}$, and it rises only as the square root of surface tension. Surfactant concentration enters through $\gamma$ by way of the Gibbs adsorption relation,

and temperature enters through the Eötvös relation,

which makes surface tension fall as the liquid warms. Prior experimental work has characterised the collapse of isolated catenoids and hemicatenoids under varying viscosity and surfactant loading, including the role of Plateau border friction. Our contribution is not the physics of a single film but a systematic parameter sweep on a connected stack, treating rupture position as a controlled variable rather than an experimental detail. Networks of films are well studied in the foam literature, so the novelty here is specific: a one-dimensional connected chain with a known thickness gradient along it.

Apparatus and timing

A plastic slinky of diameter 8.1 cm was taped to the base of a container holding one litre of water, so that it could not be lifted free and burst from below. Soap was added by syringe in 20 ml steps, giving soap-to-water volume ratios from 0.10 to 0.20. The slinky was immersed, rotated ten times to wet every coil, then drawn upward by hand to a height read off a metre rule clamped vertically beside the container. Puncturing the film at one end started the collapse, which was recorded by a camera fixed at constant distance and height.

Front speed was obtained by dividing the extension by the interval between the puncture and the last film collapsing, read from the video. That definition matters for interpreting the extension results, so we return to it below. For the heated runs the container sat in an electronic water bath, weighted down to stop it floating, with the metre rule remounted at the new liquid level to keep extension readings honest. Every condition was run three times, and uncertainties quoted here are one standard deviation over those three runs.

Concentration and extension

Speed falls as more surfactant is added and rises as the stack is pulled further open. Both trends hold at every level of the other variable. Fitting the individual trials, the concentration slope is negative at all five extensions and steepens from 27.6 to 140.0 cm/s per unit ratio as the stack is extended, with $p \le 1.7 \times 10^{-3}$ throughout.

The shape of the effect is as informative as its size. The relationship is not linear: a quadratic reduces the residual sum of squares by roughly an order of magnitude at 40 and 50 cm, so the response saturates as more surfactant is added, which is what one would expect once the interface approaches full coverage. Most of the change happens at the dilute end, and at 30 cm the trend is carried almost entirely by the lowest ratio. That same condition is the one whose films gave way before the stack could be opened past 50 cm, consistent with weaker films at low surfactant loading, and it marks the lower edge of the workable range for this geometry.

Transit time

The extension result deserves closer attention than we first gave it. Since speed was computed as extension divided by transit time, a constant transit time on its own would produce a proportional rise of speed with extension. Recovering the times shows that this is close to what happens: the interval for the collapse to cross the entire stack varies by 3 to 11 percent across a factor of 2.3 in length, and no concentration shows a significant dependence on extension ($p \ge 0.08$).

What the transit time does depend on is the solution. It rises monotonically from 2.05 s at a ratio of 0.10 to 2.83 s at 0.20. The finding is therefore that the stack collapses in a near-fixed interval set by the liquid rather than by the length of the structure, and that front speed rises with extension largely because the same interval covers more ground. Film thinning at larger extension is consistent with this through the $h^{-1/2}$ scaling, and a thickness measurement would tell us which of the two is doing the work. The answer is interesting either way. If the interval really is fixed by the solution rather than by the geometry, then some property of the liquid, its drainage rate or its viscosity, sets a characteristic timescale for the collapse of the whole stack that barely depends on how many films the stack contains.

Temperature

Heating the solution roughly doubles the propagation speed. Against the unheated room temperature of 16.5 °C, the 37.5 °C runs are faster by factors of 1.64 to 2.38, significant at every extension. The second heating step to 72 °C adds much less and is not resolvable at 50 cm, so the response appears to saturate above roughly 40 °C. Transit time falls from about 2.8 s to between 1.1 and 1.6 s.

The mechanism is worth setting out carefully, because the two terms in $v_{\mathrm{TC}}$ pull against each other as the liquid warms. The Eötvös relation makes surface tension fall with rising temperature, and since speed goes as $\gamma^{1/2}$, that term on its own predicts a slower front. The thickness term points the other way: lower viscosity and faster gravitational drainage at elevated temperature thin the films, raising $h^{-1/2}$. Our measurements are consistent with the thickness term winning, and by a wide margin. Separating the two is a well-posed next experiment, since interferometric thickness measurement at two temperatures would fix $h$ independently and let the predicted velocity be compared against the measured one rather than merely agreeing in sign.

Region of rupture

Rupture position produced the largest effect of anything we varied. Puncturing at the base rather than the top raised the speed by factors of 1.94 to 2.30 at matched concentration, from a mean of 16.2 to 35.1 cm/s. Because the puncture had to be made just above the base, the bottom-rupture runs sat about 5 cm above their nominal extension. That offset is small against the effect and can be accounted for: charging it against the measured extension slope removes 9 percent of the gap and leaves a ratio of 2.07, so extension alone does not explain the difference.

The direction of this effect is worth dwelling on, because the simplest thickness argument gets it backwards. Drainage leaves more liquid at the base, and it is tempting to read that abundance as the reason a collapse starting there runs faster. Taylor–Culick says the opposite: a thicker film recedes more slowly, which is precisely the direction the extension result depends on. What reconciles the two is that the thickness which matters is not the one at the puncture but the gradient the front travels into. A bottom-initiated front advances into successively thinner films and accelerates; a top-initiated front advances into successively thicker films and decelerates. This is the same thinning mechanism read along the axis of the stack rather than across a single film, and it makes a checkable prediction: the front should not move at constant speed, but should visibly gain or lose speed as it goes. Timing each film collapse separately, instead of taking one interval from puncture to final collapse, would test that directly and is the measurement we would most like to have.

Drainage

A separate observable, and one worth distinguishing carefully from the propagation measurements: when a film bursts the rim breaks into a fine mist of droplets, and that liquid then runs back down the spiral. We timed how long the stack took to clear at a fixed extension of 70 cm. The time rises steeply between ratios of 0.10 and 0.14, from 21.0 to 38.6 s, then flattens, gaining only 8.4 s over the rest of the range. A linear fit returns a respectable coefficient of determination but is plainly mis-specified, since its residuals are systematically curved. We report the trend as monotonic and saturating and do not fit a predictive model to it. The natural explanation is that higher surfactant loading raises the effective viscosity and thickens the interfacial layers, slowing gravitational thinning.

Analysis

Speeds were reduced from the video timings and analysed in Python. Regressions are fitted to the individual runs rather than to the average of each condition. Averaging the three runs first would smooth away how much they disagree with each other, so the fitted line would look tighter than the measurements actually support.

import numpy as np

# Frame indices of puncture and of the last film collapsing, read off the
# video, one row per trial. Speed is the extension crossed per unit time.
FPS = 30.0

def speed(extension_cm, f_start, f_end, fps=FPS):
    """Front speed in cm/s for one run."""
    dt = (f_end - f_start) / fps
    return extension_cm / dt, dt

trials = [  # extension (cm), puncture frame, final-collapse frame
    (30, 112, 195),
    (30, 128, 222),
    (30, 141, 232),
]

v = [speed(L, a, b) for L, a, b in trials]
for (s, dt), (L, _, _) in zip(v, trials):
    print(f"{L} cm: {dt:5.2f} s -> {s:5.1f} cm/s")

mean = np.mean([s for s, _ in v])
sd   = np.std([s for s, _ in v], ddof=1)
print(f"mean {mean:.1f} +/- {sd:.2f} cm/s over n = {len(v)}")

import numpy as np
from scipy import stats

# Soap-to-water volume ratio -> front speed, three trials per condition.
CONC = np.array([0.10, 0.12, 0.14, 0.16, 0.18, 0.20])
trials = {  # 30 cm extension
    0.10: (13.8, 14.6, 13.9), 0.12: (10.5, 12.3, 12.6),
    0.14: (11.9, 10.2, 12.1), 0.16: (12.1, 11.8, 10.0),
    0.18: ( 9.7, 11.4, 11.6), 0.20: (11.5,  9.8, 11.1),
}

# Fit the trials, not the means: the means hide the replicate scatter.
x = np.repeat(CONC, 3)
y = np.concatenate([trials[c] for c in CONC])
fit = stats.linregress(x, y)
print(f"slope {fit.slope:7.1f} cm/s per unit ratio  "
      f"R2 {fit.rvalue ** 2:.3f}  p {fit.pvalue:.2e}")

# Is a straight line the right model? Compare residual sums of squares.
mu  = np.array([np.mean(trials[c]) for c in CONC])
lin  = np.polyfit(CONC, mu, 1)
quad = np.polyfit(CONC, mu, 2)
for name, coef in (("linear", lin), ("quadratic", quad)):
    r = mu - np.polyval(coef, CONC)
    print(f"{name:10s} RSS {np.sum(r ** 2):6.3f}")

# Leverage check: the lowest concentration is also the least comparable
# condition, since it ruptured before longer extensions were reachable.
m = x > 0.10
sub = stats.linregress(x[m], y[m])
print(f"without 0.10: R2 {sub.rvalue ** 2:.3f}  p {sub.pvalue:.3f}")

import numpy as np
from scipy import stats

# Front speed at matched extensions, puncture at top vs at the base.
TOP = {30: (11.5,  9.8, 11.1), 40: (14.6, 12.8, 14.3),
       50: (19.2, 17.1, 18.6), 60: (22.4, 20.5, 21.9)}
BOT = {30: (20.3, 22.1, 20.6), 40: (30.5, 32.2, 30.3),
       50: (37.8, 40.1, 37.9), 60: (48.9, 51.2, 48.7)}

for L in TOP:
    t, p = stats.ttest_ind(BOT[L], TOP[L], equal_var=False)
    ratio = np.mean(BOT[L]) / np.mean(TOP[L])
    print(f"{L} cm: {ratio:.2f}x faster  t {t:5.2f}  p {p:.5f}")

tm = np.array([np.mean(TOP[L]) for L in TOP])
bm = np.array([np.mean(BOT[L]) for L in BOT])
print(f"pooled: {tm.mean():.1f} -> {bm.mean():.1f} cm/s, "
      f"ratio {bm.mean() / tm.mean():.2f}")

# The bottom runs sat ~5 cm higher than nominal. Charge that offset
# against the measured extension slope before believing the effect.
SLOPE = 0.327          # cm/s per cm, fitted at this concentration
adj = bm - 5 * SLOPE
print(f"after debiting the 5 cm offset: ratio {adj.mean() / tm.mean():.2f}")

Future work

The next phase of this study is already scoped. The results above establish which variables matter and how large each effect is. The work ahead is to measure the two quantities Taylor–Culick actually depends on, resolve the collapse film by film, and close the loop between a predicted velocity and a measured one.

Film thickness comes first, because it enters the velocity under a square root and carries the explanation of all three main results. Thin-film interferometry along the stack will give the axial thickness profile directly, which turns the drainage gradient into a measured quantity and makes the rupture-position argument quantitative: once the profile is known, an expected velocity follows at every point along the axis. Pairing that with tensiometry across the concentration and temperature series supplies $\gamma$ at each condition, and together the two convert Taylor–Culick from a way of reading the direction of an effect into a prediction that can be checked against a number. The same tensiometry will place our concentration range relative to the critical micelle concentration, which sharpens the interpretation of the concentration result considerably.

Second, resolving the collapse film by film. An order-of-magnitude estimate from Taylor–Culick puts the recession of a single film across the 4 cm slinky radius between roughly 5 and 25 ms, against a front that crosses the whole stack in about 2.8 s, so high-speed imaging in the range of thousands of frames per second will separate the two timescales cleanly. With sheet illumination to bring out the neck of each catenoid, the frame-by-frame extraction of neck radius and film half-angle established in the published catenoid literature becomes available on a connected stack rather than a single film. That turns one interval per run into a velocity profile along the axis, and it tests the central prediction of this work directly: a bottom-initiated front should accelerate as it climbs into thinner films, and a top-initiated front should decelerate.

Alongside these measurements we will refine the protocol in four ways. Making each solution up to a fixed total volume separates concentration from the quantity of liquid in the container, so the two can be varied independently. Recording the coil count and pitch fixes the per-film geometry and lets us locate the extension at which the coil gap exceeds the limit for a catenoid between two equal rings, roughly 1.33 times the ring radius; that boundary is worth reaching deliberately, since films held beyond it are no longer minimal surfaces and should behave differently. A denser temperature sweep will resolve the flattening above roughly 40 °C and test whether it marks the crossover where the falling surface tension and the thinning film balance against each other. And reading the true extension for each bottom-rupture run will let the two rupture geometries be compared at exactly matched length, alongside a dilute extension of the concentration series that tests the Gibbs mechanism in the regime where surface tension still responds strongly to added surfactant.