This repository contains the numerical code used to produce the figures in:

• PRL Letter: Spontaneous ratchet currents and transition dynamics in active wetting, N. Grodzinski, R. L. Jack, M. E. Cates (April 2026) (doi.org/10.1103/sq25-hfsb).
• PRE Companion paper: Hydrodynamic theory of wetting by active particles, N. Grodzinski, M. E. Cates , R. L. Jack (April 2026) (doi.org/10.1103/rdc2-gsqc).

If you use this code or data in academic work, please cite the above papers.

It comprises (i) a finite-volume solver for the kinetic equation of active Brownian particles (ABPs) in 1D with an external potential, and (ii) Jupyter notebooks that load saved data and reproduce each figure.

This code and data are licensed under the Creative Commons Attribution 4.0 International License (CC-BY-4.0). You are free to share and adapt the material for any purpose, including commercially, provided you give appropriate credit. See LICENSE or https://creativecommons.org/licenses/by/4.0/ for the full terms.

We solve the kinetic equation for the ABP one-particle density f(x, theta, t),

d_t f + d_x F_x + d_theta F_theta = 0,

with currents derived from the closure of the BBGKY hierarchy at second order in density (see manuscript for the full derivation). The relevant control parameters are the mean density phi, the Peclet number Pe = v_0 / sqrt(D_E D_O), and the amplitude epsilon of a localised wall potential V(x) of half-width d.

.
|-- full_dynamics.py          # FV solver (importable module)
|-- README.md
|-- data/                     # saved simulation outputs (see below)
|-- Letter_figures_1.ipynb    # PRL Fig. 1 (phase diagram, snapshots, hysteresis)
|-- Letter_figures_2.ipynb    # PRL Fig. 2 (anomalous current)
|-- Letter_figures_3.ipynb    # PRL Fig. 3 (dynamical pathways)
|-- Letter_figures_4.ipynb    # PRL Fig. 4 (minimal model)
|-- Figure2_BinodalSpinodal.ipynb     # PRE Fig. 3 (linear stability + binodals)
|-- Figure4_SurfacePhaseDiagram.ipynb # PRE Fig. 4 (surface phase diagram from sweeps)
|-- Figure5_CriticalWetting.ipynb     # PRE Fig. 5 (diverging film width, current, linear instability)
|-- Figure6_DensityCurrent.ipynb      # PRE Fig. 6 (J ~ L^{-1} scaling, J vs epsilon)
|-- Figure7_CahnCondtn.ipynb          # PRE Fig. 9 (Cahn condition cartoon)
|-- AdditionalFig_DynamicalPathway.ipynb    # Additional Figure for interest (FW <-> PW transition, bubble bursting)
|-- AdditionalFig_MinimalModel.ipynb        # Additional Figure for interest (minimal phi^4-like model, instability theory + simulation)
|-- AdditionalFig_BulkFreeEnergy.ipynb     # Additional Figure for interest (bulk free energy structure in minimal model)
|-- AdditionalFig_MinimalSolutions.ipynb   # Additional Figure for interest (minimal model PW/FW solutions)

The single solver module. The class FullSimulation wraps the grid, the adaptive time stepper, and history recording; the inner loops are JIT- compiled with Numba for speed. Key entry points:

from full_dynamics import FullSimulation

sim = FullSimulation(N_x=300, N_theta=15, v_0=10, L_x=10,
                     epsilon=2.0, hump_width=0.5)
sim.set_random(phi=0.659, delta=0.001)   # uniform + small noise IC
sim.evolve(t_f=200.0, n_records=200)
sim.plot_history()                       # rho(x, t) kymograph

After evolve(), sim.history_rho, sim.history_p1, sim.history_q1, sim.history_f and sim.history_Fx hold the time series. The numerical scheme is upwind finite-volume in x, central in theta, with adaptive Euler stepping bounded by the CFL condition. Boundary conditions in x default to periodic; closed (no-flux) BCs are implemented by subclassing and overriding update() (see GrCE_fullsim in Figure3_BinodalSpinodal.ipynb).

The closure functions d_s(rho), s(rho), Q(rho) use a polynomial approximation for the self-diffusion function $d_s(\rho)$.

Each notebook is self-contained at the level of plotting (loads from data/, generates the figure, saves nothing). The simulation cells that generated the data are kept in the notebooks but commented out (wrapped in triple-quoted strings or '''...'''). To re-run a sweep, uncomment the relevant cell; expect runtimes from minutes to hours depending on the sweep.

• Letter_figures_1.ipynb -- Fig. 1: (a) surface phase diagram in (Pe, epsilon) loaded from Letter_Fig1_a_data.npz, with the FW / PW / homogeneous regions; (b) representative snapshots in the (x, theta) plane from data/heatmap_data.npz; (c) hysteresis loops and asymmetry/thin-film-width data in epsilon at fixed Pe from data/hyseresis_curves.npz and data/profiles_hysteresis2.npy.
• Letter_figures_2.ipynb -- Fig. 2: anomalous current J ~ L^{-1} scaling (data/J_Pe10.npz) and example PW state profile (data/example_PW_state.npz).
• Letter_figures_3.ipynb -- Fig. 3: dynamical pathways between FW and PW (data/full_to_partial_wetting.npz, data/partial_to_full_wetting.npz), showing kymographs, current overlays, and asymmetry.
• Letter_figures_4.ipynb -- Fig. 4: minimal scalar field model example PW/FW evolutions (data/PW_history.npy, data/FW_history.npy).

• Figure2_BinodalSpinodal.ipynb -- spinodal from linear stability analysis (theory + numerical check via the largest eigenvalue of the linearised operator W) and binodal from full nonlinear simulations (data/binodals_data.npz, data/linear_instability.npz, data/stability.npy). Defines GrCE_fullsim (closed BCs for binodal measurement).
• Figure4_SurfacePhaseDiagram.ipynb -- surface phase diagram from a 2D sweep over (Pe, epsilon) loaded from data/combined_sweep_data.npz. Defines the asymmetry and concentration order parameters used to classify states (FW / PW / homogeneous), and produces the wetting-line cartoons.
• Figure5_CriticalWetting.ipynb -- diverging film width and the associated transverse current at the critical wetting transition; uses data/PT_current.npz and data/profiles_hysteresis2.npy.
• Figure6_DensityCurrent.ipynb -- detailed look at the PW state: density and polarisation profiles, anomalous current scaling J ~ L^{-1} (data/J_Pe10.npz, data/example_PW_state.npz).
• Figure7_CahnCondtn.ipynb -- analytical/cartoon plots illustrating the Cahn wetting condition (no simulation data -- pure plotting from closed-form expressions)
• AdditionalFig_DynamicalPathway.ipynb -- transient evolution between FW and PW including a "droplet bursting" event; uses data/full_to_partial_wetting.npz and data/partial_to_full_wetting.npz.
• AdditionalFig_MinimalModel.ipynb -- the reduced 1D scalar field model (cubic potential + current J_0). Defines run_sim for the minimal dynamics, theory line theory_Jcrit(alpha, kappa, L) for the critical current, and produces the FW/PW evolution kymographs..
• AdditionalFig_BulkFreeEnergy.ipynb -- bulk free energy structure: counts the number of solutions of the chemical-potential equation as (mu, b) vary, identifies the critical Pe for binodal coexistence (Pe_min = sqrt(2 b_min)), and plots f(rho) for a few Pe.
• AdditionalFig_MinimalSolutions.ipynb -- stationary PW / FW profiles of the minimal model and their first/second derivatives, used to check the generalised Young-Dupre construction.

The data/ directory contains the saved outputs needed to reproduce all figures without rerunning sweeps. The naming is hopefully self-explanatory from the notebook references above; in particular:

• combined_sweep_data.npz: full (Pe, epsilon) sweep used to build the surface phase diagram; contains Pe, epsilon, rho_history.
• heatmap_data.npz: f(x, theta) snapshots for the three example states (FW, PW, homogeneous).
• hyseresis_curves.npz, profiles_hysteresis2.npy: forward/backward sweeps in epsilon at fixed Pe.
• J_Pe10.npz, example_PW_state.npz: anomalous current scaling data.
• full_to_partial_wetting.npz, partial_to_full_wetting.npz: long-time evolutions through the dynamical pathway.
• PW_history.npy, FW_history.npy, AS_series.npz, final_asymmetries_fig5b.npy: minimal-model data.
• binodals_data.npz, linear_instability.npz, stability.npy: bulk spinodal and binodal data.
• PT_current.npz: critical wetting current.
• Letter_Fig1_a_data.npz: classification data.

Note: Letter_Fig1_a_data.npz (in the repo root) contains the classification array and critical epsilon values used directly in PRL Fig. 1(a); it is derived from combined_sweep_data.npz via the cells in Figure4_SurfacePhaseDiagram.ipynb.

• 2D simulations / methods which are mentioned in the manuscript (as these are not central).
• Auxiliary scripts that were used during exploration but did not feed into the published figures.

• Python >= 3.9
• numpy, scipy, matplotlib, numba, tqdm
• LaTeX installation (for text.usetex=True in figure rendering); set mpl.rcParams['text.usetex'] = False if you do not have one.

A minimal environment.yml and requirements.txt would be a sensible addition (see suggestions below).

jupyter lab Letter_figures_1.ipynb

and run all cells. All notebooks execute in under a minute when loading from data/; full re-simulation (uncommenting the generation cells) is much slower. Some figures have been cosmetically edited/assembled in InkScape, but the data/plots is the same.

If you use this code, please cite the PRL Letter and PRE paper above. BibTeX files are attached.

Noah Grodzinski, DAMTP, University of Cambridge -- [email protected].