When a Microbubble Burns Like a Star — Keller–Miksis and the Secondary Bjerknes Force#

A 5 μm air bubble meets an ultrasound wave, swells, and is then crushed to a tenth of its size. In that collapse the gas inside reaches thousands of degrees. Sonoluminescence, the star-like flash in a cup of water, comes from that short interval. An ultrasonic cleaner and a lithotripter use the same collapse.

Keller–Miksis is the equation that pushes the radius. Following Nagy and Hegedűs (2025), one bubble is integrated, then the secondary Bjerknes force between two bubbles.

“

Paper: D. Nagy, F. Hegedűs, "Assessing the accuracy of the coupled-spherical-bubble approach for bubble pairs in an acoustic field", Ultrasonics Sonochemistry 123 (2025) 107651. DOI: 10.1016/j.ultsonch.2025.107651

Push the amplitude slider toward 1.2 bar and the radius swells once, hard; the model toggle only changes the height of the rebound after collapse.

The Rayleigh–Plesset rebound is the taller one. On Keller–Miksis the gap shrinks. The liquid sound speed written in the code is 1481 m/s.

The Equation Governing a Single Bubble#

Treat the bubble as a sphere and the only unknown is the radius $R(t)$. Integrate momentum in an incompressible liquid and Rayleigh–Plesset falls out.

$\rho_L$ is the liquid density. Dots are time derivatives. The left side is the inertia of liquid pushed out radially. The right side is the pressure difference on the wall.

Liquid pressure at the wall, $p_L$, has three terms.

$p_G$ is the gas pressure inside, $\sigma$ surface tension, $\mu_L$ liquid viscosity. The second term pulls the bubble shut. The third resists the wall motion.

A polytropic closure shuts the gas.

$R_E$ is the equilibrium radius, $p_0$ the ambient pressure, $n_G$ the exponent. The code uses 1.4. Halve the radius and the factor is $2^{3\times 1.4}=2^{4.2}$, about 18. That rebound throws the collapse back out.

Driving rides on the far-field pressure.

With the code defaults $p_A=1.2\times 10^5$ Pa and $p_0=1.0\times 10^5$ Pa, the floor of $p_\infty$ is $p_0-p_A=-2\times 10^4$ Pa. In the half-cycle where absolute pressure goes negative, a polytropic air law does not stand in for a vapor cavity. The cavitation number defined from vapor pressure is the scale for that interval. The bubble swells while pressure is low and is crushed while pressure is high.

The Liquid Compresses Too — the Keller–Miksis Correction#

The weak spot in Rayleigh–Plesset is the incompressible liquid. When the wall speed at collapse approaches the speed of sound, that assumption breaks. Keller–Miksis keeps a finite wall sound speed $c_L$ through first order.

When $\dot{R}/c_L$ is small the parentheses collapse to 1 and Rayleigh–Plesset returns. As the wall speed grows, those terms radiate collapse energy as sound and hold the amplitude down. Nagy and Hegedűs grade the models by wall Mach number. Rayleigh–Plesset is order zero, Keller–Miksis order one, Gilmore order two.

The rebound gap on the toggle above is that sound. $c_L=1481$ m/s is the water sound speed in the code.

Integrating the Collapse Directly#

Drop the state to $[R,\dot{R}]$ and integrate with fourth-order Runge–Kutta. The collapse is sharp enough that the step is a nanosecond.

import numpy as np
 
P0, RHO, SIGMA, MU, C, KAPPA = 1.0e5, 998.0, 0.0725, 1.0e-3, 1481.0, 1.4
 
def bubble_rhs(y, t, RE, pA, f):
    R, V = y
    w = 2 * np.pi * f
    pG0 = P0 + 2 * SIGMA / RE
    pG = pG0 * (RE / R) ** (3 * KAPPA)
    pL = pG - 2 * SIGMA / R - 4 * MU * V / R
    pInf = P0 - pA * np.sin(w * t)
    # d/dt(pL - pInf), with the R-ddot term isolated to the left side
    dP = (-3 * KAPPA * pG * V / R + 2 * SIGMA * V / R**2
          + 4 * MU * V**2 / R**2 + pA * w * np.cos(w * t))
    num = (-1.5 * (1 - V / (3 * C)) * V**2
           + (1 + V / C) * (pL - pInf) / RHO
           + R * dP / (RHO * C))
    den = (1 - V / C) * R + 4 * MU / (RHO * C)
    return np.array([V, num / den])
 
def integrate_km(RE=5e-6, pA=1.2e5, f=60e3, dt=5e-10, steps=200000):
    y = np.array([RE, 0.0])
    t, Rlog = 0.0, []
    for _ in range(steps):
        k1 = bubble_rhs(y, t, RE, pA, f)
        k2 = bubble_rhs(y + 0.5 * dt * k1, t + 0.5 * dt, RE, pA, f)
        k3 = bubble_rhs(y + 0.5 * dt * k2, t + 0.5 * dt, RE, pA, f)
        k4 = bubble_rhs(y + dt * k3, t + dt, RE, pA, f)
        y = y + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
        y[0] = max(y[0], 0.1 * RE)   # guard against singular collapse
        t += dt
        Rlog.append(y[0])
    return np.array(Rlog)
 
R = integrate_km()
print(f"R_max / R_E = {R.max() / 5e-6:.2f}")
print(f"R_min / R_E = {R.min() / 5e-6:.2f}")
# R_max / R_E = 3.71
# R_min / R_E = 0.10

The comments record $R_{\max}/R_E=3.71$ and $R_{\min}/R_E=0.10$. The bubble reaches 3.7 equilibrium radii, then a tenth. At maximum radius the factor runs the other way. $3.71^{4.2}$ is about 250, so the gas pressure falls to about 1/250 of $p_{G0}$, near 500 Pa. The bubble is almost empty, and the collapse is the outside pressure falling into that empty volume.

Put the code constants into that 0.10 and the pressure closes. With $\sigma=0.0725$ and $R_E=5\times 10^{-6}$ m, $2\sigma/R_E=2.9\times 10^4$ Pa. Then $p_{G0}=p_0+2\sigma/R_E=1.29\times 10^5$ Pa, about 1.29 atmospheres. At minimum radius the volume ratio is $(R_E/R)^{3 n_G}=10^{4.2}\approx 1.6\times 10^4$. Gas pressure is $1.29\times 10^5\times 1.6\times 10^4\approx 2\times 10^9$ Pa, about 20,000 atm. Surface tension at that same radius, $2\sigma/R$, is ten times larger than at equilibrium, $2.9\times 10^5$ Pa, and still a trace next to 20,000 atm. The rebound is the gas, not the film.

The temperature ratio on the same exponent is $(R_E/R)^{3(n_G-1)}=10^{1.2}\approx 16$. "Thousands of degrees" matches a room-temperature start times that factor. "A few thousand atmospheres" does not match the pressure formula. The formula sits near 20,000 atm.

The integration window is written in the call. $\Delta t=5\times 10^{-10}$ s and 200,000 steps cover $10^{-4}$ s. At 60 kHz the period is $1/60000$ s, so the window is exactly six cycles.

Bubbles Push and Pull on Each Other#

The paper is about a pair. A vibrating bubble radiates a pressure wave. Under an incompressible assumption the pressure from bubble $j$ at distance $r$ scales with the volume acceleration.

That pressure adds to the neighbor's $p_\infty$ and couples the two radius equations. The time average is the secondary Bjerknes force.

$V_i$ is volume, $D$ the center distance. In phase, they attract. Out of phase, they repel. Phase follows the drive frequency against each bubble's Minnaert frequency.

With the code values $n_G=1.4$, $p_0=10^5$ Pa, and $\rho_L=998$, $\sqrt{3 n_G p_0/\rho_L}\approx 20.5$. At $R_E=5$ μm, $f_0\approx 6.53\times 10^5$ Hz, about 650 kHz. A 60 kHz drive is about an eleventh of that resonance. Opposite phase opens when the other bubble's $f_0$ lies below 60 kHz, which means a radius above about 54 μm. Put 60 kHz between a 5 μm bubble and that larger one, and the small bubble is below resonance while the large one is above it.

Spread the two radius sliders and place the drive between their Minnaert frequencies, and the arrow flips from attraction to repulsion.

Similar radii stay in phase at any frequency here. Green is attraction, red is repulsion. Only a partner above about 54 μm lets 60 kHz fall between the two resonances.

Where the Spherical Model Breaks Down#

Keller–Miksis and Gilmore both assume a sphere. Nagy and Hegedűs set that assumption against direct numerical simulation with ALPACA.

For an isolated bubble in a mild collapse the spherical model is surprisingly close. The collapse pressure still tracks even when the bubble is not a perfect sphere.

When the wall Mach number nears 1, Gilmore stays closer to the DNS than Keller–Miksis does. It pushes liquid compressibility through an equation of state.

When a pair collapses hard at short range, a jet forms, liquid punching through one side. The sphere breaks, and a spherical model overstates the internal pressure. That case belongs to DNS.

The spherical model is cheap, and it is usually enough for a cloud's void fraction or a sonoreactor. When a bubble hits a wall or a neighbor and jets, that limit comes before the number.

A jet, and the 20,000 atm was the sphere's number#

The 20,000 atm from $R_{\min}/R_E=0.10$ and $n_G=1.4$ is the polytropic value while the bubble stays a sphere. Once a jet punches the wall, that sphere is already gone. The pressure is an overestimate, and a close pair's collapse is a DNS problem.

On the same 5 μm bubble, $2\sigma/R_E=2.9\times 10^4$ Pa, 0.29 of an atmosphere. Even while the acoustic correction holds, it is still open whether this film pressure tightens an explicit step harder than Keller–Miksis does. The time window opened by implicit surface tension is where that question sits.