1. Introduction
Successive drops of coloured inks mixed with surfactant create intricate patterns by Marangoni spreading in the Japanese art technique of Suminagashi (see figure 1a). A surfactant–ink drop is gently deposited at the surface of a thin layer of water, which may have a small initial concentration of endogenous surfactant due to normal environmental contamination. It then spreads outwards and equilibrates before reaching the edges of the container. Then successive drops deposited at different locations of the liquid surface form the intricate patterns. During pattern creation, the artist can blow on the surface with a straw after drop equilibrations to further deform the pattern. Eventually, the pattern is captured on pieces of paper placed onto the surface (Chambers Reference Chambers1991). Rouwet & Iorio (Reference Rouwet and Iorio2017) noticed similar patterns occurring in volcanic crater lakes, and hypothesised that similar physics were responsible: thermal gradients in the lake create Marangoni flows, and wind action creates a blowing effect, resulting in marbling patterns of the adsorbed multicoloured sediments. In this study, we seek to understand the Lagrangian trajectories of material points and curves on a surface during the spreading of surfactants on a confined surface, and thus the dynamics of any adsorbed passive tracer, similar to the advection of ink by surfactant in the Suminagashi technique.
Figure 1. (a) Pictures showing the Japanese art technique of Suminagashi (Rouwet & Iorio Reference Rouwet and Iorio2017). Successive drops of a mixture of coloured ink and surfactant are deposited on the surface of a thin film of water to create a multicoloured pattern. Blowing on the surface then creates further intricate patterning. (b) Picture of a Suminagashi pattern of ink on water created by artist Bea Mahan (Reference Mahan2011). (c) Schematic of the model problem. Circular deposits of insoluble exogenous surfactant (red) spread on the surface 
$\varOmega$ of a thin layer of viscous liquid (blue) of mean height 
$h$ confined in a rectangular region of dimensions 
$L_1$ and 
$L_2$, where the surface contains an initially uniform endogenous surfactant (green). We assume that the ratio of vertical to horizontal length scales is small enough, and that the Bond number (ratio of gravitational to surface tension forces) is large enough, for height deflections caused by spreading to be negligible, confining spreading to the flat plane of the surface 
$\varOmega$.
In addition to the cultural importance of Suminagashi, which has been part of Japanese art since the 12th century, and similar practices in China for even longer (Ishii & Muro Reference Ishii and Muro1989), understanding Marangoni-driven surface motion can help us to better understand various industrial and biological applications involving surfactants carrying passive, adsorbed material. For example, Deng et al. (Reference Deng, Zheng, Bai, Wang, Zhao and Huang2018) showed how small amounts of surfactant added to perovskite (a calcium titanium oxide mineral) can suppress the formation of islands during the drying phase of blade coating by creating Marangoni flows that keeps the solution coating even. Many other coating processes involve Marangoni flows induced by trace amounts of surfactants. Some methods of drug delivery in lungs mix pharmaceutical substances with exogenous surfactant (Haitsma, Lachmann & Lachmann Reference Haitsma, Lachmann and Lachmann2001), so that the surfactant acts as a carrier to spread the drug through the airways. In particular, surfactant replacement therapies have been used successfully in lungs of neonates affected with respiratory distress syndrome (Avery & Mead Reference Avery and Mead1959; Jobe Reference Jobe1993; Rodriguez Reference Rodriguez2003; Halliday Reference Halliday2008). The surfactant-driven spreading in the complex and confined tree-like geometry of the lungs acts against its natural endogenous surfactant (Espinosa et al. Reference Espinosa, Shapiro, Fredberg and Kamm1993; Jensen, Halpern & Grotberg Reference Jensen, Halpern and Grotberg1994; Grotberg, Halpern & Jensen Reference Grotberg, Halpern and Jensen1995; Halpern, Jensen & Grotberg Reference Halpern, Jensen and Grotberg1998; Temprano-Coleto et al. Reference Temprano-Coleto, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2018; Mcnair et al. Reference Mcnair, Temprano-Coleto, Peaudecerf, Gibou, Luzzatto-Fegiz, Jensen and Landel2023). These methods of delivery can help to overcome difficulties such as poor solubility of the pharmaceuticals (Hidalgo, Cruz & Pérez-Gil Reference Hidalgo, Cruz and Pérez-Gil2015).
Molecules and substances that act as surfactants are ubiquitous in the environment. They can cause unexpected fluid flows that have confounded scientists and engineers, as described by Manikantan & Squires (Reference Manikantan and Squires2020), who discussed the ‘hidden’ variables related to surfactant dynamics in many fluid flows. The present study addresses insoluble surfactant spreading into pre-existing, endogenous surfactant on a thin film of a bounded Newtonian viscous liquid, allowing us to use lubrication theory to approximate the Stokes flow in the liquid film. Lubrication theory for insoluble surfactant-driven flows has its origins with the work of Borgas & Grotberg (Reference Borgas and Grotberg1988) who derived coupled partial differential equations (PDEs) describing the leading-order evolution of the liquid film height and surfactant concentration. The work was then extended theoretically and experimentally by Gaver & Grotberg (Reference Gaver and Grotberg1990, Reference Gaver and Grotberg1992). Thess, Spirn & Jüttner (Reference Thess, Spirn and Jüttner1997) and Jensen & Halpern (Reference Jensen and Halpern1998) showed that the coupled equations in the limit of large Bond number could be combined into a single nonlinear diffusion (or ‘porous medium’) equation governing surfactant concentration evolution as a function of space and time. The effect of gravity is to suppress deflections of the surface, removing the functional dependence of the spreading on the dynamic film height.
In this paper, we explore a link between the theory of surfactant spreading and the theory of optimal transport. This theory was initiated by Monge (Reference Monge1781), who was trying to find the optimal way to transport mounds of soil under some cost function. The theory was extended into its modern formulation by Kantorovich (Reference Kantorovich1942, Reference Kantorovich2006). Most of its current uses are found in machine learning and image analysis (Kolouri et al. Reference Kolouri, Park, Thorpe, Slepcev and Rohde2017). A powerful result, enabling significant simplification of optimal transport problems, occurs when the cost function takes a quadratic form, yielding the quadratic Monge–Kantorovich optimal transport problem (qMK). For such cost functions, solutions for the optimal map of material from initial to final location can be shown to be the gradient of a convex function that satisfies a so-called Monge–Ampère equation. A variety of approaches have been taken to find solutions of this nonlinear equation (Froese & Oberman Reference Froese and Oberman2011; Benamou, Froese & Oberman Reference Benamou, Froese and Oberman2012, Reference Benamou, Froese and Oberman2014). Otto (Reference Otto2001), building on work by Jordan, Kinderlehrer & Otto (Reference Jordan, Kinderlehrer and Otto1998) and Benamou & Brenier (Reference Benamou and Brenier2000), showed that porous medium equations (a class of equations to which the Jensen & Halpern (Reference Jensen and Halpern1998) surfactant equation that we use in this study belongs) have the variational structure of a gradient flow on a Riemannian manifold measured by the quadratic Wasserstein distance. The square of the Wasserstein distance, which is defined as the minimiser of a functional, doubles as qMK, which suggests that solutions to the surfactant-induced transport problem may be approximated by solving the Monge–Ampère equation under certain conditions. We explore these conditions in this paper, and consider whether the Monge–Ampère equation could be an efficient tool to determine equilibrium solutions to this complex confined transport problem.
The primary aim of this study is to understand the underlying physics behind surfactant-induced Marangoni dynamics in a confined environment when the surface contains an initial endogenous concentration of surfactant, which is the case for most environmental fluids. The spreading of multiple exogenous deposits is particularly considered; this was investigated experimentally and with COMSOL
$^{\circledR}$ models recently by Iasella et al. (Reference Iasella, Sharma, Garoff and Tilton2024), showing how adjacent droplets interact and deform. A Lagrangian framework, which has been adopted in the analysis of other transport problems with nonlinear diffusive character (Meĭrmanov, Pukhnachev & Shmarev Reference Meĭrmanov, Pukhnachev, Shmarev, Kegel, Maslov, Neumann and Wells1997), enables us to compute efficiently individual surface particle trajectories and equilibrium states as functions of initial distributions. Moreover, the Lagrangian framework reveals underpinning flow phenomena such as stretching, compression and rotational motion that govern the particle trajectories. While there have been limited investigations of Lagrangian surfactant dynamics in one spatial dimension (Grotberg et al. Reference Grotberg, Halpern and Jensen1995), there is none (to our knowledge) in higher dimensions, despite the potential relevance to a variety of applications. Furthermore, while some authors have exploited the gradient flow structure of thin-film evolution equations (Thiele, Archer & Pismen Reference Thiele, Archer and Pismen2016; Henkel, Snoeijer & Thiele Reference Henkel, Snoeijer and Thiele2021), we are not aware of prior studies linking thin-film flows to optimal transport. We show how we can exploit this link for practical purposes. In particular, we describe a procedure to reproduce the intricate patterns of Suminagashi art, through resolution of the Monge–Ampère equation associated with the surfactant transport model. These results appear to capture, at least qualitatively, the dominant physics behind Suminagashi art, suggesting a powerful tool for other applications where surface transport is dominated by surfactants in confined environments.
In § 2.1, we use a two-dimensional extension of the model of Jensen & Halpern (Reference Jensen and Halpern1998) (derived in Appendix A) to describe transport of material particles on a surface. We outline a physical problem in Eulerian coordinates in a confined rectangular domain, implementing initial conditions that represent multiple deposits of exogenous surfactant spreading on a surface with an initially uniform endogenous surfactant concentration. We solve the particle-tracking problem using a finite-difference method by first solving for the evolution of surfactant concentration, and then interpolating the gradient of this solution onto a second Lagrangian grid where we integrate the surface velocity to find the trajectories of surface particles initially located at each grid point. In § 2.2, we reformulate the problem in Lagrangian coordinates, and show how it can be reduced from three to two scalar PDEs, enabling the same calculation without the intermediate step of finding the evolution of the surfactant concentration, and without the need to interpolate concentration gradients from an Eulerian to a Lagrangian grid. We solve the resulting scheme using a finite-element method. In § 2.3, we show how to approximate the equilibrium locations of surface particles as a function of their initial locations via a Monge–Ampère equation, without having to compute their intermediate trajectories. In § 3.1, we show consistency between the Eulerian and Lagrangian methods, and describe dynamical phenomena not normally associated with spreading surfactants, such as drift and flow reversals due to confinement. In § 3.2, we show that solutions of the Monge–Ampère equation approximate the equilibrium solution well when the endogenous and exogenous concentrations are of comparable magnitude, and also provide a credible approximation when the endogenous concentration is much smaller. We show how, in the limit of small endogenous concentration, the boundaries of the deposits become almost polygonal with self-similar structures at the corners, resembling a two-dimensional foam. We discuss subtle discrepancies between the Monge–Ampère solution, and a solution computed with the Eulerian particle-tracking method, indicating that surfactant transport can be considered almost, but not exactly, optimal. We analyse the two-dimensional mapping between the initial surfactant distribution and its equilibrium distribution, and discuss how the divergence and curl of the mapping can reveal regions of stretching, compression and rotational motion. Finally, we show that successive solutions of the Monge–Ampère equation, combined with divergence-free maps to mimic blowing, can be used to create a computational Suminagashi marbling pattern, illustrating the power of the optimal-transport approximation. Additional results are shown in supplementary material available at https://doi.org/10.1017/jfm.2024.334 to provide further evidence supporting the main findings and discussion presented in this paper.
2. Model and methods
2.1. The Eulerian particle-tracking problem
2.1.1. The problem and derivation of the model
We investigate the trajectories of particles on the surface of a viscous Newtonian liquid advected by surface tension gradients caused by a non-uniform concentration profile of insoluble surfactant, which is assumed to have negligible molecular diffusivity. Concentration gradients are caused by deposits of exogenous surfactant added to a uniform concentration field of endogenous surfactant. We assume that both species of surfactant have the same material properties, which combine to create a single concentration field that has a linear relationship with surface tension. A typical length scale is found from the initial size of an exogenous deposit, which is much greater than the initial height of the film. The thickness of the film is assumed to remain approximately uniform during the spreading, as we assume that any large vertical deflections are suppressed by gravity (in a large-Bond-number limit). The spreading takes place in a closed region with rectangular horizontal cross-section 
$\varOmega$, given in non-dimensional Cartesian coordinates as 
$0\leq x\leq L_1$, 
$0 \leq y\leq L_2$ confined by impermeable walls. Surfactant concentrations are scaled by the maximum initial concentration of one of the deposits. As explained in Appendix A, the surfactant is transported from its initial profile to its final equilibrium state via the nonlinear diffusion equation, which describes the evolution of the surfactant concentration as a function of space and time:
\begin{equation} \varGamma_t = \tfrac{1}{4}\,\boldsymbol{\nabla}_{\boldsymbol{x}}\boldsymbol{\cdot} (\varGamma\, \boldsymbol{\nabla}_{\boldsymbol{x}}\varGamma), \end{equation}
where 
$\boldsymbol {\nabla }_{\boldsymbol {x}}$ is the gradient operator in the 
$\boldsymbol {x}= (x,y)$ plane of the Eulerian coordinates, and 
$\varGamma _t$ is the derivative of surfactant concentration with respect to non-dimensional time 
$t$; here, time is scaled by the ratio of liquid viscosity to maximum surface tension gradient (Appendix A). We impose a no-flux boundary condition at the periphery of the domain:
\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{x}}\varGamma \boldsymbol{\cdot} \boldsymbol{n}_b = 0 \quad \text{on } \partial \varOmega, \end{equation}
where 
$\partial \varOmega$ is the boundary of the domain 
$\varOmega$, and 
$\boldsymbol {n}_b$ is a unit normal vector to the boundary of the domain. Comparison of (2.1) with the non-dimensional surface transport equation 
$\varGamma _t+\boldsymbol {\nabla }_{\boldsymbol {x}}\boldsymbol {\cdot }(\boldsymbol {u}_s\varGamma )=0$ for a flat surface and non-diffusive surfactant shows that the surface velocity is
\begin{equation} \boldsymbol{u}_s(\boldsymbol{x},t) ={-}\tfrac{1}{4}\,\boldsymbol{\nabla}_{\boldsymbol{x}} \varGamma. \end{equation}
Since we impose no flux of surfactant at the boundaries of 
$\varOmega$, as time goes to infinity, concentration gradients vanish to reach an equilibrium or steady state, so the initial concentration profile of surfactant 
$\varGamma _0(x,y)$ spreads to a uniform state with concentration 
$\bar {\varGamma }>\delta >0$, where 
$\delta$ is the initial endogenous concentration. We do not consider the singular limit 
$\delta =0$, which is beyond the scope of this study. In that case, spreading at the edges of the deposits would continue until the edges meet a solid boundary or the edges of another deposit. The final concentration relates to the initial concentration profile by
\begin{equation} \bar{\varGamma}=\frac{\int_{\varOmega}\varGamma_0(x,y) \,\mathrm{d}\kern0.7pt x\, \mathrm{d}y}{\int_{\varOmega} \,\mathrm{d}\kern0.7pt x\, \mathrm{d}y} = \frac{1}{L_1L_2}\int_{\varOmega} \varGamma_0(x,y) \,\mathrm{d}\kern0.7pt x \,\mathrm{d}y . \end{equation}
Equation (2.1) represents a natural generalisation of the spatially one-dimensional nonlinear diffusion equation derived in Jensen & Halpern (Reference Jensen and Halpern1998), and aligns with the two-dimensional formulation of Thess et al. (Reference Thess, Spirn and Jüttner1997). In stepping from one to two dimensions, an extra degree of freedom must be considered: any surface velocity field for which 
$\boldsymbol {u}_s\varGamma$ has zero divergence will not change surface concentrations but will nevertheless transport surface particles. This is illustrated in Appendix A by considering the influence of an imposed surface stress, as might arise from external blowing on the liquid film. For a monolayer close to equilibrium, the divergence of the stress field is area-changing; this is resisted by Marangoni effects (A10). However, the curl of the stress field in this simple model generates a flow that can redistribute surfactant (i.e. surface material elements carrying either endogenous or exogenous surfactant) without inducing surface tension gradients (A11). As well as being exploited by Suminagashi artists, this feature highlights a potential degeneracy in (2.1): namely, that the energetic cost of any flow that preserves concentrations of surface material elements is not captured by the evolution equation.
We now introduce a Lagrangian coordinate system 
$(x_0,y_0,\tau )$ to complement the Eulerian system 
$(x,y,t)$. We define a mapping 
$\boldsymbol {X}=(X,Y)$ between them, such that particles starting on the interface at 
$\boldsymbol {x}_0=(x_0,y_0)\in \varOmega$ at 
$t=0$ are advected at time 
$t=\tau$ to
\begin{equation} x = X(x_0,y_0,\tau), \quad y = Y(x_0,y_0,\tau). \end{equation}Since surfactant transport is purely advective under (2.1), the mapping satisfies
\begin{equation} \frac{\partial \boldsymbol{X}(x_0,y_0,\tau)}{\partial \tau} ={-}\frac{1}{4}\,\boldsymbol{\nabla}_{\boldsymbol{x}} \varGamma(\boldsymbol{X}(x_0,y_0,\tau),\tau). \end{equation}
The mapping function 
$\boldsymbol {X}(x_0,y_0,\tau )$ from initial to current particle location is the main quantity that we seek throughout this study. The initial conditions for each simulation that we perform in this study will be of the form
\begin{equation} \varGamma_0(x_0,y_0) = \begin{cases} \delta + \mathcal{F}(x_0,y_0), & \text{in } \varOmega',\\ \delta, & \text{in } \varOmega -\varOmega', \end{cases} \end{equation}
where 
$\delta =\min {\varGamma _0(x_0,y_0)}$ for all 
$(x_0,y_0)$ in 
$\varOmega$ represents the initially uniform endogenous surfactant, and 
$\mathcal {F}$ is a function describing the initial distribution of exogenous surfactant deposited in 
$\varOmega '$, a region of 
$\varOmega$. In this study, we consider only non-overlapping depositions of exogenous surfactant that are axisymmetric about their own centre, and with a radially decreasing concentration profile. Although we have studied various initial distributions for the exogenous surfactant deposits (see the supplementary material), we focus on quadratic distributions, which we denote as
\begin{equation} \mathcal{C}_q(\boldsymbol{x}_0;\boldsymbol{x}_c,r,\varGamma_{0,c}-\delta) = \begin{cases} (\varGamma_{0,c}-\delta)\left(1- \dfrac{|\boldsymbol{x}_0-\boldsymbol{x}_c|^2}{r^2}\right), & |\boldsymbol{x}_0-\boldsymbol{x}_c|\leq r,\\ 0, & |\boldsymbol{x}_0-\boldsymbol{x}_c|> r,\end{cases} \end{equation}
which is centred at 
$\boldsymbol {x}_c=(x_c,y_c)$, where the initial concentration has a local maximum 
$\varGamma _{0,c}$, with deposit radius 
$r$. The concentration profile 
$\varGamma _0$ is continuous when added to the endogenous field, and the Euclidean distance is given by 
$|\boldsymbol {x}_0-\boldsymbol {x}_c| \equiv \sqrt {(x_0-x_c)^2+(\kern0.7pt y_0-y_c)^2}$. The subscript 
$q$ in (2.8) refers to the quadratic nature of the initial concentration profile. In Appendix F and in § S5 of the supplementary material, we consider circular concentration profiles with other functional forms.
2.1.2. Scenarios studied
We have investigated scenarios involving one, two or three distinct deposits (i.e. 
$\varOmega '$ is constituted of one, two or three disconnected regions in 
$\varOmega$). The different configurations studied for the one- and two-deposit problems are presented in the supplementary material (see table S1). These two problems are helpful to understand basic dynamical features and the impact of the relevant non-dimensional parameters, as we will discuss briefly in § 3. However, the one- and two-deposit problems miss topological features that appear only with three or more exogenous deposits, such as internal corners where the edges of the deposits meet away from the domain boundaries. As we will discuss in § 3, internal corners display self-similar features. For the sake of simplicity and to enable analytical progress, we focus mainly on the three-deposit problem for the rest of this paper. Nevertheless, we anticipate that many of the results found with the three-deposit problem will also apply to problems involving more deposits. Therefore, we devise a model problem where 
$\mathcal {F}(x_0,y_0)$ consists of three circular regions of different radii (
$r_1=1$, 
$r_2$ and 
$r_3$ in non-dimensional variables; see figure 1c) containing exogenous surfactant with quadratic concentration profiles, with differing non-dimensional maximum values 
$1$, 
$\varGamma _2$ and 
$\varGamma _3$ in the different regions (the number in the subscript corresponds to the region). Deposit 1, the smallest, is centred at 
$(x_1,y_1)$; the second largest circular deposit is centred at 
$(x_2,y_2)$; the largest is centred at 
$(x_3,y_3)$. Therefore, using our notation for circular deposits (2.8), we have
\begin{equation} \mathcal{F}(x_0,y_0) = \mathcal{C}_q(x_1,y_1,1,1-\delta) + \mathcal{C}_q(x_2,y_2,r_2,\varGamma_2-\delta) + \mathcal{C}_q(x_3,y_3,r_3,\varGamma_3-\delta). \end{equation}
For the three-deposit problem, we choose 
$r_2=2$, 
$r_3=3$, 
$\varGamma _2=1$ and 
$\varGamma _3=2$. For every problem tackled in this paper and in the supplementary material, we choose 
$L_1=13$ and 
$L_2=11$. The centres of the deposits are chosen to be 
$(x_1,y_1) = (6,2)$, 
$(x_2,y_2) = (10,5)$ and 
$(x_3,y_3)=(4,7)$ for most of the solutions presented, unless otherwise stated.
2.1.3. Numerical scheme for the Eulerian particle-tracking problem
A finite-difference approximation of (2.1) and (2.6) is calculated using two rectangular grids. The first grid is used to solve for an approximation of (2.1) subject to boundary conditions (2.2) and initial conditions (2.7) and (2.9) in an Eulerian reference frame, which is accomplished using a second-order central differencing system in space, and a first-order forward Euler method in time (choosing a sufficiently small time step to ensure stability). This is solved simultaneously with a forward Euler approximation of (2.6) for the dynamics of the particle paths on a second grid in the Lagrangian reference frame. At each time step, the concentration gradient is approximated on the Eulerian grid, and interpolated onto the Lagrangian grid at the current particle locations using a linear interpolation method, meaning that the method as a whole is first-order in space and time.
The simulation is computed from 
$t=0$ to a large time 
$t=t_f$ when the solution approximates the steady state. The value of 
$t_f$ is found by considering the analysis in Appendix B, which shows how to ensure that the map is within a small tolerance vector 
$[X_{tol},Y_{tol}]^{\rm T}$ of the steady state everywhere (we set 
$[X_{tol},Y_{tol}]^{\rm T}=[10^{-3},10^{-3}]^{\rm T}$).
2.2. The Lagrangian particle-tracking problem
2.2.1. Derivation of the Lagrangian method
Rather than solving the three scalar PDEs in (2.1) and (2.6) in an Eulerian framework, it is sufficient to solve only two PDEs by adopting a Lagrangian framework, as we now demonstrate, by calculating 
$\boldsymbol {X}(x_0,y_0,\tau )$ without the intermediate step of determining surfactant concentrations. We present a Lagrangian scheme reminiscent of that presented by Carrillo, Matthes & Wolfram (Reference Carrillo, Matthes, Wolfram, Bonito and Nochetto2021) for a general Wasserstein gradient flow. The chain rule combined with (2.5a,b) yields the material derivative 
$\partial /\partial \tau |_{x_0,y_0} = \partial /\partial t |_{x,y}+ \boldsymbol {u}_s\boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x}}$, where 
$\boldsymbol {u}_s = \boldsymbol {X}_{\tau }$, with the 
$\tau$ subscript meaning the partial derivative with respect to 
$\tau$. It is also the case that
\begin{equation} \begin{pmatrix}
\dfrac{\partial}{\partial x_0} \\ \dfrac{\partial}{\partial
y_0} \end{pmatrix} = \begin{pmatrix} X_{x_0} & Y_{x_0} \\
X_{y_0} & Y_{y_0} \end{pmatrix} \begin{pmatrix}
\dfrac{\partial}{\partial x} \\ \dfrac{\partial}{\partial
y} \end{pmatrix}, \quad \text{or} \quad
\boldsymbol{\nabla}_{\boldsymbol{x_0}} =
(\boldsymbol{\nabla}_{\boldsymbol{x_0}} \boldsymbol{X}
)^{\rm T}\,\boldsymbol{\nabla}_{\boldsymbol{x}}.
\end{equation}
We define tensor calculus operators as
\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{x_0}} \begin{pmatrix} a_1 \\ a_2 \end{pmatrix} = \begin{pmatrix} a_{1x_0} & a_{1y_0} \\ a_{2x_0} & a_{2y_0} \end{pmatrix}, \quad \boldsymbol{\nabla}_{\boldsymbol{x_0}} \boldsymbol{\cdot} \begin{pmatrix} a_1 & a_2 \\ a_3 & a_4 \end{pmatrix} = \begin{pmatrix} a_{1x_0} + a_{3y_0} & a_{2x_0} +a_{4y_0} \end{pmatrix}. \end{equation}The Jacobian of the mapping (2.5a,b),
\begin{equation} \alpha \equiv \det(\boldsymbol{\nabla}_{\boldsymbol{x_0}}\boldsymbol{X}) = X_{x_0}Y_{y_0}-X_{y_0}Y_{x_0}, \end{equation}
quantifies how area elements are deformed by the map between initial and current particle positions, such that area elements 
$\mathrm {d} A_{\boldsymbol {x_0}}$ and 
$\mathrm {d} A_{\boldsymbol {x}}$ are related by 
$\mathrm {d} A_{\boldsymbol {x}} = \alpha \,\mathrm {d} A_{\boldsymbol {x_0}}$. By conservation of mass, we can equate integrals of the surfactant concentration over the Lagrangian and Eulerian domains, respectively:
\begin{equation} \int_{\boldsymbol{X}^{{-}1}(\Delta \varOmega)} \varGamma_0(x_0,y_0) \,\mathrm{d} A_{\boldsymbol{x_0}} = \int_{\Delta \varOmega} \varGamma(\boldsymbol{X},t) \, \mathrm{d} A_{\boldsymbol{x}} , \end{equation}
where 
$\boldsymbol {X}^{-1}(\Delta \varOmega )$ is the pre-image of any subset 
$\Delta \varOmega$ of the Eulerian domain 
$\varOmega$, and there is a one-to-one mapping between the domains. Using the Jacobian of the mapping, we can change variables on the right-hand side of (2.13) to give
\begin{align} \int_{\boldsymbol{X}^{{-}1}(\Delta \varOmega)} \varGamma_0(x_0,y_0) \,\mathrm{d} A_{\boldsymbol{x_0}} = \int_{\boldsymbol{X}^{{-}1}(\Delta \varOmega)} \varGamma(\boldsymbol{X}(x_0,y_0,\tau),\tau)\, \alpha(\boldsymbol{X}(x_0,y_0,\tau),\tau) \,\mathrm{d} A_{\boldsymbol{x_0}} . \end{align}
We are now integrating over the same space with respect to the same variables, and as 
$\Delta \varOmega$ is arbitrary, the integrands must be equal, yielding
\begin{equation} \varGamma(\boldsymbol{X}(x_0,y_0,\tau),\tau)\,\alpha(\boldsymbol{X}(x_0,y_0,\tau),\tau) = \varGamma_0(x_0,y_0). \end{equation}
This is the main statement of mass conservation in 
$\varOmega$, valid for any 
$\tau \geq 0$, and is key for our analysis in this subsection and the next.
The choice of Lagrangian coordinate system is arbitrary, and in the rest of this subsection, we choose a spatially non-uniform coordinate system 
$(\xi, \eta )$. This coordinate system, non-uniform in 
$\varOmega$, also defines a geometric transformation of the domain 
$\varOmega$, which is achieved by deforming 
$\varOmega$ such that 
$(\xi, \eta )$ become regularly spaced Cartesian coordinates. We call this new domain the deformed Lagrangian domain, with coordinates 
$(\xi,\eta )$ replacing 
$(x_0,y_0)$. In § 2.3 we will revert back to 
$(x_0,y_0)$, which there will refer to regular Cartesian coordinates in an undeformed copy of the Eulerian domain such that 
$(x,y)=(x_0,y_0)$ at 
$\tau =0$. (These two domains will be referred to as the deformed and undeformed Lagrangian domains, respectively.) For now, however, we choose a coordinate system 
$(\xi,\eta )$ such that the initial surfactant concentration is uniform in the deformed domain, with 
$\varGamma _0(\xi,\eta )=1$ everywhere. This new coordinate system 
$(\xi, \eta )$ defines a geometric transformation of the rectangular domain 
$\varOmega$, such that surface areas are stretched or compressed until the concentration per unit area in the deformed system is 
$1$ everywhere. To illustrate, if a region of unit area has an initial uniform concentration of 
$0.25$ in the undeformed domain, then in the deformed domain it would have an area of 
$0.25$ and therefore an initial uniform concentration of 
$1$. In the coordinate system of the deformed domain, (2.15) becomes
\begin{equation} \alpha(\boldsymbol{X}(\xi,\eta,\tau),\tau)\,\varGamma(\boldsymbol{X}(\xi,\eta,\tau),\tau) = 1. \end{equation}
With this choice, and using (2.5a,b), it follows that 
$\boldsymbol {\nabla }_{\boldsymbol {x}} (\alpha \varGamma ) = \alpha \, \boldsymbol {\nabla }_{\boldsymbol {x}} \varGamma + \varGamma \, \boldsymbol {\nabla }_{\boldsymbol {x}}\alpha = 0$, and so
\begin{equation} \alpha\,\boldsymbol{\nabla}_{\boldsymbol{x}} \varGamma ={-}\varGamma\, \boldsymbol{\nabla}_{\boldsymbol{x}}\alpha ={-}\varGamma \left(\boldsymbol{\nabla}_{\boldsymbol{\xi}} \boldsymbol{X}\right)^{{-}{\rm T}} \boldsymbol{\nabla}_{\boldsymbol{\xi}}\alpha, \end{equation}
where 
$\alpha = \det (\boldsymbol {\nabla }_{\boldsymbol {\xi }}\boldsymbol {X})$, and 
$\boldsymbol {\nabla }_{\boldsymbol {\xi }} = [ \partial /\partial \xi, \partial /\partial \eta ]^{\rm T}$. The particle velocity (2.6) is given by 
$\boldsymbol {X}_{\tau } = - \boldsymbol {\nabla }_{\boldsymbol {x}} \varGamma /4$, so (2.12), (2.16) and (2.17) give
\begin{equation} 4\alpha^2\left(\boldsymbol{\nabla}_{\boldsymbol{\xi}} \boldsymbol{X}\right)^{\rm T} \boldsymbol{X}_{\tau} = \boldsymbol{\nabla}_{\boldsymbol{\xi}}\alpha. \end{equation}This expresses the time evolution of material particle locations in Eulerian coordinates as a function of the deformed Lagrangian coordinates. We can expand (2.18) as the system
$$\begin{gather} X_{\tau} = \frac{1}{4 \alpha^3} \left(\alpha_{\xi}Y_{\eta} - \alpha_{\eta}Y_{\xi}\right), \end{gather}$$
$$\begin{gather}Y_{\tau} = \frac{1}{4 \alpha^3} \left(\alpha_{\eta} X_{\xi} - \alpha_{\xi}X_{\eta}\right), \end{gather}$$
with 
$\alpha = X_{\xi } Y_{\eta }-X_{\eta }Y_{\xi }$. In turn, (2.19) can be rewritten as
\begin{equation} \boldsymbol{X}^{\rm T}_{\tau} ={-}\frac{1}{8}\,\boldsymbol{\nabla}_{\boldsymbol{\xi}}\boldsymbol{\cdot} \left(\frac{1}{\alpha^2} \begin{pmatrix} Y_{\eta} & -X_{\eta} \\ -Y_{\xi} & X_{\xi} \end{pmatrix} \right), \end{equation}which is an equation in divergence form that is easier to solve than (2.18) or (2.19) when using a finite-element method. Initial conditions are imposed via (2.16), so
\begin{equation} \frac{1}{\varGamma(\boldsymbol{X}(\xi,\eta,0),0)} = X_{\xi}(\xi,\eta,0)\,Y_{\eta}(\xi,\eta,0) - X_{\eta}(\xi,\eta,0)\,Y_{\xi}(\xi,\eta,0) . \end{equation}
We choose 
$Y_{\eta }(\xi,\eta, 0 )=1$ and 
$Y_{\xi }(\xi,\eta, 0 )=0$, so that 
$X_{\xi }(\xi,\eta,0)= 1/\varGamma (\boldsymbol {X}(\xi,\eta,0),0)$. This yields a purely one-dimensional transformation, as illustrated in figure 2, from the undeformed to the deformed Lagrangian domain, simplifying the calculation of the deformed geometry. The initial conditions for 
$\xi$ are therefore obtained through
\begin{equation} \int_0^x \varGamma_0(x',y) \,\mathrm{d}\kern0.7pt x' +C(\kern0.7pt y) = \xi(x,y,0), \quad \xi(0,y,0)=0. \end{equation}
Here, 
$C(\kern0.7pt y)$ is an arbitrary piecewise function chosen such that 
$\xi$ is continuous, and 
$1/X_{\xi } = \xi _x$ because we have fixed 
$Y$ and 
$t$. After finding the indefinite partial integral (2.22), we substitute 
$y=\eta$ and 
$x=X$, and invert (2.22) (numerically if needed) to find 
$X$ as an explicit function of 
$(\xi,\eta )$. Calling this solution 
$G(\xi,\eta )$, the initial conditions can be summarised as
\begin{equation} Y(\xi,\eta,0) = \eta, \quad X(\xi,\eta,0) = G(\xi,\eta). \end{equation}
Figure 2. (a) The Eulerian domain of the dynamic Lagrangian problem presented in § 2.2. This domain is broken into nine different regions (denoted R1 to R9) to compute the piecewise continuous definition of 
$\xi$, given in § S1 of the supplementary material. The red circles are the locations of the initial deposits of exogenous surfactant. (b) The deformed Lagrangian domain, calculated such that (2.16) holds for the Eulerian initial conditions (2.7) and (2.9) with the parameter choices taken in § 2.1.3. This is the domain in which we compute the numerical solution of (2.20) with boundary conditions (2.26).
2.2.2. The three-deposit problem
We illustrate the Lagrangian method introduced in § 2.2.1 by solving the model problem with the parameters outlined in § 2.1.2. For the initial conditions (2.7) and (2.9), the solution of (2.22) is
\begin{align} \xi = \begin{cases} x-
\left(\dfrac{(x-x_1)^3}{3}+x(\kern0.7pt
y-y_1)^2\right)\left(1-\delta\right) +C_1(\kern0.7pt y), &
|\boldsymbol{x}-\boldsymbol{x}_1| \leq 1,\\ \varGamma_2x-
\dfrac{1}{r_2^2}\left(\dfrac{(x-x_2)^3}{3}+x(\kern0.7pt
y-y_2)^2\right)\left(\varGamma_2-\delta\right)
+C_2(\kern0.7pt y), & |\boldsymbol{x}-\boldsymbol{x}_2|
\leq r_2,\\ \varGamma_3 x-
\dfrac{1}{r_3^2}\left(\dfrac{(x-x_3)^3}{3}+x(\kern0.7pt
y-y_3)^2\right)\left(\varGamma_3-\delta\right)
+C_3(\kern0.7pt y), & |\boldsymbol{x}-\boldsymbol{x}_3|
\leq r_3,\\ \delta x + C_4(\kern0.7pt y), &
\text{everywhere else in } \varOmega. \end{cases}
\end{align}
Here, 
$C_1(\kern0.7pt y)$, 
$C_2(\kern0.7pt y)$, 
$C_3(\kern0.7pt y)$ and 
$C_4(\kern0.7pt y)$ are determined for the choice 
$\varGamma _2=1$, 
$\varGamma _3=2$, 
$(x_1,y_1) = (6,2)$, 
$(x_2,y_2) = (10,5)$ and 
$(x_3,y_3)=(4,7)$ in § S1 of the supplementary material, along with the definition of the Lagrangian coordinates of the three circles. Finding this initial condition involves breaking the Lagrangian domain into nine regions, as shown in figure 2. By imposing 
$Y=\eta$, and imposing that the line 
$X=0$ corresponds to 
$\xi =0$, only the right-hand side of the Lagrangian domain, which we call 
$\partial \varOmega _R$ (defined for this problem in equation (S1.2) of the supplementary material), is not a straight line. We substitute 
$\eta =y$ into (2.24) and then invert (2.24) numerically to find the initial expression for 
$X$ as an explicit function of 
$\xi$ and 
$\eta$.
The boundary conditions for (2.20), and for the steady-state problem presented in the next subsection, are derived from the dynamic boundary condition (2.2). Analysis in Appendix C reveals that for corner angles less than 
${\rm \pi}$, such as we have in the domain that we consider, a particle that begins on one of the four edges of the rectangle must stay on that edge for all time, and the appropriate boundary conditions accompanying (2.27) in the undeformed Lagrangian domain are the Dirichlet conditions
\begin{equation} X = 0, L_1 \text{ on } x_0=0,L_1 \quad \text{and} \quad Y=0, L_2 \text{ on } y_0=0,L_2, \end{equation}which ensures that (2.2) is satisfied. This means that in the deformed Lagrangian domain,
\begin{equation} X(0,\eta,\tau) = 0, \quad X(\xi,\eta,\tau) = L_1 \text{ on } \partial \varOmega_R, \quad Y(\xi,0,\tau) = 0, \quad Y(\xi,L_2,\tau) = L_2. \end{equation}2.2.3. Numerical solution
Having inverted (2.24) numerically to find the initial conditions (2.23a,b), we use these initial conditions to solve (2.20) subject to boundary conditions (2.26), from 
$\tau =0$ to a final time taken to approximate the steady-state 
$\tau =t_f$, in the Lagrangian domain shown in figure 2(b), using COMSOL
$^{\circledR}$. For reproducibility purposes we provide the details of the COMSOL
$^{\circledR}$ settings chosen: we use the Mathematics suite, using the coefficient form PDE set-up that is designed to handle PDEs in divergence form such as (2.20). We discretise using standard COMSOL
$^{\circledR}$ triangulation method, and we use quadratic Lagrange basis functions with 
$314\,198$ degrees of freedom plus 
$16\,578$ internal degrees of freedom, and set the relative tolerance to 
$10^{-9}$. We store the solution at every 
$2$ time units.
2.3. The steady-state problem
2.3.1. Formulation of the problem
We now consider the problem of approximating the equilibrium locations of surface particles (their locations as 
$t\to \infty$) as a function of their initial locations directly, i.e. without any intermediate calculation of surfactant concentrations, or intermediate calculation of particle trajectories. We return to the coordinate systems used in § 2.2; however, here we revert to calling the Lagrangian coordinates 
$(x_0,y_0)$ to indicate that the Lagrangian domain is a copy of the Eulerian domain, defined as 
$0\leq x_0\leq L_1$ and 
$0\leq y_0\leq L_2$, and 
$(X,Y) = (x_0,y_0)$ at 
$t=\tau =0$. Using these variables, at steady state, (2.15) becomes
\begin{equation} X_{x_0}Y_{y_0} - X_{y_0}Y_{x_0} = \frac{\varGamma_0(x_0,y_0)}{\bar{\varGamma}}, \end{equation}
which is a PDE describing the mapping function 
$(X,Y)$ to the spatial coordinates 
$(x,y)$ for particles starting at 
$(x_0,y_0)$, in the limit 
$t\to \infty$. Equation (2.27) needs to be solved subject to boundary conditions (2.25a,b).
For one-dimensional problems, e.g. 
$\varGamma _0=\varGamma _0(x_0)$, we can impose 
$Y=y_0$ and (2.27) has a unique solution. However, in two dimensions, (2.27) constitutes only one equation for the two unknowns 
$(X,Y)$, and therefore does not have a unique solution, so we turn to the Helmholtz decomposition theorem to make progress. By this theorem, we know that we can write the map 
$\boldsymbol {X}=[X,Y]^{\rm T}$ in terms of two scalar potentials 
$\phi (x_0,y_0)$ and 
$\psi (x_0,y_0)$ such that
\begin{equation} [X,Y]^{\rm T} = \boldsymbol{\nabla}_{\boldsymbol{x}_0}\phi + \boldsymbol{\nabla}_{\boldsymbol{x}_0}\times\boldsymbol{\psi}, \end{equation}
where 
$\boldsymbol {\psi }$ is a vector of magnitude 
$\psi$ pointing out of the plane (in the 
$z$-direction), with 
$\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot } \boldsymbol {X}=\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi$ and 
$\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times \boldsymbol {X} = - \boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\boldsymbol {\psi }$. To make this Helmholtz decomposition unique up to constants, we impose the boundary conditions
\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{x}_0}\phi \boldsymbol{\cdot} \boldsymbol{n}_b = [x_0,y_0]^{\rm T}\boldsymbol{\cdot} \boldsymbol{n}_b, \quad \boldsymbol{\nabla}_{\boldsymbol{x}_0}\times \boldsymbol{\psi}\boldsymbol{\cdot} \boldsymbol{n}_b = 0 \quad \text{on } \partial \varOmega , \end{equation}which satisfies (2.25a,b).
The map at time 
$t$ is generated by (2.6), the right-hand-side of which is an Eulerian gradient of the instantaneous surfactant concentration. Thus the map remains irrotational with respect to the Eulerian coordinates. Now we investigate whether the map at time 
$t$ can be approximated by a map that is irrotational with respect to the Lagrangian coordinates, as this would allow us to remove the indeterminacy in (2.27), since 
$\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times [X,Y]^{\rm T} =\boldsymbol {0}$ yields 
$\psi$ equal to a constant, reducing the problem (2.27) to finding a solution for a single scalar potential 
$\phi$. We summarise the statement that we want to test as that, for all time 
$t$,
\begin{equation} |\boldsymbol{\nabla}_{\boldsymbol{x}_0}\phi| \gg |\boldsymbol{\nabla}_{\boldsymbol{x}_0}\times\boldsymbol{\psi}|. \end{equation}
In effect, we test the idea that because the Eulerian curl of 
$\boldsymbol {u}_s$ is zero, and material particles on boundaries are not allowed to traverse corners, (2.30) might hold for all time, at least when the rearrangement of the surface is small. We will test this hypothesis a posteriori in § 3.
Assuming that the map (2.28) is given by 
$[X,Y]^{\rm T} = \boldsymbol {\nabla }_{\boldsymbol {x}_0}\phi$, (2.27) and boundary conditions (2.25a,b) reduce to the Monge–Ampère equation
\begin{equation} \phi_{x_0 x_0}\phi_{y_0 y_0} - \phi_{x_0 y_0}^2 = \frac{\varGamma_0(x_0, y_0)}{\bar{\varGamma}} \quad \text{on } 0\leq x_0 \leq L_1, 0\leq y_0 \leq L_2 , \end{equation}subject to
\begin{equation} \phi_{x_0} = x_0 \text{ on } x_0 = 0,L_1, \quad \phi_{y_0}=y_0 \text{ on } y_0=0,L_2, \quad \phi(0,0)=0. \end{equation}
The last boundary condition is necessary to close the problem, as 
$\phi$ is unique only up to a constant. The Monge–Ampère equation arises often in the theory of optimal transport, a connection that we will discuss further in § 4.
2.3.2. Numerical method
We solve (2.31) subject to the boundary conditions (2.32) for the initial concentration profile of surfactant (2.7) and (2.9) using an iterative Newton–Raphson scheme for a finite-difference approximation of the solution, the full details of which are in Appendix D. The Newton–Raphson scheme converges to the desired solution only if the initial guess is in the basin of attraction of the desired solution, which for a nonlinear problem such as (2.31) and (2.32) is difficult to determine a priori. We surmount this problem with the following continuation scheme. Using a parameter 
$\beta _j \in [0,1]$, we take advantage of the fact that the PDE
\begin{equation} \phi_{x_0 x_0}\phi_{y_0 y_0} - \phi_{x_0 y_0}^2 - 1 + \beta_j \left(1 - \frac{\varGamma_0(x_0,y_0)}{\bar{\varGamma}} \right) = 0, \end{equation}
subject to boundary conditions (2.32), has a known solution when 
$\beta _j = 0$, namely 
$\phi = x_0^2/2 + y_0^2/2$; when 
$\beta _j=1$, we have the desired solution to (2.31) and (2.32). We step from 
$\beta _0=0$ to 
$\beta _J=1$, in steps of some fixed quantity 
$\Delta \beta =1/J$ (where 
$J$ is an integer), solving (2.33) and (2.32) each time. Starting from 
$\beta _0=0$ and 
$\phi _0=x_0^2/2+y_0^2/2$, we find 
$\phi _{j+1}$ by using 
$\phi _j$ as a guess solution for (2.33) and (2.32), where 
$\beta _{j+1}=\beta _j +\Delta \beta$. If we choose 
$\Delta \beta$ to be small enough, then we ensure that we stay inside the basin of attraction of solutions, finding the desired solution to (2.31) and (2.32) when 
$j=J$.
We use this process to solve (2.31) and (2.32) for intermediate and low values of endogenous surfactant, 
$\delta =0.25$ and 
$\delta = 0.002$. We solve for the larger value of 
$\delta$ using a grid with grid points spaced uniformly 0.05 units apart in MATLAB
$^{\circledR}$, using the software's ‘sparse’ variable type to handle the large sparse matrices, and its efficient algorithms for finding solutions to linear systems such as (D3) with a direct LU factorisation scheme. This solution is obtained by using 
$\Delta \beta =0.1$. For the solution with the smaller value of 
$\delta$, we use grid points spaced evenly 0.05 units apart. We need 
$\Delta \beta =0.0025$ for this second solution, which means that the computational cost is increased. The convergence of the numerical scheme is presented in § D.2. In addition, we present a method for creating a computational Suminagashi picture in Appendix E.
To quantify how well the Monge–Ampère method approximates the solution found by the Eulerian particle-tracking method at 
$t=t_f$ (assumed to be an accurate solution of the steady state), we define metrics that characterise the difference between solutions found using the two methods for the same initial conditions. We define the Euclidean distance between final particle locations 
$X_{EU}$ and 
$X_{MA}$ predicted by both methods and normalised by the longest side of the domain,
\begin{equation} \frac{1}{L_1}\,|{X}_{EU}-{X}_{MA}| \equiv \frac{1}{L_1}\,\sqrt{(X_{EU}-X_{MA})^2+(Y_{EU}-Y_{MA})^2}, \end{equation}which we call the normalised absolute error between the two methods for a given initial particle location. Statistics of the error are then obtained by analysing distributions for a large number of the initial particle locations, particularly the median, the upper quartile, the 90th percentile and the maximum values of (2.34).
3. Results
Table 1 summarises all of the simulations and their parameters that are presented in the results section, with a key with which we refer to each simulation.
Table 1. A table presenting a summary of the simulations presented in § 3, together with parameters used, and a key with which we refer to each simulation. The methods used are the Eulerian particle-tracking method (2.1) and (2.6), the Lagrangian particle-tracking method (2.20), and the Monge–Ampère method (2.31). For all these simulations, we choose 
$r_2=2$, 
$r_3=3$, 
$\varGamma _2=1$ and 
$\varGamma _3=2$.
3.1. Particle-tracking solutions
The results for the Eulerian (Eul[0.25]) and Lagrangian (Lag[0.25]) particle-tracking methods (presented in §§ 2.1.3 and 2.2.3) with 
$\delta =0.25$ are shown in figures 3(a) and 3(b), respectively, and also as supplementary movies 1 and 2, where each thin coloured line represents a particle trajectory, terminating at a black dot at 
$t=t_f$. The trajectories shown in figure 3(a) represent 1 in every 225 trajectories calculated, selected such that their initial locations are evenly spaced. The data obtained by the solution for the Lagrangian method presented in figure 2(b) are spaced irregularly, with each data point corresponding to a node of the mesh used in COMSOL
$^{\circledR}$ to discretise the deformed Lagrangian domain; we display 1 in every 50 particles from the data list obtained from the simulation, so the density of particles shown is not significant.
Figure 3. Solution to the example problem of three circular deposits of exogenous surfactant spreading together with 
$\delta =0.25$. (a) The results of Eul[0.25] (see supplementary movie 1). The initial boundaries of the exogenous surfactant circular deposits are the black circles, and the final locations are the thick red lines. The points described by the green curves map to the black circles at 
$t=t_f$. Individual particle trajectories are plotted using thin coloured lines terminating at black points. (b) The results of Lag[0.25] with the same colour scheme as in (a) (see supplementary movie 2). The particles represent 
$1/50$ of all the particle trajectories calculated, which are chosen at random, so the density of particles shown is not significant. (c) Graph showing the results of MA[0.25] overlaid onto Eul[0.25] and Lag[0.25]. The steady-state boundaries of the three deposits and the curves found by the inverse map (which spread from and to the black circles in the steady state, respectively) are given by the colour scheme shown in the figure legend.
In figures 3(a) and 3(b), the largest deposit spreads out through Marangoni stresses, and compresses the other two deposits. Flow reversals (sharp turns of particle trajectories of more than 
$90^{\circ }$) arise in several areas for two reasons. First, reversals in the top left-hand corner are due to confinement. Early outward spreading is into a region containing endogenous surfactant at low concentration 
$\delta$; later reversals arise once the surfactant concentration in this region is much larger due to non-local compression of the endogenous material. Second, points that begin on the edges of the smaller two deposits nearest the centre of the domain first spread into the centre, but soon the effect of the largest deposit spreading is felt, and these points reverse their trajectories. The final shapes of the smallest deposits are non-trivial oval shapes, the centres of which are shifted away from their initial locations. Some particles to the top left of the centre of the largest deposit traverse distances close to 
$1$ unit in length and then move approximately the same distance back, close to where the particles started. Particles compress into the top and bottom right-hand corners. A variety of trajectories are evident: for example, particles in the top and left have trajectories that involve straight lines and sharp turns, whereas particles towards the bottom right describe gentle arcs. In figure S1 in § S2 of the supplementary material, we present an overlay of the contour plots of 
$X(x_0,y_0)$ and 
$Y(x_0,y_0)$ for the solutions at 
$t=t_f$ obtained from the Eulerian and Lagrangian particle-tracking methods, respectively. The methods find the same particle locations to within a distance of 
$0.05$ almost everywhere, apart from the locations of small oscillations in the Lagrangian solution that appear to be an artefact of the domain deformation as discussed in § S2, and much closer than that in most places. Some of the small discrepancies that do exist can be explained partly by the fact that small errors arise by interpolating the Lagrangian solution onto a regular, rectangular grid to make the comparison, and errors occur in the Eulerian solution by the interpolation of the gradient of the evolving concentration shown in § S3 (figure S2) of the supplementary material at every time step in that solution.
3.2. The steady-state solution
We investigated the steady-state solution of a variety of configurations involving one and two deposits (see table S1 in the supplementary material). We compare the results obtained between Eulerian (EU) particle tracking and the Monge–Ampère (MA) method in § S5 of the supplementary material. We tested how the final equilibrium shape of the deposits is influenced by the proximity of the domain boundaries. In the case of a single initial deposit, we find only small differences in the discrepancy between the EU and MA results, quantified using (2.34) for deposit locations at various distances from the domain boundaries. The median normalised absolute error is approximately 
$10^{-4}$, and the maximum error is bounded by 
$2\times 10^{-3}$. In the case of two initial deposits, the median normalised absolute error is approximately 
$5\times 10^{-4}$, and the maximum error is bounded by 
$5\times 10^{-3}$. The median discrepancy between the two methods tends to be inversely correlated with the symmetry of the initial configuration, whereas the upper quartile, 90th percentile and maximum discrepancy are much noisier for both the one-deposit and two-deposit problems studied. Discrepancies between the EU and MA methods increase with an increase in the number of deposits, and with a decrease in 
$\delta$ (the normalised initial endogenous surfactant concentration), as shown in Appendix F. As stated previously, we choose to focus on the three-deposit case.
3.2.1. The three-deposit problem with 
$\delta =0.25$
The solution for the approximation of the edges of the three deposits in the steady state (MA[0.25]) found by the Monge–Ampère method (outlined in § 2.3.2) is presented in figure 3(c) for 
$\delta = 0.25$, with the approximations of the steady state found from the Eulerian (Eul[0.25]) and Lagrangian (Lag[0.25]) particle-tracking solutions overlaid. The final edges of the deposits predicted by MA[0.25] are almost indistinguishable except in a few places, which supports assumption (2.30). Predictions of Eul[0.25] and Lag[0.25] are indistinguishable to the naked eye in figure S1 of the supplementary material, providing a reliable benchmark against which to test the prediction of MA[0.25]. We also use the inverse maps to calculate the contours that map to the initial drop boundaries under the spreading in figure 3; again, only very small discrepancies between MA[0.25] and Eul[0.25] are evident.
A comparison of the global behaviour of the Monge–Ampère approximation of the map from initial to final particle configuration (MA[0.25]) with the map calculated from the particle-tracking solution (Eul[0.25]) is given by contour plots in figure 4 (see also a colour map of the absolute error between the two predictions for the final particle location in § S4, figure S3, of the supplementary material). Dense contours in figures 4(a,b) indicate that surface areas starting at these locations are stretched by the mapping, and similarly large gaps between contours indicate that the map compresses the surface. Conversely, in figures 4(c,d), dense contours of the inverse maps indicate that surface areas finishing at these locations have been compressed by the spreading, and large gaps between contours indicate that the spreading has stretched the surface. The Monge–Ampère approximation agrees with the particle-tracking solution in most places, although some noticeable discrepancies exist, such as in the left half of the smallest deposit most notably. The median error across the solution is approximately 0.25 % of the domain length, and the error for every particle is within 
$1.5$% of the domain length, as shown in figure S7 of the supplementary material.
Figure 4. Contour plots of solutions for the map from initial configuration to steady state found from MA[0.25] and Eul[0.25]. (a) 25 evenly spaced contours of 
$X_{MA}$ taken from MA[0.25] (red) overlaid with the same-valued contours of 
$X_{EU}$ taken from Eul[0.25] (blue). (b) 25 contours of 
$Y_{MA}$ (red) overlaid with 
$Y_{EU}$ (blue). (c) The inverse map 
$X_{MA}^{-1}$ (red) overlaid with 
$X_{EU}^{-1}$ (blue), with the same contour scheme. (d) The same for 
$Y_{MA}^{-1}$ (red) overlaid with 
$Y_{EU}^{-1}$ (blue).
To illustrate where and how the discrepancies arise, figure 5 shows the divergence and curl of the map found from Eul[0.25], shown in both Lagrangian and Eulerian coordinate systems. From the Helmholtz decomposition (2.28), figures 5(a,c) show 
$\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot } (\boldsymbol {X}-\boldsymbol {x}_0)=\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi -2$, and figures 5(b,d) show 
$(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp } = -\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi$, where 
$(\cdot )_{\perp }$ means the 
$z$-component, perpendicular to the plane of the solution. The vector field 
$\boldsymbol {X}-\boldsymbol {x}_0$, which points from initial to final particle locations, is an easier quantity to interpret physically than 
$\boldsymbol {X}$ itself. Given boundary conditions (2.29) (the boundary condition for 
$\psi$ can be taken to be equivalent to the Dirichlet condition 
$\psi =0$ on all four boundaries), the fact that 
$|\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi |$ is an order of magnitude greater than 
$| \boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi |$ almost everywhere is further evidence justifying our assumption (2.30) (the ratio 
$\nabla ^2_{x_0}\psi /\nabla ^2_{x_0}\phi$ is plotted in § S4, figure S4, of the supplementary material), although it is certainly not the case that the curl of the map vanishes. In figures 5(a,c), 
$\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot } (\boldsymbol {X}-\boldsymbol {x}_0)<0$ represents surface areas with net compression, and 
$\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot } (\boldsymbol {X}-\boldsymbol {x}_0)>0$ represents areas with net expansion by the map. Areas within the deposits expand, as do area elements connecting the largest deposit with the smaller deposits, and the corner regions compress. Saddle-like area elements directly between each of the initial deposit locations, and between each deposit and the nearest boundary, are compressed and expanded in orthogonal directions. In figures 5(b,d), positive values of 
$(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp }$ refer to anticlockwise net local rotation (twist) by the map, and negative values for clockwise twist. It is notable that weak twisting motions arise where interfaces spread towards a nearby boundary, or near another drop interface, with regions of oppositely oriented twisting typically appearing in pairs. The most intense twisting appears to be confined to regions immediately outside the boundaries of the exogenous surfactant drops.
Figure 5. Contour plots showing the divergence and curl of the vector field from initial to final particle location taken from Eul[0.25]. (a) Plot of 
$\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot }(\boldsymbol {X}-\boldsymbol {x}_0)= \boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi -2$ in Lagrangian coordinates. (b) Plot of 
$(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp } = -\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi$ in Lagrangian coordinates. (c) Plot of 
$\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot }(\boldsymbol {X}-\boldsymbol {x}_0)= \boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi -2$ in Eulerian coordinates. (d) Plot of 
$(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp } = -\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi$ in Eulerian coordinates.
3.2.2. Weak endogenous surfactant
Many fluids in the environment have low levels of contaminant surfactant, and in some cases (especially in controlled laboratory conditions) vanishingly small endogenous surfactant concentrations, so we would like to understand the Lagrangian motion of surface particles in the limit of very small 
$\delta$. The Lagrangian dynamic method (§ 2.2.3) is not capable of handling small values of 
$\delta$, because the deformed Lagrangian domain becomes extremely stretched. The Eulerian particle-tracking method (§ 2.1.3), and the Monge–Ampère approximation (§ 2.3.2), however, can both find well-behaved results with extremely small 
$\delta$.
We present the solution to the Monge–Ampère approximation MA[0.002] overlaid with the solution for the Eulerian particle-tracking method Eul[0.002] in figure 6. Figures 6(a,b) show the overlaid contour plots of the approximations for the steady-state solutions for 
$X$ and 
$Y$ in the Lagrangian coordinates, and figures 6(c,d) show the inverse mapping in the Eulerian coordinates. Figure 6(e) shows an overlay of the prediction for the final boundary locations for the three deposits using both methods.
Figure 6. Plots showing a comparison between MA[0.002] and Eul[0.002]. (a) Contour plots of 
$X_{MA}$ taken from MA[0.002] (red) overlaid with the same-valued contours from 
$X_{EU}$ (blue) taken from Eul[0.002]. (b) The same for 
$Y_{MA}$ (red) and 
$Y_{EU}$ (blue). (c) The same scheme for the inverse maps 
$X_{MA}^{-1}$ (red) and 
$X_{EU}^{-1}$ (blue). (d) Similarly for 
$Y_{MA}^{-1}$ (red) and 
$Y_{EU}^{-1}$ (blue). (e) An overlay of the final deposit boundaries from MA[0.002] (red) and Eul[0.002] (blue). ( f) Particle trajectories, each given by a thin coloured line terminating at a black dot, from Eul[0.002] (see supplementary movie 3). The three red dashed ellipses each contain a complete particle trajectory that involves two sharp changes of direction.
Figure 6 shows that when 
$\delta$ is small, the Monge–Ampère method does less well in approximating the solution, which is expected as particles spread further, leading to large discrepancies between the Lagrangian and Eulerian curls of the maps; however, it is still credible as a first-order approximation. In particular, for particles that are initially located within the three deposits, the approximation is still accurate, with the largest discrepancy occurring for particles of endogenous surfactant. This is shown clearly when comparing figures 6(c,d), which show the inverse map in the Eulerian coordinates, to figures 6(a,b), which show the map in the Lagrangian coordinates. (This is also shown clearly in the plot of the absolute and relative errors between the predictions for final particle locations in § S4, figures S3(b,d), of the supplementary material.) Contours coincide over an appreciably greater region of the final configuration in comparison to the initial configuration. In Appendix F, we present box-and-whisker plots of the absolute error between final particle locations predicted by the Monge–Ampère and Eulerian particle-tracking methods for several values of 
$\delta$, which show how the maximum discrepancy grows as 
$\delta \to 0$ (10 % of the domain length for 
$\delta =0.002$), but the median discrepancy remains an order of magnitude or more smaller than the maximum discrepancy for 
$\delta =0.002$.
The final locations of particles beginning on the region occupied by endogenous surfactant are compressed into effectively three lines, as shown by the blue curves in figure 6(e) for Eul[0.002]. The Monge–Ampère approximation (MA[0.002]) for the location of the lines (in red) is relatively accurate in most places. Figure 6( f) and supplementary movie 3 show the particle trajectories from the Eulerian particle-tracking method. The variety of trajectories is remarkable, with many particles having two sharp changes in direction during the spreading (see, for example, the three trajectories highlighted by red dashed ellipses). This likely reflects the fact that particle trajectories can be influenced by different deposits at different times.
Figure 7 shows the divergence and curl of the map computed from Eul[0.002]. The fact that the Monge–Ampère approximation is less accurate for small 
$\delta$ is reflected in the fact that 
$|\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi |$ and 
$|\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi |$ are shown to be of the same order of magnitude in certain places within the initial configuration (the ratio 
$\nabla ^2_{x_0}\psi /\nabla ^2_{x_0}\phi$ is plotted in figure S4 in § S4 of the supplementary material). However, figures 7(c,d), plotted in the final configuration, once again show that, for particles that start within the deposits, it is still the case that (2.30) holds. Once again, 
$\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot } (\boldsymbol {X}-\boldsymbol {x}_0)<0$ represents area elements with net compression, and 
$\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot } (\boldsymbol {X}-\boldsymbol {x}_0)>0$ represents area elements with net expansion by the map. Areas within the deposits expand, as do area elements connecting the largest deposit with the smaller deposits, and we also now see expansion for areas connecting the deposits with the boundaries. In figure 7(b), the patterns created by 
$(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp }$ are similar to the case where 
$\delta =0.25$, although the modulus of twist is greater in magnitude.
Figure 7. Contour plots showing the divergence and curl of the vector field from initial to final particle location taken from Eul[0.002]. (a) Plot of 
$\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot }(\boldsymbol {X}-\boldsymbol {x}_0)= \boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi -2$ in Lagrangian coordinates. (b) Plot of 
$(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp } = -\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi$ in Lagrangian coordinates. (c) Plot of 
$\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot }(\boldsymbol {X}-\boldsymbol {x}_0) = \boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi -2$ in Eulerian coordinates. (d) Plot of 
$(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp } = -\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi$ in Eulerian coordinates.
3.2.3. Self-similarity of corner regions
Computationally, the Monge–Ampère solution is significantly cheaper to solve compared to either of the particle-tracking solutions. We take advantage of this to analyse the shape of the corners of deposits as 
$\delta \to 0$ using finely discretised Monge–Ampère approximations. We note that the errors found between the Eulerian and Monge–Ampère predictions are small near the corner regions, even with 
$\delta =0.002$ (see figure S3 in the supplementary material). This suggests that predictions from the Monge–Ampère solution remain accurate in the corner regions in the limit 
$\delta \to 0$. Examples of these corner regions are numbered in green in figure 8(a), which shows the results of MA[0.04]. To compute the deposit edges as sufficiently smooth curves to accurately analyse curvature of the edges was not possible from the particle-tracking solutions with available computational power, but is possible with the Monge–Ampère method. During spreading, endogenous surfactant at initial concentration 
$\delta$, occupying an 
$O(1)$ area, is compressed into narrow threads (between pairs of drops) and into what resemble seven Plateau borders (at corners between drops, and at the domain boundary), at final concentration of order unity. Assuming that the bulk of the endogenous surfactant is driven into these seven regions, and noting that each has an 
$O(1)$ aspect ratio, we can expect each drop corner to curve over a length scale of 
$O(\delta ^{1/2})$. This scaling is motivated by the observation that an area with an 
$O(1)$ length scale and with an initial concentration of endogenous surfactant 
$\delta$ (giving a total mass of 
$O(\delta )$) is ultimately compressed into an area with final concentration 
$\bar {\varGamma }\sim 1$ over an 
$O(\delta ^{1/2})$ length scale (preserving the total mass of 
$O(\delta )$).
Figure 8. Scaled curvature plots for a selection of the corners of the boundaries of the three circular deposits at the steady state for small values of 
$\delta$. (a) A graph of the solution from MA[0.04] that shows in green our numbering system for corners inside each deposit. (b–f) Plots of the natural logarithm of the curvature scaled by 
$\delta ^{1/2}$ against the arc length scaled by 
$\delta ^{-1/2}$, where 
$s=0$ identifies the vertex of the corner in each case.
Figure 8 shows how the equilibrium shapes of the deposits calculated by Monge–Ampère become more and more polygonal in the limit 
$\delta \to 0$, with corner regions adopting a self-similar form. For several values of small 
$\delta$, we calculate the curvature 
$K$ of the equilibrium boundary of each deposit as a function of arc length 
$s$. To show how each of a selection of corners sharpens as 
$\delta \to 0$, we set 
$s=0$ at the local curvature maximum, and plot 
$\log {(\delta ^{1/2}K)}$ against 
$s\delta ^{1/2}$. Figures 8(b–f) show collapse of the data, illustrating how, as 
$\delta \to 0$, the curvature of corners of the polygons (numbered in green in figure 8a) is proportional to 
$\delta ^{-1/2}$, with this curvature occurring over arc lengths proportional to 
$\delta ^{1/2}$, and the edges of the deposits becoming effectively straight lines away from the corners. We see that in small regions away from 
$s=0$, the quantity 
$\ln {(\delta ^{1/2}K)}$ becomes linear, indicating a functional dependence such as 
$\delta ^{1/2}K \sim \exp (-\lambda \delta ^{-1/2}|s|)$ for some constant 
$\lambda$. In practice, the boundaries between the deposits have a small curvature in the Eulerian particle-tracking solution (figure 6e), so the deposits instead resemble slightly distorted polygons, with corner regions slightly distorted from the Monge–Ampère prediction. As discussed above, the self-similarity observed for the internal corners, and the appearance of a characteristic length scale of order 
$\delta ^{1/2}$, are based on the fundamental physical principle of mass conservation. Since this principle is independent of the choice of method, we expect that the results observed using the Monge–Ampère method will remain valid with the Eulerian and Lagrangian methods. This provides further evidence that the self-similarity behaviour observed in figure 8 is physically grounded.
Figure 9 shows the Monge–Ampère approximation of the three-deposit problem with 
$\delta = 0.005$ (MA[0.005]Alt1–Alt6), using multiple different centres of the deposits (given in the caption), revealing a variety of near-polygonal structures. Each deposit approaches a triangle, a quadrilateral, or in the case of the largest deposit (figures 9b,d), a hexagon. The final configuration approaches a structure that can be characterised by up to four coordinates, i.e. the locations of corners shared by multiple deposits, which we call the characteristic coordinates of a given initial configuration. As in figure 8, we expect the Monge–Ampère approximation to provide a good first-order estimate of equilibrium drop shape predicted by the full dynamic problem. However, at the present time we are unable to offer any simple strategy for determining the characteristic coordinates directly from the initial conditions.
Figure 9. Plots showing the steady-state mapping and inverse mapping for 
$\delta =0.005$ with varying initial locations for some of the circular deposits 1 and 2 (deposit 3, as shown in figure 8a, remains fixed). Initial locations of the boundaries of the circular deposits are in blue, with the final locations of those boundaries in red. The green curves map to the blue circles under the same map. (a–f) The results for MA[0.005]Alt1 to MA[0.005]Alt6, respectively, from the key shown in table 1.
4. Discussion
In this study, we have evaluated the trajectories of passive surface particles confined in a rectangular domain under the action of surfactant spreading on a thin film in a large-gravity and large-Péclet-number limit using two separate methods (Eulerian particle tracking, outlined in § 2.1, and Lagrangian particle tracking, outlined in § 2.2), the results of which corroborate each other. We have also identified a direct method (using the Monge–Ampère equation outlined in § 2.3) for approximating the equilibrium configuration as 
$t\to \infty$. The solutions to this problem show how confinement and drop–drop interactions can lead to drift of surface particles (in addition to spreading), and how drop–boundary interactions lead to transient flow reversals. It is striking the way that some particles move a significant distance from their initial location, before moving back close to where they started, and that in figure 6( f) some particles sharply change directions twice, which is reminiscent of the multiple regimes identified for multiple surfactant sources noted by Iasella et al. (Reference Iasella, Sharma, Garoff and Tilton2024). Predicting the equilibrium location of the interface between exogenous and endogenous surfactants is straightforward for spreading in one spatial dimension (relying on mass conservation arguments). However, in two-dimensions (as here) it becomes a non-trivial task, even for simple scenarios with one or two initial deposits, particularly when symmetry is broken. With three initial deposits, new topological features appear in the form of internal corners between deposits. Nevertheless, equilibrium shapes possess asymptotic structures, such as the self-similar geometry of the internal corners in the limit of small initial endogenous surfactant concentration (figures 8 and 9). Such observations might offer a route to explicit predictions.
In order to compute the particle trajectories involved in the surfactant-driven spreading, we constructed a finite-difference numerical method on a regular grid. We used an interpolated concentration gradient to determine the local velocity of the particles. The surfactant concentration gradients were calculated using a finite-difference solution from a second grid. In addition, we reformulated the Eulerian surfactant transport equation (2.1) and the Lagrangian particle transport equation (2.6) into a single Lagrangian vector equation (2.20), which describes the dynamics of material particles. This was achieved by choosing a Lagrangian coordinate system for which the surfactant concentration is initially uniform in the deformed Lagrangian domain (i.e. having uniform mass per unit Lagrangian area), and choosing this domain such that the initial conditions are simple to compute from (2.21). Although the calculation of surfactant concentration is bypassed by this approach, the concentration can still be recovered at any time from the Jacobian via (2.16).
From a computational point of view, the Eulerian and Lagrangian numerical methods used to compute the particle trajectories each have some advantages and drawbacks. They are similar in computational expense; however, the Eulerian method introduces additional errors at every time step by first approximating 
$\boldsymbol {\nabla }_{\boldsymbol {x}}\varGamma$ on the Eulerian grid, and then interpolating this quantity onto the current position of surface particles on the Lagrangian grid. The Eulerian method also requires storage of two solutions (Eulerian surfactant concentration and Lagrangian particle position) on two separate grids, whereas the Lagrangian method needs storage of only a single solution. Some of the drawbacks of the Lagrangian method are that: the deformed domain needs to be calculated beforehand (although this needs to be done only once); the Lagrangian domain has an irregular shape, which makes numerical resolution of (2.20) difficult; and the domain can become so deformed for small 
$\delta$ that numerical resolution of (2.20) is not practical. For moderate 
$\delta$, however, the Lagrangian formulation is a mathematically elegant approach that finds the solution without needing to introduce interpolation errors, and the close agreement of the results between the two methods (figures 4 and S1) provides corroboration of the particle-tracking method, lending support to our calculations of Eulerian solutions for small 
$\delta$.
Since equilibrium drop shapes, obtained in the limit of large time, can be computationally expensive to calculate with the Eulerian and Lagrangian methods, we have also explored how a Monge–Ampère equation can be used to find approximations of the equilibrium configurations. This was achieved by approximating the map between initial and final configuration as the gradient of a scalar potential, and neglecting the rotational part of the map. The computational benefits of this approach are beyond dispute, as only a single, scalar, time-independent PDE needs to be solved. For example, the Monge–Ampère solution shown in figure 4 can be computed for approximately 
$2.3\times 10^7$ grid points in a shorter time than a solution for approximately 
$3\times 10^5$ grid points using the Eulerian particle-tracking method from § 2.2.3 for the same parameters. Beyond computational issues, the Monge–Ampère method reveals a counterintuitive feature of surfactant-spreading dynamics in two dimensions: despite always being transported by an instantaneously irrotational flow (2.3), in the Eulerian sense, the mapping of particles from their initial to their current state can accumulate a weak rotational component with respect to Lagrangian coordinates, as illustrated in figures 5(b,d) and 7(b,d). This small Lagrangian rotational component can lead to weak distortions of final equilibrium configurations between the Monge–Ampère predictions and the Eulerian and Lagrangian predictions (figure 4). In the supplementary material, we explored a variety of alternate exogenous configurations to analyse further the accuracy of the Monge–Ampère predictions. These show how symmetry of the initial conditions reduces the median error between the Monge–Ampère approximation and the Eulerian solution. In general, the error between the Monge–Ampère approximation and the Eulerian solution reduces with larger initial endogenous surfactant concentration 
$\delta$, and with having fewer initial deposits. It is worth noting that when spreading is one-dimensional, the Monge–Ampère equation gives the exact solution. This is due to the fact that topologically, no rearrangement of the particles can be performed in one dimension without an energetic cost captured by the evolution equation. In contrast, in two dimensions, particles can rearrange through a divergence-free stress field (such as blowing by a Suminagashi artist after reaching equilibrium), which has no energetic cost for the evolution equation (2.1).
Furthermore, the Monge–Ampère formulation makes a connection between the theory of surfactant-driven transport and the theory of optimal transport. Monge–Ampère equations arise for optimal transport problems where the transport satisfies the so-called quadratic Monge–Kantorovich optimal transport problem (qMK). If surfactant were transported optimally according to qMK, then the map 
$\boldsymbol {X}$ would satisfy
\begin{equation} \min_{\boldsymbol{X}:\varGamma_0\to \bar{\varGamma}} \int_{\varOmega}|\boldsymbol{X}-\boldsymbol{x}_0|^2\varGamma_0 \, \mathrm{d} A_{\boldsymbol{x}_0}, \end{equation}
where 
$\min _{\boldsymbol {X}:\varGamma _0\to \bar {\varGamma }}$ means that we choose the minimum over all possible maps 
$\boldsymbol {X}$ that transport the initial surfactant concentration profile 
$\varGamma _0$ to the final uniform profile 
$\bar {\varGamma }$ (equivalently, from all possible maps 
$\boldsymbol {X}$ that satisfy (2.27)). It is known that such maps 
$\boldsymbol {X}$, which satisfy (4.1), are the gradients of (convex) scalar potentials that satisfy the Monge–Ampère equation (Rockafellar Reference Rockafellar1970; Caffarelli & McCann Reference Caffarelli and McCann2010). Unfortunately, the surfactant problem (2.1) does not satisfy (4.1) precisely, as we can see by the discrepancies in figure 6, for example, which shows a different prediction for the final configuration compared to the Eulerian particle-tracking method. However, Otto (Reference Otto2001) showed that solutions to equations of porous medium type (a class of equations to which (2.1) belongs) are gradient flows on the function space 
$\mathcal {M}$ of possible solutions; this function space can be shown to satisfy the definition of a Riemannian manifold when distances on the manifold are measured by the Wasserstein-2 distance, defined variationally as
\begin{equation} W_2(\varGamma_0,\bar{\varGamma}) = \sqrt{ \min_{\boldsymbol{X}:\varGamma_0\to \bar{\varGamma}} \int_{\varOmega}|\boldsymbol{X}-\boldsymbol{x}_0|^2\varGamma_0 \, \mathrm{d} A_{\boldsymbol{x}_0} }. \end{equation}
The Wasserstein-2 distance (4.2) is simply the square root of qMK in (4.1). Hence optimal maps that satisfy (4.1) must follow the shortest path between 
$\varGamma _0$ and 
$\bar {\varGamma }$ on 
$\mathcal {M}$ (a geodesic). The question of how closely the Monge–Ampère method approximates the correct steady-state solution for surfactant spreading is therefore equivalent to the question of how far the gradient flow deviates from the geodesic on 
$\mathcal {M}$. This leads to the possibility that a quantification of how well the Monge–Ampère method approximates the correct solution could be made with variational analysis. We do not pursue this further here, except to recall the degeneracy in the evolution equation revealed in Appendix A, namely that (2.1) does not account for the dissipation associated with any flow that preserves surface concentration; this raises the possibility that such flows may account for differences between equilibria predicted by the Monge–Ampère and Eulerian descriptions. Similarly, we leave as open the question of whether the solutions as 
$\delta \to 0$ can be used to infer the behaviour of the final state with an initially clean interface (
$\delta =0$).
In the limit of vanishing endogenous surfactant concentration 
$\delta \to 0$, using the Monge–Ampère approximation, we find that exogenous deposits approach equilibrium shapes with polygonal boundaries (figure 9). Moreover, the internal corners between deposits present self-similar geometries, with a typical spatial extent scaling as 
$\delta ^{1/2}$ (figure 8). This self-similar behaviour can be explained by mass conservation. In the limit of 
$\delta \to 0$, we also observe that the topology at equilibrium is reduced to a small finite number of characteristic coordinates that are the intersections between the polygonal final shapes of the deposits. Predicting these coordinates a priori remains an open problem.
In figure 10 and supplementary movie 4, we explore how the Monge–Ampère method could be used to create a pattern resembling Suminagashi art. We have used divergence-free maps to mimic the blowing process used by Suminagashi artists. (The details for the creation of this picture are in Appendix E.) Such a computation, which involves 20 consecutive calculations of the Monge–Ampère equation, would be extremely expensive to run using a full dynamic solution. The choice of a divergence-free map found by time-stepping (E7) for the blowing step used by artists is credible, as gentle blowing on a surfactant-laden surface would likely deform the surface in such a way as to not create concentration gradients. Indeed, at equilibrium, the surfactant prevents the formation of concentration gradients (Manikantan & Squires Reference Manikantan and Squires2020), because the Marangoni force opposes them (Appendix A). Moreover, we exploited the facts that gravitational forces in the Suminagashi process are sufficiently strong to suppress deformations of the gas–liquid interface, and that surface diffusion is sufficiently weak compared to advection to avoid smearing of the edge deposits. In other contexts, Marangoni flows can induce surface deformations. Despite the fact that the coupled evolution equations between the bulk and the surface for such problems can be written in a gradient flow formulation (Thiele et al. Reference Thiele, Snoeijer, Trinschek and John2018), the link between such flows and optimal transport is not clear; the bulk flow that transports fluid is distinct from the surface flow that transports surfactant, such that distinct maps may be needed to characterise such problems. The numerical results from the Monge–Ampère method presented in figure 10 are only suggestive in their reproduction of Suminagashi art; nevertheless, we hope that they will encourage experimentalists to explore this link further.
Figure 10. Successive solutions of the Monge–Ampère equation are shown in (a–e), as detailed in Appendix E, showing a final pattern that is reminiscent of a Suminagashi pattern in ( f) (see also supplementary movie 4). The image was created by 20 solutions of the Monge–Ampère equation with a stirring or blowing step after the release of every four deposits. (a–c) Creation of the pattern after four depositions, after the first stirring step, and after 
$12$ depositions, respectively. (d–f) The solution just before the final stirring step after 20 depositions, the final result where a monochrome colour scheme is added, and a Suminagashi pattern made by Bea Mahan (Reference Mahan2011) for qualitative comparison.
Many applications involving surfactants are concerned with the transport of passive solutes by the Marangoni effect in confined geometries, ranging from pharmaceutical delivery to the human lung to the creation of artistic patterns using Suminagashi techniques. We have presented methods for finding the dynamic behaviour of material surface particles, and for finding an approximation of the equilibrium location of material particles without having to solve for transient dynamics. The equilibrium approximation was achieved by showing a connection between surfactant dynamics and the theory of optimal transport, a research area at the forefront of modern mathematics, through a Monge–Ampère equation associated with the surfactant-driven transport. We hope that this connection, and the methods presented here, will spark the imaginations of researchers interested in both fundamental understanding and practical work related to surfactant-induced Marangoni flows carrying solutes.
Supplementary material
Supplementary material and movies are available at https://doi.org/10.1017/jfm.2024.334.
Acknowledgements
The authors thank all anonymous reviewers for constructive comments, and R.M. wishes to thank T. Bickel and M. Heil for constructive discussions. For the purpose of Open Access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.
Funding
J.R.L. and O.E.J. acknowledge financial support from EPSRC grant EP/T030739/1. R.M. acknowledges funding from EPSRC's PhD studentship and Doctoral Prize Fellowship schemes.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of the nonlinear diffusion equation (2.1)
We consider a liquid layer of density 
$\rho ^*$, viscosity 
$\mu ^*$ and uniform thickness 
$h^*$, sitting on the horizontal plane 
$z^*=0$, subject to a restoring force (provided, for example, by a strong vertical gravitational field 
$g^*$) that suppresses out-of-plane deflections of the gas/liquid interface at 
$z^*=h^*$ (where stars denote dimensional quantities). We therefore disregard the normal stress condition in the equations below, but provide a condition on the relevant Bond number that permits this approximation. An insoluble surfactant with concentration 
$\varGamma ^*(\boldsymbol {x}^*_\parallel,t^*)$ occupies the interface, where 
$\boldsymbol {x}^*_\parallel$ denotes horizontal or in-plane coordinates such that 
$\boldsymbol {x}^*_\parallel \equiv (x^*,y^*)$. The surfactant lowers surface tension via a linear equation of state, so that in-plane surface tension gradients are 
$-\mathcal {A}^*\,\boldsymbol {\nabla }^*_{\parallel } \varGamma ^*$, with 
$\mathcal {A}^*>0$ being the surface activity of the surfactant (Manikantan & Squires Reference Manikantan and Squires2020), 
$\boldsymbol {\nabla }^*_\parallel \equiv \hat {\boldsymbol {x}}\,\partial ^*_x +\hat {\boldsymbol {y}}\, \partial ^*_y$, and where 
$(\hat {\boldsymbol {x}},\hat {\boldsymbol {y}},\hat {\boldsymbol {z}})$ are unit vectors in the three Cartesian coordinate directions. Adopting lubrication theory, we assumed that the horizontal velocity field 
$\boldsymbol {u}^*_\parallel (\boldsymbol {x}^*_\parallel, z^*,t^*)$ in the liquid layer satisfies the Stokes equation 
$\mu ^* \boldsymbol {u}^*_{\parallel,z^*z^*}=\boldsymbol {\nabla }^*_\parallel p^*$ in 
$0\leq z^*\leq h^*$, where 
$p^*(\boldsymbol {x}^*_\parallel,t^*)$ is the leading-order pressure field. The horizontal volume flux 
$\boldsymbol {q}^*(\boldsymbol {x}^*_\parallel,t^*)=\int _0^{h^*}\boldsymbol {u}^*_\parallel \,\mathrm {d}z^*$ satisfies 
$\boldsymbol {\nabla }^*_\parallel \boldsymbol {\cdot } \boldsymbol {q}^*=0$, in order to maintain uniformity of the layer thickness. Integrating the momentum equation, applying a no-slip boundary condition 
$\boldsymbol {u}^*_\parallel =\boldsymbol {0}$ on 
$z^*=0$, and the tangential stress condition
\begin{equation} \mu^* \boldsymbol{u}^*_{{\parallel},z^*}={-}\mathcal{A}^*\,\boldsymbol{\nabla}^*_\parallel \varGamma^* +\boldsymbol{\tau}^*\quad\mathrm{at}\ z^*=h^*, \end{equation}
where 
$\boldsymbol {\tau }^*(\boldsymbol {x}^*_\parallel,t^*)$ is an imposed shear stress from the gas phase, we obtain 
$\mu ^*\boldsymbol {u}^*_\parallel =-\boldsymbol {\nabla }^*_\parallel p^* z^*(h^*-z^*/2)+(\boldsymbol {\tau }^*-\mathcal {A}^*\,\boldsymbol {\nabla }^*_\parallel \varGamma ^*)z^*$, yielding the following expressions for the flux and surface velocity:
$$\begin{gather} \mu^*\boldsymbol{q}^* ={-}\tfrac{1}{3} h^{^*3}\,\boldsymbol{\nabla}^*_\parallel p^* +\tfrac{1}{2}h^{^*2} (\boldsymbol{\tau}^*-\mathcal{A}^*\,\boldsymbol{\nabla}^*_\parallel \varGamma^*), \end{gather}$$
$$\begin{gather}\mu^* \boldsymbol{u}^*_s={-}\tfrac{1}{2}h^{^*2}\,\boldsymbol{\nabla}^*_\parallel p^* +h^*(\boldsymbol{\tau}^*-\mathcal{A}^*\,\boldsymbol{\nabla}^*_\parallel\varGamma^*). \end{gather}$$
Material particles at the interface are transported via 
$\mathrm {d}\boldsymbol {x}^*_\parallel /\mathrm {d}t^*=\boldsymbol {u}^*_s(\boldsymbol {x}^*_\parallel, t^*)$.
Using Helmholtz decomposition, we may write 
$\boldsymbol {\tau }^*=\boldsymbol {\nabla }^*_\parallel \varphi ^*+\boldsymbol {\nabla }^*_\parallel \times (\kappa ^*\hat {\boldsymbol {z}})$ for some scalar potentials 
$\varphi ^*(\boldsymbol {x}^*_\parallel, t^*)$ and 
$\kappa ^*(\boldsymbol {x}^*_\parallel,t^*)$, where 
$\boldsymbol {\nabla }^*_\parallel \boldsymbol {\cdot } \boldsymbol {\tau }^*=\boldsymbol {\nabla }^{*2}_\parallel \varphi ^*$ and 
$\boldsymbol {\nabla }^*_\parallel \times \boldsymbol {\tau }^*=\boldsymbol {\nabla }^{*2}_\parallel \kappa ^*$. Then the mass conservation constraint 
$\boldsymbol {\nabla }^*_\parallel \boldsymbol {\cdot }\boldsymbol {q}^*=0$ implies
\begin{equation} 0=\boldsymbol{\nabla}^{*2}_\parallel\big[ -\tfrac{1}{3}h^{*3} p^{*} +\tfrac{1}{2}h^{*2} ( \varphi^{*}- \mathcal{A}^{*} \varGamma^{*})\big]. \end{equation}
We impose no-flux (Neumann) conditions on the pressure, surfactant concentration and stress potential 
$\varphi ^*$ at the periphery of the domain, so that on the boundary, 
$\boldsymbol {n}_b\boldsymbol {\cdot } \boldsymbol {\nabla }^{*}_\parallel p^{*}=\boldsymbol {n}_b\boldsymbol {\cdot } \boldsymbol {\nabla }^{*}_\parallel \varphi ^{*} =\boldsymbol {n}_b\boldsymbol {\cdot } \boldsymbol {\nabla }^{*}_\parallel \varGamma ^{*}=0$, where 
$\boldsymbol {n}_b$ is the unit outward normal. Integrating (A3), the pressure gradient satisfies
\begin{equation} h^{*}\,\boldsymbol{\nabla}^{*}_\parallel p^{*}=\tfrac{3}{2}\,\boldsymbol{\nabla}^{*}_\parallel\left( \varphi^{*} - \mathcal{A}^{*} \varGamma^{*}\right ). \end{equation}Thus the surface velocity field (A2b) becomes
\begin{equation} \mu^{*} \boldsymbol{u}^{*}_s= \tfrac{1}{4}h^{*}(\boldsymbol{\nabla}^{*}_\parallel \varphi^{*} - \mathcal{A}^{*}\,\boldsymbol{\nabla}^{*}_\parallel \varGamma^{*}) + h^{*}\, \boldsymbol{\nabla}^{*}_\parallel \times (\kappa^{*} \hat{\boldsymbol{z}}), \end{equation}which can be inserted into the surfactant transport equation
\begin{equation} \varGamma^{*}_t+\boldsymbol{\nabla}^{*}_\parallel \boldsymbol{\cdot} (\boldsymbol{u}^{*}_s \varGamma^{*})=0 \end{equation}(neglecting surface diffusion). We highlight two special cases of the resulting evolution equation for the surfactant concentration.
First, in the absence of an imposed shear stress (
$\boldsymbol {\tau }^*=0$), (A5) and (A6) yield a generalisation to two dimensions of the nonlinear diffusion equation derived in Jensen & Halpern (Reference Jensen and Halpern1998):
\begin{equation} \varGamma^{*}_t=\frac{\mathcal{A}^{*}h^{*}}{4\mu^{*}}\,\boldsymbol{\nabla}^{*}_\parallel \boldsymbol{\cdot} (\varGamma^{*}\,\boldsymbol{\nabla}^{*}_\parallel \varGamma^{*}). \end{equation}As assumed previously, the condition for the interface to remain flat is that horizontal pressure gradients generated by interfacial deflections (through gravity forces, for instance) dominate those arising from Marangoni effects in (A4), namely that the Bond number is large:
\begin{equation} 1 \ll \frac{\rho^{*} g^{*} h^{*2}}{\mathcal{A}^{*}\,\Delta \varGamma^{*}}, \end{equation}
where 
$\Delta \varGamma ^{*}$ represents characteristic surfactant concentration differences driving the flow. To neglect the effects of surface diffusivity 
$D^*$ from (A7), the surface Péclet number should be large, or
\begin{equation} \frac{D^* \mu^*}{\mathcal{A}^*\,\Delta \varGamma^*\,h^*}\ll 1. \end{equation} Second, retaining the imposed shear stress 
$\boldsymbol {\tau }^*$ but decomposing the surfactant concentration such that 
$\varGamma ^{*}=\bar {\varGamma }^{*}+\hat {\varGamma }^{*}(\boldsymbol {x}^{*}_\parallel,t^{*})$, and linearising about a uniform state 
$\bar {\varGamma }^{*}>0$ (assuming 
$\vert \hat {\varGamma }^{*}\vert \ll \bar {\varGamma }^{*}$), (A5) and (A6) become
\begin{equation} \hat{\varGamma}^{*}_t- \frac{\mathcal{A}^{*}h^{*}\bar{\varGamma}^{*}}{4\mu^{*}}\,\boldsymbol{\nabla}^{*2}_\parallel \hat{\varGamma}^{*} ={-} \frac{h^{*}\bar{\varGamma}^{*}}{4\mu^{*}}\,\boldsymbol{\nabla}^{*}_\parallel \boldsymbol{\cdot} \boldsymbol{\tau}^{*}. \end{equation}
The rotational component of the shear stress (associated with 
$\kappa ^{*}$) moves material particles (via (A2b)) but does not change surface concentration in (A10) (in this linear approximation). Thus, based on (A5), there is unrestricted transport of material elements under
\begin{equation} \frac{\mathrm{d}\boldsymbol{x}^{*}_\parallel}{\mathrm{d}t^{*}}=\frac{h^{*}}{\mu^{*}} \big[\boldsymbol{\nabla}^{*}_\parallel \times (\kappa^{*} \hat{\boldsymbol{z}})\big]\big\vert_{\boldsymbol{x}^{*}_\parallel}, \end{equation}a feature that we will exploit to create Suminagashi patterns in Appendix E. In contrast, shear stress with non-zero divergence acts as a forcing to the transport equation governing surfactant concentration in (A10), which responds diffusively in the linear approximation.
The governing non-dimensional surfactant transport equation (2.1) of the problems studied in this paper (see § 2) is based on the dimensional equation (A6). Considering the initial condition illustrated in figure 1, we non-dimensionalise (A6) using a characteristic horizontal length scale 
$r_1^*$, representing the initial radius of a deposit of exogenous surfactant. We impose that the liquid height 
$h^*$ is small compared to this length scale, with the ratio given as the small parameter 
$\epsilon = h^*/r_1^* \ll 1$. We impose that the ratio of horizontal lengths 
$L_1^*/L_2^*$ is 
$O(1)$ with respect to 
$\epsilon$ (this is to ensure that the ratio of length scales does not affect the asymptotics). Surface tension 
$\gamma ^*=\gamma _0^*-\mathcal {A}^* \varGamma ^*$ is non-dimensionalised by 
$\gamma = (\gamma ^*-\gamma _c^*)/S^*$, where 
$S^* = \gamma _0^*-\gamma _c^*= \mathcal {A}^* \varGamma _c^*$, with 
$\gamma _c^*$ the surface tension when 
$\varGamma ^*=\varGamma ^*_c$, a characteristic concentration used to non-dimensionalise all surfactant concentrations, and 
$\gamma _0^*$ the nominal surface tension when 
$\varGamma ^*=0$. The surface velocity scale is obtained from the viscous–Marangoni stress balance in the dynamic boundary condition at the liquid surface (A1), which gives characteristic surface velocity 
$\epsilon S^*/\mu ^*$. The scale for the vertical velocity component 
$w^*$ is found from non-dimensionalising the continuity equation 
$\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}$ (where 
$\boldsymbol {\nabla }$ and 
$\boldsymbol {u}$ are the full three-dimensional nabla operator and velocity field) by the horizontal velocity and length scales, yielding 
$\epsilon ^2 S^*/\mu ^*$. The time scale is found by combining a horizontal velocity scale with the length scale 
$r_1^*$, yielding 
$\mu ^*r_1^*/(\epsilon S^*)$. The pressure scale that we use is 
$S^*/(\epsilon r_1^*)$. In summary, we relate unstarred dimensionless variables to starred dimensional variables by
\begin{equation} \left. \begin{gathered}
x = \frac{x^*}{r_1^*}, \quad y = \frac{y^*}{r_1^*}, \quad z
= \frac{z^*}{\epsilon r_1^*}, \quad L_i =
\frac{L_i^*}{r_1^*} \text{ for } i=1,2,\\ u = \frac{u^*
\mu^* }{\epsilon S^*}, \quad v = \frac{v^* \mu^* }{\epsilon
S^*}, \quad w = \frac{w^* \mu^* }{\epsilon^2 S^*}, \quad
\gamma = \frac{\gamma^*-\gamma_c^*}{S^*},\\ t =
\frac{\epsilon t^* S^*}{\mu^*r_1^*}\quad \varGamma =
\frac{\varGamma^*}{\varGamma_c^*}, \quad \delta =
\frac{\delta^*}{\varGamma_c^*}, \quad \bar{\varGamma} =
\frac{\bar{\varGamma}^*}{\varGamma_c^*}, \quad p =
\frac{p^* \epsilon r_1^*}{S^*}. \end{gathered} \right\} \end{equation}
Using the scales above, we can non-dimensionalise the governing surfactant transport equation (A7) to obtain (2.1).
Appendix B. Late-time correction for the mapping ${X}$
We run the Eulerian and Lagrangian methods for simulating particle trajectories from 
$t=0$ to some final time 
$t_f$, where 
$t_f$ needs to be chosen such that it gives an accurate approximation of the steady state within a certain tolerance. To find such 
$t_f$, we analyse the error between the mapping calculated numerically at 
$t_f$, and an estimate of the mapping computed theoretically using a linearised version of the governing equation as 
$t \to \infty$. At late time, the nonlinear diffusion equation (2.1) can be approximated by the linear diffusion equation
\begin{equation} \varGamma_t = \frac{\bar{\varGamma}}{4}\,\boldsymbol{\nabla}_{\boldsymbol{x}}^2\varGamma, \end{equation}
with the no-flux boundary condition (2.2) and where 
$\bar {\varGamma }$ is defined in (2.4). Using separation of variables, the solution of (B1) can be found analytically as the double series
\begin{equation} \varGamma(x,y,t) = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \sigma_{mn} \cos{\left(\frac{m {\rm \pi}x}{L_1}\right)} \cos\left(\frac{n {\rm \pi}y}{L_2}\right) {\rm e}^{-\omega_{mn}(t-t_f)}, \end{equation}
for a series of constants 
$\sigma _{mn}$, where
\begin{equation} \omega_{mn} = \frac{\bar{\varGamma}}{4} \Bigg(\frac{m^2{\rm \pi}^2}{L_1^2} + \frac{n^2{\rm \pi}^2}{L_2^2}\Bigg). \end{equation}Thus the surface velocity field is given by
\begin{equation} \boldsymbol{u}_s ={-}
\frac{1}{4}\,\boldsymbol{\nabla}_{\boldsymbol{x}}\varGamma
= \begin{pmatrix} \displaystyle\sum_{m=0}^{\infty}
\displaystyle\sum_{n=0}^{\infty} \displaystyle\frac{m {\rm \pi}
\sigma_{mn}}{4L_1} \sin{(m{\rm \pi} x/L_1)}\cos{(n {\rm \pi}y/L_2)}\,{\rm
e}^{-\omega_{mn}(t-t_f)}\\
\displaystyle\sum_{m=0}^{\infty}\displaystyle\sum_{n=0}^{\infty}\displaystyle\frac{n
{\rm \pi} \sigma_{mn}}{4L_2} \cos{(m{\rm \pi} x/L_1)}\sin{(n {\rm \pi}y/L_2)}\,{\rm
e}^{-\omega_{mn}(t-t_f)} \end{pmatrix}.
\end{equation}
We integrate this expression from 
$t=t_f$ to 
$t\to \infty$, giving a linear correction for the map of particle trajectories:
\begin{equation} \boldsymbol{X}_{cr} = \begin{pmatrix} \displaystyle\sum_{m=1}^{\infty}\displaystyle\sum_{n=0}^{\infty}\displaystyle\frac{m {\rm \pi}\sigma_{mn}}{4L_1 \omega_{mn}} \sin{(m{\rm \pi} x/L_1)}\cos{(n {\rm \pi}y/L_2)}\\ \displaystyle\sum_{m=0}^{\infty} \displaystyle\sum_{n=1}^{\infty} \displaystyle\frac{n {\rm \pi}\sigma_{mn}}{4L_2 \omega_{mn}} \cos{(m{\rm \pi} x/L_1)}\sin{(n {\rm \pi}y/L_2)} \end{pmatrix}. \end{equation}
Hence the equilibrium mapping between 
$t=0$ and 
$t\to \infty$ can be approximated by 
$\boldsymbol {X} \approx \boldsymbol {X}_{t_f}+\boldsymbol {X}_{cr}$, where 
$\boldsymbol {X}_{t_f}$ is the mapping solution from 
$t=0$ to a large time 
$t_f$, and the coefficients 
$\sigma _{mn}$ are given by
\begin{align} \sigma_{mn} =
\begin{cases}
\displaystyle\frac{1}{L_1L_2}\int_{x=0}^{L_1}\int_{y=0}^{L_2}
\varGamma(x,y,t_f) \,\mathrm{d}y\, \mathrm{d}\kern0.7pt x =
\bar{\varGamma}, & \text{for } m,n = 0,\\
\displaystyle\frac{2}{L_1L_2}\int_{x=0}^{L_1}\int_{y=0}^{L_2}
\cos\left(\displaystyle\frac{n {\rm \pi}
y}{L_2}\right) \varGamma(x,y,t_f) \,\mathrm{d}y
\,\mathrm{d}\kern0.7pt x , & \text{for } n>0, m=0,\\
\displaystyle\frac{2}{L_1L_2}\int_{x=0}^{L_1}\int_{y=0}^{L_2}
\cos\left(\displaystyle\frac{m {\rm \pi}
x}{L_1}\right) \varGamma(x,y,t_f) \,\mathrm{d}y
\,\mathrm{d}\kern0.7pt x , & \text{for } m>0, n=0,\\
\displaystyle\frac{4}{L_1L_2}\int_{x=0}^{L_1}\int_{y=0}^{L_2}
\cos\left(\displaystyle\frac{m {\rm \pi}
x}{L_1}\right)
\cos\left(\displaystyle\frac{n {\rm \pi}
y}{L_2}\right) \varGamma(x,y,t_f) \,\mathrm{d}y\,
\mathrm{d}\kern0.7pt x , & \text{for } m,n>0. \end{cases}
\end{align}
It is relevant to our arguments for using the approximation (2.30) to point out that the late-time mapping approximation (B5) is exactly the gradient of the scalar potential
\begin{equation} \phi_{cr} = \left(\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\frac{ \sigma_{mn}}{4 \omega_{mn}} \cos{(m{\rm \pi} x/L_1)}\cos{(n {\rm \pi}y/L_2)} \right) - \frac{ \sigma_{00}}{4 \omega_{00}}. \end{equation} At a large time 
$t=t_f$, we expect the three leading order terms for 
$\varGamma (x,y,t_f)$ to be
\begin{equation} \varGamma(x,y,t_f) \approx \bar{\varGamma} + \sigma_{01}\cos{\left(\frac{{\rm \pi} y}{L_2}\right)} + \sigma_{10}\cos{\left(\frac{{\rm \pi} x}{L_1}\right)}, \end{equation}
as higher-order modes are subject to exponential decay over a smaller time scale through (B3). (Taking 
$t_f=1047.8$ with 
$\delta =0.25$, for example, we found that the next largest coefficient 
$\sigma _{mn}$ was a factor of 
$O(10^{-3})$ smaller than either 
$\sigma _{01}$ or 
$\sigma _{10}$, justifying (B8).) If we write 
$\varGamma _a$ for the concentration at 
$(0,L_2)$ and 
$\varGamma _b$ for the concentration at 
$(L_1,L_2)$, then
\begin{equation} \varGamma_a \approx \bar{\varGamma} + \sigma_{10} - \sigma_{01}, \quad \varGamma_b \approx \bar{\varGamma} - \sigma_{10} - \sigma_{01} , \end{equation}which means
\begin{equation} \sigma_{10} \approx \frac{1}{2}\,(\varGamma_a-\varGamma_b), \quad \sigma_{01} \approx \frac{1}{2}\,(2\bar{\varGamma}-\varGamma_a-\varGamma_b) . \end{equation}
The leading-order terms for the correction to the mapping 
$\boldsymbol {X}_{cr}$ that maps particles from 
$t_f$ to 
$t\to \infty$ can therefore be approximated as
\begin{equation} \boldsymbol{X}_{cr} \approx \begin{pmatrix} \dfrac{ L_1(\varGamma_a-\varGamma_b)}{2{\rm \pi} \bar{\varGamma}} \sin{({\rm \pi} x/L_1)}\\ \dfrac{ L_2(2\bar{\varGamma}-\varGamma_a-\varGamma_b)}{2{\rm \pi} \bar{\varGamma}} \sin{( {\rm \pi}y/L_2)} \end{pmatrix}. \end{equation}
We run a given Eulerian particle-tracking simulation to a time 
$t_f$ where the leading-order amplitude of the correction is below a pair of chosen tolerance levels 
$[X_{tol},Y_{tol}]$, or in other words, until
\begin{equation} \left| \frac{ L_1(\varGamma_a-\varGamma_b)}{2{\rm \pi} \bar{\varGamma}} \right| < X_{tol}, \quad \left| \frac{ L_2(2\bar{\varGamma}-\varGamma_a-\varGamma_b)}{2{\rm \pi} \bar{\varGamma}} \right| < Y_{tol}. \end{equation}Appendix C. Behaviour of solutions to the nonlinear diffusion equation (2.1) near sharp corners
We wish to understand the behaviour of the two-dimensional nonlinear diffusion equation (2.1) subject to boundary conditions (2.2), near a sharp corner of a wedge-shaped domain. We introduce a polar coordinate system 
$(r,\theta )$ with the origin at the corner, with one boundary located at 
$\theta =0$, and the other at 
$\theta =\varPhi$. The surfactant concentration 
$\varGamma$ must neither diverge nor go to zero at the corner, so we look for expansions of the form
\begin{equation} \varGamma_{(m)} = A_{m,0}(\theta,t) + A_{m,1}(\theta,t)\,r^{a_{m,1}}+ A_{m,2}(\theta,t)\,r^{a_{m,2}}+\cdots, \end{equation}
where 
$0< a_{m,1}< a_{m,2}<\cdots$, and 
$m$ is an index as we anticipate multiple expansions, each indexed by a different value of 
$m$. The result will be a sum of asymptotic series, a technique used to derive a corner expansion for surfactant-induced flow in Mcnair, Jensen & Landel (Reference Mcnair, Jensen and Landel2022). Substituting an expansion for a given 
$m$ into (2.1), we obtain
\begin{align} &4\left( \frac{\partial A_{m,0}}{\partial t} + r^{a_{m,1}}\,\frac{\partial A_{{m,1}}}{\partial t} + r^{a_{m,2}}\,\frac{\partial A_{{m,2}}}{\partial t}+ \cdots \right)\nonumber\\ &\quad= \frac{1}{2r^2}\,\frac{\partial^2 A_{m,0}^2}{\partial \theta^2} + a_{m,1}^2A_{m,0}A_{{m,1}}r^{a_{m,1}-2} + r^{a_{m,1}-2}\,\frac{\partial^2}{\partial \theta^2} (A_{m,0} A_{{m,1}}) \nonumber\\ &\qquad + \frac{1}{2}\,r^{2a_{m,1}-2}\,\frac{\partial^2 A_{m,1}^2 }{\partial \theta^2} + A_{m,1}^2 (2a_{m,1})^2 r^{2a_{m,1}-2} + r^{a_{m,2}-2}\,\frac{\partial^2}{\partial \theta^2}(A_{m,0}A_{m,2})\nonumber\\ &\qquad + a_{m,2}^2A_{m,0}A_{m,2}r^{a_{m,2}-2}+ \cdots . \end{align}
As 
$a_{m,1}>0$, the leading-order equation is 
$\partial ^2 A_{m,0}^2/\partial \theta ^2 = 0$, which, when solved subject to boundary conditions 
$\partial A_{m,0}/\partial \theta =0$ (from (2.2)), gives 
$A_{m,0}=A_{m,0}(t)$. The balance at the next order is
\begin{equation} 4\,\frac{\partial A_{m,0}}{\partial t} = a_{m,1}^2A_{m,0}A_{m,1}r^{a_{m,1}-2} + r^{a_{m,1}-2}\,\frac{\partial^2}{\partial \theta^2} (A_{m,0} A_{m,1}). \end{equation}
One possible exponent, which we index with 
$m=0$, in the expansion is 
$a_{0,1}=2$, yielding an expansion driven by 
$\partial A_{0,0}/\partial t$,
\begin{equation} \varGamma_{(0)} = A_{0,0}(t) +\frac{\partial A_{0,0}/\partial t}{A_{0,0}}\,r^2+ O(r^4), \end{equation}
representing a purely radial flow that drives surfactant into or out of the corner, leading to changes in the corner concentration 
$A_{0,0}(t)$. The fact that the series goes in powers of 
$r^{2n}$ can be obtained by examining (C2) at the next order.
We now look for other possible expansions indexed by 
$m=1,2,3,\ldots$, and to avoid duplication of the primary flow contribution (C4), we set 
$A_{1,0}=A_{2,0}=\cdots =0$. (When 
$\varPhi ={\rm \pi} /2$, again assuming 
$a_{0,1}=2$, (C3) also possesses a homogeneous solution 
$A_{0,1}=f_{0,1}(t) \cos (2\theta )$, representing a stagnation point flow in the corner, with strength 
$f_{0,1}(t)$ determined by conditions far from the corner.) More generally, (C3) possesses homogeneous solutions satisfying, for 
$m=1,2, \ldots$,
\begin{equation} a_{m,1}^2A_{m,1} + \frac{\partial^2}{\partial \theta^2} ( A_{m,1}) = 0 . \end{equation}
This means that 
$A_{{m,1}} = f_{{m,1}}(t)\cos {(a_{m,1}\theta )}$, after applying 
$\partial A_{m,1}/\partial \theta =0$ at 
$\theta =0$. Applying 
$\partial A_{m,1}/\partial \theta =0$ at 
$\theta =\varPhi$ means that
\begin{equation} a_{m,1} = \frac{m{\rm \pi}}{\varPhi}, \end{equation}
for any integer 
$m$, and this provides an infinite set of possible exponents. (When 
$\varPhi ={\rm \pi} /2$, (C6) with 
$m=1$ recovers the stagnation point flow cited above with amplitude 
$f_{0,1}$.) By examining the next-order balance in (C2), it must also be the case that 
$a_{m,2} = m{\rm \pi} /\varPhi +2$, and so on. Therefore, the full solution is given by the sum of asymptotic series
\begin{equation} \varGamma = \varGamma_{(0)} +\varGamma_{(1)}+\varGamma_{(2)}+\varGamma_{(3)}+\cdots, \end{equation}
where 
$\varGamma _{(0)}$ is given by (C4), and
\begin{gather} \varGamma_{(1)} = f_{1,1}(t)\cos{({\rm \pi} \theta/\varPhi)}\,r^{{\rm \pi}/\varPhi}+ A_{1,2}(\theta,t)\,r^{{\rm \pi}/\varPhi+2}+ A_{1,3}(\theta,t)\,r^{{\rm \pi}/\varPhi+4}+ \cdots, \end{gather}
\begin{gather} \varGamma_{(2)} = f_{2,1}(t)\cos{(2{\rm \pi} \theta/\varPhi)}\,r^{2{\rm \pi}/\varPhi}+ A_{2,2}(\theta,t)\,r^{2{\rm \pi}/\varPhi+2}+ A_{2,3}(\theta,t)\,r^{2{\rm \pi}/\varPhi+4}+ \cdots, \end{gather}
\begin{gather} \varGamma_{(3)} = f_{3,1}(t)\cos{(3{\rm \pi} \theta/\varPhi)}\,r^{3{\rm \pi}/\varPhi}+ A_{3,2}(\theta,t)\,r^{3{\rm \pi}/\varPhi+2}+ A_{3,3}(\theta,t)\,r^{3{\rm \pi}/\varPhi+4}+ \cdots , \end{gather}and so on. The series can be summarised as
\begin{equation} \varGamma(r,\theta,t) = \sum_{n=0}^{\infty}\sum_{m=0}^{\infty}A_{m,n}(\theta,t)\,r^{m{\rm \pi}/\varPhi+2n}. \end{equation}
The velocity field is 
$\boldsymbol {u}_s = -\boldsymbol {\nabla }_{\boldsymbol {x}}\varGamma /4$, so as 
$r\to 0$, it is a combination of a radial and a stagnation point flow:
\begin{align} \boldsymbol{u}_s\approx \hat{\boldsymbol{r}}\left(2\,\frac{\partial A_{0,0}/\partial t}{A_{0,0}}\,r+\frac{\rm \pi}{4\varPhi}\,f_{1,1}(t)\cos\left(\frac{{\rm \pi}\theta}{\varPhi}\right) r^{({\rm \pi}/\varPhi)-1}\right)-\hat{\boldsymbol{\theta}}\,f_{1,1}(t)\, \frac{\rm \pi}{4\varPhi}\sin\left(\frac{{\rm \pi}\theta}{\varPhi}\right) r^{({\rm \pi}/\varPhi)-1}. \end{align}
For internal angles in a convex domain (when 
$\varPhi <{\rm \pi}$), the velocity is proportional to 
$r$ to a positive power, so it goes to zero at the corner. In particular, when 
$\varPhi ={\rm \pi} /2$, a particle on the boundary 
$\theta =0$ or 
$\theta ={\rm \pi} /2$ has an inward radial velocity bounded above by 
$Fr$ for some finite 
$F>0$. A particle starting at 
$r=r_0$ at 
$t=0$ will then lie in 
$r>r_0\exp (-F t )$ for 
$t>0$, never reaching the corner in finite time, supporting the use of (2.25a,b). For wedge angles 
$\varPhi <{\rm \pi} /2$, the radial flow is dominant as 
$r \rightarrow 0$, and (2.25a,b) again applies. For wedge angles satisfying 
$\varPhi \in ({\rm \pi} /2,{\rm \pi} )$, the stagnation point flow dominates as 
$r\rightarrow 0$. In this case, although the velocity field along a boundary vanishes as 
$r\rightarrow 0$, a sustained inward flow 
$\mathrm {d}r/\mathrm {d}t=-Fr^{({\rm \pi} /\varPhi )-1}$ for constant 
$F>0$ has the potential to drive particles to the corner in finite time, with 
$r=[F(t_0-t)(2-({\rm \pi} /\varPhi ))]^{1/(2-({\rm \pi} /\varPhi ))}$ as 
$t\rightarrow t_0$ for some 
$t_0$, calling into question the validity of (2.25a,b) in this case. For even larger wedge angles (
$\varPhi >{\rm \pi}$), regularity of the velocity field demands 
$f_{1,1}=0$. This highlights the influence of any smoothing over a length scale of order 
$\delta \ll 1$ that might be applied to computations of flows around such corners in non-convex domains, where velocities of magnitude 
$1/\delta ^{1-({\rm \pi} /\varPhi )}$ can be anticipated.
Appendix D. The numerical scheme to solve the Monge–Ampère equation
D.1. Outline of the scheme
We approximate (2.33) and (2.32) on an 
$(M+2)\times (N+2)$ rectangular grid. The points on the four boundaries we consider as known, as they can be expressed as functions of the interior points by the boundary conditions using one-sided second-order differences. This leaves us with an 
$M\times N$ grid of unknowns. We create an 
$MN\times 1$ vector of the values of the solution 
$\phi$ at these grid points by stacking the rows of the grid of unknowns into a column vector 
$\boldsymbol {\phi }$, starting from the top and proceeding downwards. So, for example, 
$\phi _{i+N}$ refers to the grid point directly below 
$\phi _i$, and 
$\phi _{i+1}$ refers to the grid point directly to the right of 
$\phi _i$ unless 
$i=pN$ for some integer 
$p$.
We approximate (2.33) using second-order differences at the 
$i$th grid point as
\begin{align} f_i(\boldsymbol{\phi})&= \frac{(\phi_{i+1}-2\phi_{i}+\phi_{i-1})(\phi_{i+N}-2\phi_{i}+\phi_{i-N})}{\Delta x^2\,\Delta y^2}\nonumber\\ &\quad - \left(\frac{-\phi_{i+N+1}-\phi_{i-N-1}+\phi_{i+N-1}+\phi_{i-N+1}}{4\,\Delta x\,\Delta y}\right)^2 -\mathcal{G}_{i} = 0, \end{align}
which we define to be the 
$i$th component of the vector 
$\boldsymbol {f}(\boldsymbol {\phi })$. Here, 
$\Delta x$ and 
$\Delta y$ refer to the grid spacing in the 
$x_0$ and 
$y_0$ coordinate directions, respectively, and we define
\begin{equation} \mathcal{G} = 1 - \beta_j \left(1 - \frac{\varGamma_0(x_0,y_0)}{\bar{\varGamma}} \right). \end{equation} The 
$k$th iteration of the Newton scheme for the vector 
$\boldsymbol {\phi }$ of the values of the function 
$\phi$ approximated at the grid points is
\begin{equation} \boldsymbol{\phi}^{k} = \boldsymbol{\phi}^{k-1} - \big(\boldsymbol{\nabla}_{\boldsymbol{\phi}}\boldsymbol{f}(\boldsymbol{\phi}^{k-1})\big)^{{-}1}\boldsymbol{f}(\boldsymbol{\phi}^{k-1}) , \end{equation}
where the Jacobian 
$\boldsymbol {\nabla }_{\boldsymbol {\phi }}\boldsymbol {f}(\boldsymbol {\phi })$ is the sparse matrix of derivatives of the components 
$f_i$ with respect to solution at each grid point 
$\phi _j$, so that the element 
$\partial f_i /\partial \phi _j$ is zero, except for nine diagonals given by
$$\begin{gather} \frac{\partial f_i}{\partial \phi_{i}} = \frac{2\left(4\phi_{i}-\phi_{i+1}-\phi_{i-1}-\phi_{i+N}-\phi_{i-N}\right) }{\Delta x^2\,\Delta y^2} , \end{gather}$$
$$\begin{gather}\frac{\partial f_i}{\partial \phi_{i+1}} = \frac{\phi_{i+N}-2\phi_{i}+\phi_{i-N}}{\Delta y^2\,\Delta x^2}, \end{gather}$$
$$\begin{gather}\frac{\partial f_i}{\partial \phi_{i-1}} = \frac{\phi_{i+N}-2\phi_{i}+\phi_{i-N}}{\Delta y^2\,\Delta x^2}, \end{gather}$$
$$\begin{gather}\frac{\partial f_i}{\partial \phi_{i+N}} = \frac{\phi_{i+1}-2\phi_{i}+\phi_{i-1}}{\Delta y^2\,\Delta x^2}, \end{gather}$$
$$\begin{gather}\frac{\partial f_i}{\partial \phi_{i-N}} = \frac{\phi_{i+1}-2\phi_{i}+\phi_{i-1}}{\Delta y^2\,\Delta x^2}, \end{gather}$$
$$\begin{gather}\frac{\partial f_i}{\partial \phi_{i+N+1}} = \frac{-\phi_{i+N+1}-\phi_{i-N-1}+\phi_{i+N-1}+\phi_{i-N+1}}{8\,\Delta x^2\,\Delta y^2}, \end{gather}$$
$$\begin{gather}\frac{\partial f_i}{\partial \phi_{i-N-1}} = \frac{-\phi_{i+N+1}-\phi_{i-N-1}+\phi_{i+N-1}+\phi_{i-N+1}}{8\,\Delta x^2\,\Delta y^2}, \end{gather}$$
$$\begin{gather}\frac{\partial f_i}{\partial \phi_{i+N-1}} ={-} \frac{-\phi_{i+N+1}-\phi_{i-N-1}+\phi_{i+N-1}+\phi_{i-N+1}}{8\,\Delta x^2\,\Delta y^2}, \end{gather}$$
$$\begin{gather}\frac{\partial f_i}{\partial \phi_{i-N+1}} ={-} \frac{-\phi_{i+N+1}-\phi_{i-N-1}+\phi_{i+N-1}+\phi_{i-N+1}}{8\,\Delta x^2\,\Delta y^2}. \end{gather}$$
The definitions (D1) and (D4) are valid for every interior point of the 
$M\times N$ grid of unknowns; however, for the values 
$i=1$ to 
$N$, 
$i=(M-1)N+1$ to 
$MN$, 
$(p-1)N+1$ and 
$pn$ for every integer 
$1< p< M$, the elements of 
$f$ and the Jacobian are slightly different to this as they incorporate the boundary conditions, but we omit them here for conciseness. We iterate (D3) until the sum of absolute differences of the components of 
$\boldsymbol {\phi }$ from one iteration to the next falls below a small tolerance value, which we choose to be 
$10^{-6}$.
We use this numerical method to perform the calculation using the method presented in § 2.3.2 for a relatively coarse grid until 
$\beta _J=1$, and then increase the refinement by interpolating the solution onto a more refined grid, and re-solving (2.33) and (2.32) with 
$\beta _j=1$ (a multi-grid method). We repeat this second process until we have a sufficiently refined solution.
D.2. Convergence of the numerical scheme
Figure 11 demonstrates the convergence of the numerical scheme presented in detail in § D.1 for the solution to the problem set out in § 2.3.2 with 
$\delta =0.25$, 
$0.05$ and 
$0.002$. For the purpose of illustration, we have picked a point 
$(x_0,y_0)=(5,5)$ and displayed the error between an approximation for the correct solution and solutions for 
$\phi (5,5)$, 
$X=\phi _{x_0}(5,5)$ and 
$Y=\phi _{y_0}(5,5)$ at various values for a uniformly spaced grid with 
$h=\Delta x=\Delta y$. The approximation for the correct solution is found by selecting 
$h\equiv h_{end}$ smaller than all of the other values of 
$h$ used, and using this as the grid spacing to solve for 
$\phi (5,5)\equiv \phi _{end}(5,5)$, 
$X(5,5)=\phi _{x_0}(5,5)\equiv X_{end}(5,5)$ and 
$Y(5,5)=\phi _{y_0}(5,5)\equiv y_{end}(5,5)$, for each value of 
$\delta$.
Figure 11. Graph showing the convergence of the finite-difference approximation to the solution of the Monge–Ampère equation (2.31). To illustrate convergence, we choose a point on the Lagrangian domain 
$(x_0,y_0)= (5,5)$ and calculate the solution for the three-deposit problem outlined in § 2.3.2 with 
$\delta =0.25$, 
$0.05$ and 
$0.002$ (each 
$\delta$ is assigned a different line style – see legend) for multiple values of a uniform grid spacing 
$h$ (data points are shown as circles). A set of values for each 
$\delta$, which we call 
$\phi _{end}(5,5)$, 
$X_{end}(5,5)$ and 
$Y_{end}(5,5)$, is calculated for a discretisation parameter 
$h_{end}$ that is smaller than the rest of the values of 
$h$ used. We use this solution as our approximation to the correct solution. The curves show the difference between 
$X(5,5)$ (blue), 
$Y(5,5)$ (orange) and 
$\phi (5,5)$ (purple), and the approximated correct solution, showing convergence at a rate of order 
$h^2$ (dotted black line).
The solutions for all three values in figure 11 converge approximately with slope 
$2$ on the log scale indicated by a dotted black line, which is expected for our second-order finite-difference scheme. The graphs are only approximately showing this convergence, as there are numerous sources of noise, e.g. the approximation that we have used for the correct solution. Also, the code stops iterating once the sum of the absolute difference of each grid point in a solution between one iteration and the next falls below a certain small tolerance value. The exact locations of the edges of the deposits do not coincide precisely with grid points, and this is a further source of error as the initial conditions (and therefore the right-hand side of (2.31)) for differently discretised grids are effectively slightly different. The solutions for 
$X$ and 
$Y$ also involve calculating a further derivative using a second-order finite-difference approximation, introducing further error. However, the dominant error is clearly the square of the discretisation parameter 
$h$ used in the approximation (2.31). Solutions with 
$\delta =0.002$ are significantly more expensive computationally, which is why the curves terminate at larger values of 
$h$ than the other curves. The fact that we have used a larger value of 
$h$ to calculate 
$\phi _{end}(5,5)$ and its derivatives for 
$\delta =0.002$ might also partly explain why some of these curves are noisier than for larger values of 
$\delta$.
Appendix E. Simulating Suminagashi patterns
We create a Suminagashi pattern using successive solutions of the Monge–Ampère equation. As described in the Introduction (§ 1), the creation of Suminagashi art patterns consists of repeated deposition of an ink–surfactant mixture onto a liquid surface, followed by the artist blowing on the surface between ink depositions. We simulate this procedure using the Monge–Ampère equation (2.31) to simulate the final locations of ink depositions in the spreading steps, and a divergence-free map to create the effects of the blowing steps. We show in (A11) how the rotational component of a surface shear stress can redistribute surfactant without changing its surface concentration.
The physical assumptions that we consider in this study for the Monge–Ampère method are appropriate for Suminagashi. Exploiting conditions (A8) and (A9): for a typical 1 cm layer of water with a surface tension reduction of 
$1\unicode{x2013} 10\,\textrm {g}\,\textrm {s}^{-2}$, we find that the Bond number is large and of the order of 
$10^{2}$–
$10^{3}$; surface diffusivity for most surfactants is of the order of 
$D^*= 10^{-10}$–
$10^{-9}\,\mathrm {m}^2\,\textrm {s}^{-1}$ (Chang & Franses Reference Chang and Franses1995), which leads to a very large surface Péclet number of the order of 
$10^{7}$–
$10^{9}$; furthermore, Suminagashi patterns can be transferred to paper within a minute, during which any boundary thickens under molecular (surface) diffusivity by less than 
${0.3}\,\mathrm {mm}$.
E.1. Spreading steps
Each spreading step starts from initial conditions (2.7) with 
$\mathcal {F}(x_0,y_0) = \mathcal {C}_q(x_1,y_1,0.5,1-\delta )$, and by imposing 
$L_1=13$, 
$L_2=11$, so that this is now a single deposit of exogenous surfactant spreading into endogenous surfactant. For the first spreading step, we set 
$\delta =0.01$, and solve (2.31) and (2.32). We store the final location of the boundary of this initial deposit, which we call 
$C_{1,(1)}$, where the first index labels each deposit edge, and the second index refers to how many total deposits have been released. Next, we set up a new problem by replacing 
$\delta$ by 
$\bar {\varGamma }$ from the previous problem, and we shift 
$(x_3,y_3)$ to a displaced position (using an algorithm described below as (E4)), and solve again (2.31) and (2.32), finding the edge of the new deposit, and the new locations of the edges of previously released deposits, repeating the above process 
$J$ times. This process is summarised below.
E.2. Summary of the algorithm for the spreading steps
For 
$j=1$ to 
$J$: solve (2.31) and (2.32) subject to initial conditions (2.7) with 
$\mathcal {F}(x_0,y_0) = \mathcal {C}_q(x_1,y_1,0.5,1-\delta )$, using the solution method outlined in § 2.3.2 with 
$L_1=13$, 
$L_2=11$, 
$\varGamma _2=\varGamma _3=\delta$, 
$r_1=0.5$, 
$r_2=0$, 
$r_3=0$, with 
$(x_1,y_1)=(x_{1,j},y_{1,j})$ and 
$\delta =\delta _j$, where this last quantity is found using
\begin{equation} \delta_{j} = \bar{\varGamma}_{j-1}, \quad \bar{\varGamma}_0 = 0.05, \end{equation}
and where at each step 
$\bar {\varGamma }_j$ is computed using (2.4). We find all curves at step 
$j$ by
\begin{equation} C_{i,(j+1)} = \boldsymbol{X}_{j}(C_{i,(j)} ) \quad \text{for} \ i=1,2,\ldots, j, \end{equation}
where 
$\boldsymbol {X}_j$ is the solution of (2.31) and (2.32) at step 
$j$, and one new evolving curve is introduced at each step 
$j$ by
\begin{equation} C_{j,(j)} = \{ (x,y)\,|\,(x-x_{1,j})^2+(\kern0.7pt y-y_{1,j})^2=0.5\}. \end{equation}The centres of the new deposits are chosen such that
\begin{equation} (x_{1,j+1},y_{1,j+1}) = (x_{c,j}-0.75,y_{c,j}-0.2), \quad \text{where } (x_{c,j},y_{c,j}) = \max_x{(C_{j,(j)}}), \end{equation}
until 
$\max _{x}(C_{j,(j)})>12.9$, and then
\begin{equation} (x_{1,j+1},y_{1,j+1}) = (x_{c,j}-0.05,y_{c,j}+0.75), \quad \text{where } (x_{c,j},y_{c,j}) = \min_y{(C_{j,(j)}}), \end{equation}with
\begin{equation} (x_{1,1},y_{1,1}) = (3.5,8.5). \end{equation}
We vary the locations of the centres of the new deposits, as we have observed Suminagashi artists to do this. The final picture is created by plotting 
$C_{i,(J+1)}$ for all 
$i=1$ to 
$J$.
E.3. Blowing steps
Between some of the spreading steps, we introduce blowing steps, which are computed by finding a divergence-free map 
$\boldsymbol {X}_b$ by time-stepping
\begin{equation}
\frac{\mathrm{d}}{\mathrm{d}t } \boldsymbol{X}_b =
\sum_{n_s=1}^{\infty}A_{n_s}\begin{pmatrix}
\dfrac{n_{s}{\rm \pi}}{11}\sin{\left(\dfrac{n_{s}{\rm \pi}
x}{13}\right)}\cos{\left(\dfrac{n_{s}{\rm \pi} y}{11} \right)}\\
- \dfrac{n_{s}{\rm \pi}}{13}\cos{\left(\dfrac{n_{s}{\rm \pi}
x}{13}\right)}\sin{\left(\dfrac{n_{s}{\rm \pi} y}{11} \right)}
\end{pmatrix}, \end{equation}
for some integer 
$n_{s}$, between 
$t=0$ and some 
$t_{s}$, where 
$A_{n_s}$ are the amplitudes of the modes, and where 
$\boldsymbol {X}_b=[x_0,y_0]$ at 
$t=0$. The resulting map 
$\boldsymbol {X}_b$ is divergence-free (such that concentrations do not change) and satisfies the required boundary conditions. The functions on the right-hand side of (E7) form a basis, such that any divergence-free map can be obtained by choosing a set of amplitudes 
$A_{n_s}$ and some choice of 
$t_s$. The new location of deposit edge 
$i$, after spreading step 
$j$, is given by 
$\boldsymbol {X}_b(C_{i,(j)})$.
E.4. Sequential depositions: the Suminagashi patterns
To illustrate the utility of the Monge–Ampère approximation, a solution of the Suminagashi algorithm is presented in figure 10. Here, for (E7), we choose 
$A_{ns}=0$ for every 
$n_s$ except 
$N_s=4$, and we choose 
$A_4=1$. We run 20 spreading steps, with a blowing step imposed after every four spreading steps with 
$t_s=0.15$. The creation of figure 10( f) is shown after four spreading steps, after four spreading steps and one blowing step, after 
$12$ spreading steps, after 
$12$ spreading steps and three blowing steps, and finally after 
$20$ spreading steps before the final blowing step. The final result is given in a monochrome colour scheme in figure 10(e).
Appendix F. Difference in the prediction between Eulerian particle tracking and Monge–Ampère for different values of $\delta$
Figure 12 shows a series of box and whisker plots of the normalised absolute error (2.34), measuring the Euclidean distance between predictions for final particle location between Monge–Ampère and Eulerian particle-tracking, for multiple values of 
$\delta$. Except for the solution with 
$\delta =0.002$, which has quadratic initial conditions (2.7) and (2.9), the rest of the solutions have cosine-shaped exogenous deposit profiles such that
\begin{align}
\mathcal{C}_c(\boldsymbol{x}_0;\boldsymbol{x}_c,r,\varGamma_{0,c}-\delta)
= \begin{cases}
\left(\dfrac{\varGamma_{0,c}-\delta}{2}\right)\left(1+\cos\left(
\dfrac{{\rm \pi}(|\boldsymbol{x}_0-\boldsymbol{x}_c|^2)}{r^2}
\right) \right), & |\boldsymbol{x}_0-\boldsymbol{x}_c|\leq
r, \\ 0, & |\boldsymbol{x}_0-\boldsymbol{x}_c| >
r,\end{cases} \end{align}
replaces 
$\mathcal {C}_q(\boldsymbol {x}_0;\boldsymbol {x}_c,r,\varGamma _{0,c}-\delta )$ in (2.7) and (2.9). The shape of the exogenous deposit profiles does not have a significant effect on results as shown in figure S7 of § S5 of the supplementary material, but we change the profile shape here for the sake of variety.
Figure 12. Graphs showing the normalised absolute error (2.34) for several different values of 
$\delta$. The initial conditions are (2.7) and (2.9), with 
$\mathcal {C}_q$ replaced by 
$\mathcal {C}_c$ defined in (F1) except for the first solution, which is taken from Eul[0.002] and MA[0.002] and has a quadratic initial profile. (a) The logarithm of the Euclidean distance between the predictions of final location, for each initial particle location from a 
$221\times 261$ grid, represented as a data point. The whiskers show 1.5 times the interquartile range above and below the quartiles, and particles outside this range are considered outliers and plotted as a cloud of points in blue, with the 90th percentile plotted as an orange line. (b) The same data as in (a), plotted in the form of a cumulative distribution function (CDF), with 
$\delta$ increasing in the direction of the arrow. Horizontal dashed lines indicate the median, the 75th percentile and the 90th percentile, respectively.
For moderate 
$\delta$, the error between different methods remains small, but it increases rapidly as 
$\delta$ tends to zero. However, for 
$\delta$ as small as 0.075, the normalised absolute error is below 0.01 for 75 % of particles, and even for 
$\delta =0.002$, the error is below 0.05 for 90 % of particles.
Open access
