bubble¶

Bases: Module

Standard bubble state.

Parameters:

Bases: Module, ABC

Abstract base for state-dependent (or constant) bubble properties.

Any callable eqx.Module that maps a BubbleState to a scalar array inherits from this class. The concrete subclass ConstantProperty handles the common case of a fixed value; state-dependent models, such as MarmottantSurfaceTension, override __call__ directly.

Bases: Property

A property that returns a constant value regardless of state.

Parameters:

Bases: Property

A Property backed by an Equinox neural network.

The network receives a normalised 1-D input [R / R0] and must return a 1-D output of shape (1,). You are responsible for any output transform that the physics needs. For example, wrap the final layer output with jax.nn.softplus or jnp.exp to keep a surface tension positive.

Parameters:

Examples:

>>> import equinox as eqx
>>> import jax
>>> import jax.numpy as jnp
>>> key = jax.random.PRNGKey(0)
>>> mlp = eqx.nn.MLP(in_size=1, out_size=1, width_size=32, depth=3,
...                  final_activation=jax.nn.softplus, key=key)
>>> sigma = NeuralProperty(net=mlp)

Bases: Module, ABC

Internal gas pressure model.

Computes the outward gas pressure \(p_\text{gas}\) as a function of the instantaneous state. The model reads the equilibrium configuration, R0 and P_gas0, directly from state, so gas models store only their intrinsic physical parameters, such as the polytropic exponent.

Concrete models include the polytropic law (PolytropicGas) and the van der Waals corrected gas (VanDerWaalsGas).

is_admissible(state: BubbleState) -> Array

58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

def is_admissible(self, state: BubbleState) -> jax.Array:
    """Return whether the gas law is defined at `state`.

    The default requires a positive radius, `R > 0`. A gas law with a
    smaller domain, such as
    [`VanDerWaalsGas`][jbubble.bubble.gas.VanDerWaalsGas], overrides
    this method. [`solve_eom`][jbubble.solver.solve_eom] uses it,
    through
    [`EquationOfMotion.is_admissible`][jbubble.bubble.eom.EquationOfMotion.is_admissible],
    to keep gradients finite when a rejected trial step leaves the
    domain.

    Parameters
    ----------
    state : BubbleState
        State to check. Its leaves may be NaN or infinite.

    Returns
    -------
    jax.Array
        Boolean scalar.
    """
    return state.R > 0

Bases: GasModel

Polytropic gas law.

Parameters:

Bases: GasModel

Hard-core corrected polytropic gas (van der Waals).

where \(h = h_\text{frac} R_0\) is the van der Waals hard-core radius: the radius of the bubble's gas content compressed to its excluded volume. The hard core is a property of the gas, not of the bubble, so h_frac depends on the gas species and its equilibrium density. The pressure diverges as \(R \to h\) and the law is undefined for \(R \le h\), which is_admissible reports to the solver.

Parameters:

is_admissible(state: BubbleState) -> Array

136
137
138
139
140
141
142
143
144
145
146
147
148
149
150

def is_admissible(self, state: BubbleState) -> jax.Array:
    """Return whether `R` lies outside the hard core, `R > h_frac * R0`.

    Parameters
    ----------
    state : BubbleState
        State to check. Its leaves may be NaN or infinite.

    Returns
    -------
    jax.Array
        Boolean scalar.
    """
    hard_core = self.h_frac(state) * state.R0
    return hard_core < state.R

Bases: Module, ABC

Bubble shell or coating model.

Computes the total inward stress from the shell, including:

• Laplace pressure from surface tension, \(2\sigma(R)/R\).
• Shell viscous dissipation, for example \(4\kappa_s\dot{R}/R^2\).
• Shell elastic restoring forces, for thick shells.

Every ShellModel holds a Property as its sigma field. The field's as_property converter turns a plain float into a Property.

Parameters:

p_laplace(state: BubbleState) -> Array

63
64
65

def p_laplace(self, state: BubbleState) -> jax.Array:
    r"""Return the Laplace pressure from surface tension, $2\sigma(R)/R$ [Pa]."""
    return 2.0 * self.sigma(state) / state.R

p_elastic(state: BubbleState) -> Array

67
68
69
70

@abc.abstractmethod
def p_elastic(self, state: BubbleState) -> jax.Array:
    """Return the elastic contribution from the shell [Pa]."""
    ...

p_viscous(state: BubbleState) -> Array

72
73
74
75

@abc.abstractmethod
def p_viscous(self, state: BubbleState) -> jax.Array:
    """Return the viscous contribution from the shell [Pa]."""
    ...

Bases: ShellModel

No shell coating, only Laplace pressure.

Suitable for uncoated gas bubbles. Accepts a plain float for sigma, for example 72e-3 for water.

Parameters:

Bases: ShellModel

Thin lipid shell with surface viscosity.

This is the shell model of Marmottant et al. (2005) and of its smoothed variants.

Parameters:

Bases: ShellModel

Incompressible viscoelastic shell of finite thickness (Church 1995).

The shell is a layer of incompressible Kelvin-Voigt solid with shear modulus \(G_s\), shear viscosity \(\mu_s\), and thickness \(d_s\) at the equilibrium radius \(R_0\). Integrating its radial stress across the layer gives the shell stresses of the Church (1995) model in the finite-deformation form of Qin and Ferrara (2010, Eq. 10):

where \(R\) is the outer, liquid-side radius and \(V_s\) is the shell volume divided by \(4\pi/3\). Because the shell is incompressible, \(V_s\) stays constant, and \(R^3 - V_s\) is the cube of the inner radius. The total shell pressure is

For a thin shell, \(d_s \ll R_0\), so \(V_s \approx 3 R_0^2 d_s\), and the model reduces to the thin-shell model of Hoff et al. (2000):

The viscous terms agree exactly in that limit, and the elastic terms agree to first order in the strain. For small oscillations about \(R_0\), the shell therefore adds a stiffness of \(12 G_s d_s / R_0\) per unit radial strain and a damping of \(12 \mu_s d_s / R_0^2\) per unit wall velocity, with relative corrections of order \(d_s/R_0\).

Parameters:

Bases: Property

Piecewise Marmottant surface tension law.

Three regimes, based on the radius \(R\) relative to the buckling radius \(R_b\) and the rupture radius \(R_r\):

where \(R_b\) is R_buckle_ratio * state.R0, \(\sigma_r\) is sigma_rupture, and continuity of \(\sigma\) at the elastic-to-ruptured transition gives \(R_r = R_b\sqrt{1 + \sigma_r/\chi}\).

The surface tension at \(R_0\) is \(\sigma_0 = \chi\left[(1/r_b)^2 - 1\right]\), where \(r_b\) is R_buckle_ratio, so \(r_b = (1 + \sigma_0/\chi)^{-1/2}\). A ratio of 1 or more starts the bubble buckled.

Parameters:

Bases: Property

Marmottant surface tension law with smoothed corners.

The piecewise Marmottant law clips the elastic surface tension to the interval \([0, \sigma_r]\). This law, which is jbubble's own construction, replaces the hard clip with a logistic-smoothed clamp whose width \(\varepsilon\), set by smoothing, is a fraction of \(\sigma_r\):

where \(r_b\) is R_buckle_ratio, \(\chi\) is chi, \(\sigma_r\) is sigma_rupture, and \(\operatorname{softplus}(x) = \ln(1 + e^x)\). Equivalently, \(\sigma\) is the Marmottant law averaged over a logistic spread of scale \(\varepsilon\sigma_r\) in the elastic surface tension, which rounds the buckling and rupture corners.

The law has the following properties:

• It's infinitely differentiable in \(R\), \(\chi\), \(r_b\), and \(\sigma_r\).
• It's non-decreasing in \(R\) and stays within \([0, \sigma_r]\).
• It differs from the Marmottant law by at most \(\varepsilon \ln 2 \, \sigma_r\) for every \(\chi\) and \(r_b\), and it reaches that bound at the two corners. With the default smoothing = 0.01 and \(\sigma_r = 0.072\) N/m, the bound is 0.5 mN/m.
• It converges uniformly to the Marmottant law as \(\varepsilon \to 0\).
• It's well-posed for all positive R_buckle_ratio, chi, and sigma_rupture, including bubbles that start buckled (R_buckle_ratio of 1 or more) or ruptured.

Use this law for gradient-based fitting of a lipid shell. Use MarmottantSurfaceTension for the exact piecewise law in forward simulations.

Parameters:

Raises:

Bases: Property

Marmottant-Gompertz surface tension law (Gümmer et al. 2021).

A Gompertz function replaces the corners of the piecewise Marmottant law with a smooth curve that takes the same inputs (Gümmer, Schenke, and Denner 2021, Eqs. 11 and 13-15):

where \(r_b\) is R_buckle_ratio, \(\chi\) is chi, and \(\sigma_r\) is sigma_rupture (the paper's \(\sigma_c\)). The coefficient \(b\) pins \(\sigma(R_0) = \sigma_0\), the Marmottant elastic value. The coefficient \(c\) makes the maximum slope of the curve, \(c\,\sigma_r/(e R_b)\) where \(\sigma = \sigma_r/e\), equal to the Marmottant slope at \(R = R_b\sqrt{1 + \sigma_r/(2\chi)}\). The paper takes \(\sigma_0\) as its input; R_buckle_ratio carries the same information, with \(r_b = (1 + \sigma_0/\chi)^{-1/2}\) (Eq. 11).

The class evaluates the algebraically identical form

and caps the inner exponent at 50, so a deep compression gives \(\sigma = 0\) with a finite gradient instead of an overflow.

Use this law to reproduce Gümmer et al. (2021). For gradient-based fitting, use SmoothMarmottantSurfaceTension, which converges to the Marmottant law. The Notes explain why.

Parameters:

Raises:

Bases: Module, ABC

Surrounding medium (fluid or tissue) model.

Computes the total inward viscous and elastic stresses that the surrounding medium exerts on the bubble wall. For a Newtonian liquid this is \(4\mu\dot{R}/R\). Viscoelastic media, such as Kelvin-Voigt and neo-Hookean media, add elastic restoring forces.

Subclasses implement separate methods for the viscous and elastic contributions, and the default __call__ sums them to get the total medium pressure \(p_\text{medium}(\text{state})\).

The mu field's as_property converter turns a plain float into a Property.

The medium also holds the liquid's density rho_L and speed of sound c_L, which the equations of motion and the emission models read. They are plain constants, not Property fields, because the equations of motion assume a liquid of constant density and sound speed. Both are keyword-only and default to water at 20 °C, so a subclass can add required fields after them.

Parameters:

p_viscous(state: BubbleState) -> Array

78
79
80
81

@abc.abstractmethod
def p_viscous(self, state: BubbleState) -> jax.Array:
    """Return the viscous contribution to the medium pressure [Pa]."""
    ...

p_elastic(state: BubbleState) -> Array

83
84
85
86

@abc.abstractmethod
def p_elastic(self, state: BubbleState) -> jax.Array:
    """Return the elastic contribution to the medium pressure [Pa]."""
    ...

Bases: MediumModel

Newtonian liquid medium.

Parameters:

Bases: MediumModel

Kelvin-Voigt viscoelastic medium (Yang & Church 2005).

The elastic term integrates a linear (Hookean) elastic stress through the incompressible surrounding solid. Yang and Church take the displacement field \(u(r) = (R^3 - R_0^3)/(3r^2)\) of the incompressible radial motion, so the deviatoric stress difference is \(\tau_{rr} - \tau_{\theta\theta} = -2G(R^3 - R_0^3)/r^3\), and

The model reads \(R_0\) from state.R0. The elastic term is zero at \(R = R_0\), linearises to \(4G(R - R_0)/R_0\), which is the same small-strain limit as NeoHookeanMedium, and stays bounded by \(4G/3\) under expansion. For large strains, prefer the neo-Hookean model, whose constitutive law holds at finite strain.

Parameters:

Bases: MediumModel

Neo-Hookean finite-strain viscoelastic medium.

Captures finite-strain behaviour at large oscillation amplitudes, where the Kelvin-Voigt linear approximation breaks down.

The elastic pressure comes from integrating the deviatoric stress difference \(\sigma_{\theta\theta} - \sigma_{rr}\) through the surrounding incompressible neo-Hookean solid. With the substitution \(v = r_0/r\), the integral collapses exactly to \(2G\int_{R_0/R}^{1} (1 + v^3)\,\mathrm{d}v\), which gives

Key behaviour:

• Zero at \(R = R_0\): no elastic stress at equilibrium.
• Reduces to \(4G(R - R_0)/R_0\) for small strains, identical to Kelvin-Voigt.
• Saturates to \(5G/2\) as \(R \to \infty\). Physically, the material near the bubble thins to zero, so the elastic pressure plateaus.
• Diverges to \(-\infty\) as \(R \to 0\): strong restoring force under compression.

The total medium pressure is

Parameters:

Bases: MediumModel

Power-law (generalised Newtonian) surrounding medium.

The viscosity depends on the local shear rate, \(\eta = K\dot{\gamma}^{\,n-1}\), where \(\dot{\gamma} = \sqrt{2\mathbf{D}:\mathbf{D}}\) is the rheometric shear rate, the convention under which rheometers report \(K\) and \(n\). For the incompressible radial flow \(u = R^2\dot{R}/r^2\) the shear rate is \(\dot{\gamma}(r) = 2\sqrt{3}\,|R^2\dot{R}|/r^3\), which is \(\dot{\gamma}_w = 2\sqrt{3}\,|\dot{R}|/R\) at the bubble wall. Integrating the viscous stress through the liquid,

The factor \(1/n\) comes from the radial variation of \(\eta\). For \(n = 1\) this is exactly the Newtonian result \(4K\dot{R}/R\).

The pure power law has an infinite viscosity at zero shear rate when \(n < 1\), which makes the right-hand side non-smooth whenever \(\dot{R}\) changes sign. The model regularises the wall shear rate smoothly:

Special cases:

• \(n < 1\): shear-thinning (biological tissue, blood, polymer gels).
• \(n = 1\): Newtonian; recovers NewtonianMedium with viscosity \(K\).
• \(n > 1\): shear-thickening.

Parameters:

Bases: Module, ABC, Generic

Macroscopic equation of motion for bubble dynamics.

Assembles a GasModel, a ShellModel, and a MediumModel into a complete ODE right-hand side. Concrete subclasses encode different formulations (Rayleigh-Plesset, Keller-Miksis, and others), which differ in how they relate the liquid-side boundary pressure \(p_L\) to the radial acceleration \(\ddot{R}\).

The type parameter StateType binds each equation of motion (EoM) to the state it operates on. The built-in EoMs use BubbleState; an EoM that carries extra degrees of freedom binds a BubbleState subclass, so the type checker catches a wrong state type before it fails at runtime deep inside a JAX trace.

The concrete helper p_L computes the liquid-side boundary pressure:

Because p_L is a regular method on an Equinox module, JAX differentiates through it automatically. EoMs that need \(\mathrm{d}p_L/\mathrm{d}t\), such as Keller-Miksis, use jax.grad(self.p_L) instead of hand-coding analytical derivatives for every combination of components.

The EoM receives the driving acoustic pressure as a callable p_ac_fn, so EoMs that need \(\mathrm{d}p_\text{ac}/\mathrm{d}t\) compute it with jax.grad.

initial_state seeds the equilibrium configuration, R0 and P_gas0, into the state. The ODE carries them as frozen constants, with zero time derivatives in the standard case.

Two hooks tell solve_eom about the state: is_admissible marks the states where the right-hand side is defined, and state_scale gives the magnitude of each state field, which sets the meaning of the solver tolerances. Override them in a subclass that adds state fields.

Parameters:

p_L(state: BubbleState) -> Array

105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

def p_L(self, state: BubbleState) -> jax.Array:
    r"""Liquid-side boundary pressure.

    $$
    p_L = p_\text{gas}(\text{state}) - p_\text{shell}(\text{state})
        - p_\text{medium}(\text{state})
    $$

    Parameters
    ----------
    state : BubbleState
        Current bubble state.

    Returns
    -------
    jax.Array
        Boundary pressure [Pa].
    """
    return self.gas(state) - self.shell(state) - self.medium(state)

initial_state(
    *, R: ArrayLike | None = None, R_dot: ArrayLike | None = None
) -> BubbleState

125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167

def initial_state(
    self,
    *,
    R: ArrayLike | None = None,
    R_dot: ArrayLike | None = None,
) -> BubbleState:
    r"""Return an initial state seeded with the equilibrium configuration.

    By default the bubble starts at rest at its equilibrium radius. Seeds
    `R0` and the Laplace-equilibrium gas pressure into the state:

    $$
    P_{\text{gas},0} = P_\text{amb} + \frac{2\sigma(R_0)}{R_0}
    $$

    Pass `R` or `R_dot` to start away from equilibrium, for example
    `eom.initial_state(R=1.2 * eom.R0)` for a bubble released from rest
    at 1.2 times its equilibrium radius. `R0` and `P_gas0` stay at their
    equilibrium values, so the gas pressure at `R` follows the gas law.

    Override this method for coupled systems with a larger state vector.
    An override must accept the same keyword arguments.

    Parameters
    ----------
    R : float or jax.Array, optional
        Initial radius [m]. `None` uses `R0`.
    R_dot : float or jax.Array, optional
        Initial wall velocity [m/s]. `None` uses `0`.

    Returns
    -------
    BubbleState
        State with `R`, `R_dot`, `R0`, and `P_gas0` set.
    """
    R0 = jnp.asarray(self.R0)
    P_gas0 = (
        jnp.asarray(self.P_amb)
        + 2.0 * self.shell.sigma(BubbleState(R=R0, R0=R0)) / R0
    )
    R_init = R0 if R is None else jnp.asarray(R)
    R_dot_init = jnp.zeros_like(R0) if R_dot is None else jnp.asarray(R_dot)
    return BubbleState(R=R_init, R_dot=R_dot_init, R0=R0, P_gas0=P_gas0)

is_admissible(state: StateType) -> Array

169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193

def is_admissible(self, state: StateType) -> jax.Array:
    """Return whether the right-hand side is defined at `state`.

    [`solve_eom`][jbubble.solver.solve_eom] never evaluates the equation
    of motion at an inadmissible state, so its gradients stay finite
    when an adaptive solver tries, and then rejects, a trial step that
    leaves the physical domain. The default requires `R > 0` and defers
    to [`GasModel.is_admissible`][jbubble.bubble.gas.GasModel.is_admissible],
    which also excludes, for example, the van der Waals hard core.

    Override this method in a subclass whose state has other fields with
    a restricted domain, and combine the result with
    `super().is_admissible(state)`.

    Parameters
    ----------
    state : StateType
        State to check. Its leaves may be NaN or infinite.

    Returns
    -------
    jax.Array
        Boolean scalar.
    """
    return (state.R > 0) & self.gas.is_admissible(state)

state_scale(state: StateType) -> StateType

195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247

def state_scale(self, state: StateType) -> StateType:
    r"""Return the characteristic magnitude of each state field.

    [`solve_eom`][jbubble.solver.solve_eom] integrates the dimensionless
    state `state / state_scale(state)`, so the step-size controller's
    `rtol` and `atol` mean the same thing for a 50 nm bubble as for a
    50 µm one. The default scales are:

    - `R` and `R0`: the equilibrium radius `state.R0` [m].
    - `R_dot`: the velocity scale $\sqrt{P_\text{amb}/\rho_L}$
      [m/s], with $\rho_L$ from `medium.rho_L`, about 10 m/s in water
      at atmospheric pressure. It's the
      natural velocity of bubble dynamics: $R_0\omega_0 \approx
      \sqrt{3\kappa P_\text{amb}/\rho_L}$ for the Minnaert frequency
      $\omega_0$.
    - `P_gas0`: its own value `state.P_gas0` [Pa], or `P_amb` if that
      value isn't positive.
    - Any other field of a `BubbleState` subclass: `1`, that is, SI
      units.

    Because `R0` and `P_gas0` are scaled by their own values, the solver
    sees them exactly: `(x / x) * x == x` in floating point.

    Override this method to give extra state fields a scale. The scale
    must be positive and finite; the solver treats it as a constant.

    Parameters
    ----------
    state : StateType
        Initial state, with a positive `R0`.

    Returns
    -------
    StateType
        Scale for each field, with the same structure as `state`.
    """
    P_amb = jnp.asarray(self.P_amb)
    velocity = jnp.sqrt(P_amb / jnp.asarray(self.medium.rho_L))
    R0 = jnp.asarray(state.R0)
    P_gas0 = jnp.asarray(state.P_gas0)
    pressure = jnp.where(P_gas0 > 0, P_gas0, P_amb)
    ones = jax.tree_util.tree_map(jnp.ones_like, state)
    values = (R0, velocity, R0, pressure)
    return eqx.tree_at(
        lambda s: (s.R, s.R_dot, s.R0, s.P_gas0),
        ones,
        tuple(
            jnp.broadcast_to(v, jnp.shape(leaf)).astype(jnp.result_type(leaf))
            for v, leaf in zip(
                values, (state.R, state.R_dot, state.R0, state.P_gas0), strict=True
            )
        ),
    )

Bases: EquationOfMotion

Rayleigh-Plesset equation of motion (incompressible liquid).

The simplest bubble dynamics EoM, which assumes an incompressible surrounding liquid. It takes only the parameters of EquationOfMotion, and reads \(\rho_L\) from medium.rho_L.

Bases: EquationOfMotion

Modified Rayleigh-Plesset equation with gas radiation damping.

Adds a first-order compressibility correction to the gas pressure term only, as used by Marmottant et al. (2005):

where autodiff computes \(\mathrm{d}p_\text{gas}/\mathrm{d}t = (\partial p_\text{gas}/\partial R)\,\dot{R}\).

It takes only the parameters of EquationOfMotion, and reads \(\rho_L\) and \(c_L\) from medium.rho_L and medium.c_L.

Bases: EquationOfMotion

Keller-Miksis equation of motion (first-order compressibility).

Accounts for liquid compressibility up to first order in the Mach number \(M = \dot{R}/c_L\):

JAX autodiff computes the time derivative \(\mathrm{d}p_L/\mathrm{d}t\) automatically, through the chain rule on p_L:

The \(\ddot{R}\) term moves to the left-hand side, so this EoM works with any combination of gas, shell, and medium models without hand-coded derivatives.

It takes only the parameters of EquationOfMotion, and reads \(\rho_L\) and \(c_L\) from medium.rho_L and medium.c_L.