DOI:
10.1039/C7SM02101A
(Paper)
Soft Matter, 2018,
14, 18561869
Feedback control of photoresponsive fluid interfaces†
Received
25th October 2017
, Accepted 23rd December 2017
First published on 19th February 2018
Photoresponsive surfactants provide a unique microfluidic driving mechanism. Since they switch between two molecular shapes under illumination and thereby affect surface tension of fluid interfaces, Marangoni flow along the interface occurs. To describe the dynamics of the surfactant mixture at a planar interface, we formulate diffusion–advection–reaction equations for both surfactant densities. They also include adsorption from and desorption into the neighboring fluids and photoisomerization by light. We then study how the interface responds when illuminated by spots of light. Switching on a single light spot, the density of the switched surfactant spreads in time and assumes an exponentially decaying profile in steady state. Simultaneously, the induced radial Marangoni flow reverses its flow direction from inward to outward. We use this feature to set up specific feedback rules, which couple the advection velocities sensed at the light spots to their intensities. As a result two neighboring spots switch on and off alternately. Extending the feedback rule to light spots arranged on the vertices of regular polygons, we observe periodic switching patterns for evensided polygons, where two sets of nextnearest neighbors alternate with each other. A triangle and pentagon also show regular oscillations, while heptagon and nonagon exhibit irregular oscillations due to frustration. While our findings are specific to the chosen set of parameters, they show how complex patterns at photoresponsive fluid interfaces emerge from simple feedback coupling.
1 Introduction
New microfluidic devices require a better understanding of fluids and flow control on a microscopic scale. For example, advances in digital and inertial microfluidics have resulted in useful biochemical reactors^{1} and precise sorting for micronsized objects.^{2} These devices often rely on specific geometries to control fluid flow.^{3–6} But flow can also result from gradients in surface tension through a phenomenon known as the Marangoni effect. Gradients of surface tension at the interface of fluids result from curvature^{7,8} or temperature^{9,10} gradients, or from gradients in the concentration^{11–14} of surfactants, i.e. molecules which attach to fluid interfaces and reduce surface tension. Surfactants containing one or more stereoisomers, like azobenzene,^{15} spiropyran,^{16} or stilbene^{17} are photoresponsive: they reversibly change shape under illumination with specific wavelengths of light and thereby interact differently with bordering fluids and neighboring surfactants.^{18} Consequently, their different shapes affect surface tension differently and lead to Marangoni advection between illuminated and dark regions of the interface.^{19,20} Through this mechanism, photoresponsive surfactants offer a unique way to control fluid flow by shining static or changing patterns of light on an interface.^{21,22}
Various applications for photoresponsive surfactants have been reported: they have been used to create responsive foams,^{23,24} fluids with tunable viscosity,^{25} and photoresponsive polymer brushes.^{26} Beyond tunable materials, they allow to move and position individual surfactantcovered droplets,^{27} to turn droplets into selfpropelled particles,^{20} and to pump fluid near a planar water surface.^{19} The transport mechanisms in all three examples are based on the interplay of two or more processes determining the dynamics of photoresponsive surfactants at fluid interfaces.
Here, we consider the four processes displayed schematically in Fig. 1: firstly, surfactants diffuse at the interface thermally and drift down chemical gradients. Secondly, surfactants at the interface are advected by Marangoni flow set up by gradients in surface tension in the bordering fluids. Thirdly, surfactants occasionally desorb from the interface into one of the fluids or adsorb to the interface from a fluid. And fourthly, surfactants photoisomerize or switch molecular shape, when they are exposed to light of a specific wavelength. All four processes are rationalized in two diffusion–advection–reaction equations for the surfactant densities. Earlier theoretical descriptions^{28–31} have not considered photoisomerization, or in the case of ref. 20 considered it only in the highintensity limit. However, because these processes happen concurrently, they couple to each other and produce complex behaviors: for example, photoisomerization affects surface tension, which induces advection, and also leads to concentration gradients and diffusion. By integrating all processes in our theory, we are able to investigate these complex couplings.

 Fig. 1 Schematics of all four processes governing surfactant motion at a fluid–fluid interface.  
While previous research has focused on manipulating droplets,^{20,27,32} we use the relative simplicity of a planar interface to generate and study more complex flow and light patterns in theory. In particular, we aim at demonstrating that the method of feedback control,^{33} where some inherent dynamic quantity is sensed and coupled back to the system, can be used to design fluid flow and induce novel dynamic patterns.
In this article, we specifically study the effect of pointlike spots of light on a planar fluid–fluid interface with photoresponsive surfactants for a specific set of realistic parameters. Because these surfactants usually respond to specific wavelengths of light, one can easily imagine a laser beam or similarly focused light source generating these spots.^{21}
When a single light spot is switched on, it generates a dynamic response where the density of the switched surfactant spreads in time and assumes an exponentially decaying profile in steady state. Simultaneously, the induced radial Marangoni flow reverses its flow direction from inward to outward. We use this feature to set up a specific feedback rule based on advection velocities sensed at the light spots. They couple back to the on/off states of the light spots and thereby initiate selfgenerated oscillations of two neighboring spots. We generalize our approach to regular polygons. While we mostly observe regular oscillatory patterns, seven and ninesided polygons exhibit irregular oscillations due to frustration. The simple on/off switching of the light spots is similar to thermostat control, which typically leads to oscillatory patterns^{34,35} and has been applied to microfluidic problems before.^{36,37}
The article is organized as follows: we describe our theory in Section 2 and the numerical method to solve the coupled dynamic equations in Section 3. In Section 4.1 we investigate the response to a single light spot switching on and off. Section 4.2 presents two light spots coupled by feedback control and Section 4.3 the extension to regular polygons with up to 10 edges. We close with a summary and conclusions in Section 5.
2 Theory
There is a long tradition of describing surfactant covered interfaces with diffusion–advection–reaction equations^{28–31} governing the coupled dynamics of surfactant densities ρ_{q}, 
 (1) 
Here, the diffusive currents j^{D}_{q} result from thermal motion and local interactions of the surfactants, the advection velocity v^{A} from drift due to the Marangoni effect, and the reaction terms S^{ads}_{q} describe ad and desorption to and from a fluid–fluid interface, respectively. In addition, we introduce reaction terms S^{iso}_{q} to describe photoisomerization. Our first objective is to close this system of partial differential equations by expressing all four r.h.s. terms of eqn (1) as functions of the respective densities of isomers A and B, which we denote ρ_{A} and ρ_{B}. In the following we describe all terms in more detail, formulate the equations in reduced units, and give relevant values for all material parameters.
2.1 Free energy
Motivated by the Flory–Huggins theory,^{38,39} we start from a coarsegrained Helmholtz free energy F of a compressible mixture of two components, A and B, with meanfield interactions. We define F as the integral over a free energy density f on the interface, such that . The first term of f contains entropic contributions of each component (see e.g. Chapter 3.1 in ref. 40), the second term describes meanfield interactions between the surfactants of the same or different type, and the final term is the intrinsic energy density σ_{0} of a clean interface without surfactants: 
 (2) 
Entropic terms are proportional to temperature k_{B}T and contain the thermal de Broglie wavelength λ, which does not affect dynamics in our case. Interaction terms are characterized by three energy constants E_{A}, E_{AB} and E_{B} and the saturation density ρ_{∞} of a fully covered interface.
In Appendix B we describe how eqn (2) relates to Flory–Huggins theory.
2.2 Diffusive currents
To formulate diffusive currents, we apply dynamical density functional theory (DDFT) as originally derived by Marconi and Tarazona^{41} for colloids and generalized to binary mixtures by Archer.^{42} The original derivation builds on a variant of the Nparticle Smulochowski equation for a mixture of particles of type q and density ρ_{q}.
By following the standard procedure, Archer arrives at a set of continuity equations. They are coupled via current densities j_{q}, which depend on functional derivatives of F:

 (3) 
We calculate them for each surfactant species and assume identical friction coefficients
γ_{A} =
γ_{B} =
γ. Using the free energy density
f from
eqn (2) and introducing the diffusion constant
D =
k_{B}T/
γ, we arrive at

 (4) 

 (5) 
Linear terms represent conventional thermal diffusion and nonlinear terms result from interactions between surfactants.
2.3 Advection: Marangoni effect
Another important contribution to the current densities is advection due to Marangoni flow initiated by gradients in surface tension. Here, we derive Marangoni currents for a planar fluid–fluid interface.
Surface tension σ is the thermodynamic conjugate variable to surface area, i.e. σ = ∂F/∂A. In Appendix A we show that for any free energy of the form

 (6) 
Like
f,
σ is based on the usual adiabatic approximation of local equilibrium. We substitute
f from
eqn (2) and arrive at

 (7) 
Intuitively, the first two terms mean any surfactant lowers surface tension equally but the third term takes into account their different interaction energies. For example, repulsive interactions (
E_{i} > 0) stabilize the interface by lowering surface tension.
We assume small Reynolds numbers and therefore need to solve the Stokes equations for the flow field of an incompressible fluid characterized by pressure p and velocity v = (v_{x},v_{y},v_{z}):

∇_{r}[η(z)∇_{r}v] − ∇_{r}p = −δ(z)∇_{x}σ and ∇_{r}·v = 0,  (8) 
where the viscosities of the two bounding fluids are

η(z > 0) = η_{+} and η(z < 0) = η_{−}.  (9) 
The interface is at
z = 0 and the gradient of surface tension
σ acts on both fluids like an external force. The vector
r = (
x,
y,
z) refers to any position in three dimensions, while
x = (
x,
y,0) is restricted to positions on the interface. Similarly, we define gradient operators ∇
_{r} = (∂
_{x},∂
_{y},∂
_{z}) and ∇
_{x} = (∂
_{x},∂
_{y},0). We do not allow flow through the interface and impose vanishing velocity and pressure far from the interface:
v_{z}(z = 0) = 0, v(z → ∞) = 0 and p(z → ∞) = 0 
The solution of eqn (8) can be computed using the modified Oseen tensor O(r) of Jones et al. to treat external forces acting along a fluid–fluid interface:^{43,44}

 (10) 
The tensor depends only on the average viscosity
= (
η_{+} +
η_{−})/2. The velocity field
v initiated by an arbitrary surface tension profile
σ is

 (11) 
We apply the product rule, Gauss's theorem and show in Appendix C that the resulting ring integral vanishes for a constant value
σ =
σ_{eq} at the boundary.
‡ Thus, we arrive at the velocity field

 (12) 
where the velocity response function
g(
r) is a threedimensional vector field, which decays like 
r
^{−2}:

 (13) 
Fig. 2 displays g in a plane perpendicular to the interface and containing the origin.

 Fig. 2 Fluid velocity field initiated by a positive pointlike inhomogeneity of the surface tension on a horizontal interface (green) separating two fluids. Arrows indicate direction and background brightness indicates speed. Close to the interface fluid is pulled towards the origin and pushed out along the normal.  
For our numerical solutions of eqn (1) we only need v at the interface, i.e. at z = 0. There, v(x) can be rewritten as a gradient field with δσ(x) = σ(x) − σ_{eq}:

 (14) 
This form of
v requires less computational effort than
eqn (11) and (12) because it minimizes the number of scalar convolutions we need to perform. Note that the integral only converges if
δσ(
x) decays like 
x
^{−ε},
ε > 1.
As a sidenote, the integral in eqn (14) is a Riesz potential,^{45} a form of the inverse half Laplacian (−Δ_{x})^{−1/2} operating on δσ. The operator is defined by its Fourier transform, such that

 (15) 
where
denotes the twodimensional Fourier transform, and
^{−1} its inverse.
2.4 Ad and desorption
We use the Frumkin approach to adsorption, which describes adsorbed particles as a real gas and desorbed particles as an ideal gas.^{46} Because it takes into account surface coverage and local interactions on the interface, it is an extension of the Langmuir adsorption model.^{47} Adsorption processes add two nonconservative terms to our differential eqn (1): one for adsorption, which is proportional to local free space on the interface, and one for desorption, which is proportional to the local density of each component and contains an Arrhenius factor to account for the changing surface tension or energy. Thus, the reaction terms in eqn (1) read 
 (16) 
where we also introduced adsorption and desorption rates k^{a}_{q} and k^{d}_{q} for each species. The Arrhenius energy barrier for removing particles is given by the corresponding change in surface tension, i.e. ∂_{ρq}σ.^{48} The adsorption rate k^{a}_{q} depends on the concentration of surfactants in the fluids and is assumed to be tunable in experiment. In thermal equilibrium the reaction term vanishes and one obtains the equilibrium densities ρ^{eq}_{q} by solving 
 (17) 
for given rate ratios k^{a}_{A}/k^{d}_{A} and k^{a}_{B}/k^{d}_{B}.
2.5 Photoisomerization
We describe photoisomerization using rate equations with local isomerization rates k_{i→j} in each direction: 
S^{iso}_{A}(ρ_{A},ρ_{B}) = k_{B→A}ρ_{B} − k_{A→B}ρ_{A}  (18) 

S^{iso}_{B}(ρ_{A},ρ_{B}) = k_{A→B}ρ_{A} − k_{B→A}ρ_{B}  (19) 
For our specific examples in Section 4 we only allow for photoinduced isomerization from A to B surfactants and therefore set k_{B→A} = 0. Note that total surfactant density is conserved during photoisomerization, which is quantified by S^{iso}_{A} + S^{iso}_{B} = 0.
The isomerization rates in our examples are spatially inhomogeneous. For example, a set of N pointlike light spots at positions x_{l} is described by a superposition of Diracδ functions with isomerization functions K_{l}(t):

 (20) 
Each
K_{l} is proportional to the intensity of its light spot
l. Thus, when the light spot is turned on,
K_{l} switches from 0 to a constant value
K_{A→B} specified in
Table 1.
Table 1 Our estimates for various material and photoisomerization parameters. The values for k^{a}_{A/B} are calculated using eqn (17)
Parameter 
Value 
E
_{A}

1.5k_{B}T 
E
_{AB}

1.3k_{B}T 
E
_{B}

1.1k_{B}T 
ρ
^{eq}_{A}

0.5ρ_{∞} 
ρ
^{eq}_{B}

0.0ρ_{∞} 
Mg 
1.2 × 10^{5} 
k
^{d}_{B}

1k^{d}_{A} 
k
^{a}_{A}

2.12k^{d}_{A} 
k
^{a}_{B}

1.73k^{d}_{A} 
K
_{A→B}

10D 
2.6 Closed and reduced equations
We nondimensionalize all quantities using k_{B}T for energies, desorption rate k^{d}_{A} for inverse time, saturation density ρ_{∞} for densities, and for length scales. Note that L is the length an A surfactant diffuses until it desorbs from the interface. In Appendix 10 we explain further the origin of these characteristic quantities, in particular, the rescaling of length and time.
The Marangoni number quantifies the importance of Marangoni flow by comparing the diffusion time for the characteristic length L, t_{diff} = L^{2}/D, to the time for advection by Marangoni flow, t_{Mg} = L/v_{Mg} = LΔσ/. Thus,

 (21) 
To arrive at the last expression, we set Δσ = k_{B}Tρ_{∞} and used the definition for L.
Nondimensionalizing eqn (1), (3), (14), (16) and (19) leads to reduced equations for each component q:

 (22) 
with

 (23) 
and the signs − and + of the photoisomerization term refer to the A and B surfactants, respectively.
Our numerical solutions are calculated on a square region of the interface, where we solve eqn (22) with open, Dirichletlike boundary conditions: we set the densities ρ_{q} to their equilibrium values ρ^{eq}_{q} at the boundary but also outside , so that

δσ(x ∉ ) = σ(ρ^{eq}_{A},ρ^{eq}_{B}) − σ_{eq} = 0.  (24) 
Thereby, the region of integration of the integral in
eqn (14) is restricted to
.
2.7 Parameters
The reduced equations contain parameters which are determined by material properties. Our estimates for these parameters are collected in Table 1. These estimates are not meant to replicate a specific combination of fluids and surfactants but rather to model a realistic system, which exhibits stable mixtures of surfactant species on the fluid interface.
The values for the interaction energies E_{A}, E_{B}, and E_{AB} are based on data for surface tensions of surfactant mixtures collected by Chevallier et al. for AzoTAB at a water–air interface^{19} assuming E_{AB} = (E_{A} + E_{B})/2 (for details see Appendix D). Our choice of Mg = 1.2 × 10^{5} is also based on their data, assuming D = 900 μm^{2} s^{−1}, ρ_{∞} = 5 nm^{−2}, k^{d}_{A} = 6.62 s^{−1}, and T = 300 K, and a larger value of the average viscosity = 2 mPa s.
The equilibrium surface densities ρ^{eq}_{A} and ρ^{eq}_{B} are assumed to be fully adjustable by changing the partial pressure of surfactants in the fluids, and for simplicity, we set ρ^{eq}_{B} = 0 and ρ^{eq}_{A} = 0.5ρ_{∞}. For the same reason we assume identical desorption rates k^{d}_{B} = k^{d}_{A} for both isomers.§ For reference, the aforementioned material parameters give a characteristic adsorption–diffusion length L ≈ 12 μm.
Photoisomerization is assumed to be fully tunable by adjusting laser intensity and focus. Throughout our study we set all photoisomerization rates using the value for K_{A→B} in Table 1 unless explicitely stated.
3 Numerical solver
We solve eqn (22) numerically using a combination of established algorithms implemented in the Julia programming language.^{49} There are three major steps, which can be treated separately from each other: integration in time, discretization in space and interpolation in space. We discuss all three in the following paragraphs. Readers not interested in numerical details can skip this section.
We use a finite volume discretization on a uniform grid. Accordingly, we tile the interface using roughly 56000 to 127000 hexagons (depending on the problem) and approximate all functions f by their values f_{i} at the center x_{i} of each cell i (similar to ref. 50). If we know f exactly, we approximate it by the spatial average over each cell _{i} with area A. For example, the photoisomerization field k(x) = Kδ(x − ) for a single light spot is discretized by

 (25) 
The main advantage of finite volume methods is that they reliably conserve particle number in discretized continuum equations. To discretize the diffusion term, we use the conventional central difference method, which is 2nd order accurate (see e.g.ref. 51). To treat the advective term correctly, we need an interpolation scheme to determine density values at the boundaries between cells. Here, we use an upwinddifference method developed by B. Koren (see Chapter 5 in ref. 52) which is 2nd to 3rd order accurate (depending on smoothness).¶ It is characterized by its fluxlimiting function ϕ:

 (26) 
where
r is the upwind gradient ratio between cells
i and
i + 1 of the advected field
ρ_{q}:

 (27) 
The index i + 1/2 indicates the boundary between cells i and i + 1. For a complete description of the method we refer to ref. 52.
Eqn (14) contains a nonlocal term which is a convolution of surface tension and x^{−1}. Rather than performing the direct convolution, which would need an (N^{2})operation on a mesh with N cells, we compute the integral in Fourier space. The fast Fourier transform (FFT) is an (NlogN)operation on orthogonal meshes. We use a regridding scheme to switch between hexagonal and rectangular cells and thereby take advantage of FFT (see Appendix F).
We integrate in time using the adaptive Runge–Kutta–Fehlberg (RKF45) method.^{55} RKF45 calculates 4th and 5thorder approximations of trajectories, estimates the integration error from their distance, and adapts its time step to keep it below a threshold. Integration proceeds using the 5thorder approximation.
In our case, a trajectory is characterized at every point by the vector ρ, which contains both densities, ρ = (ρ_{A},ρ_{B}), and exists at every point on the interface. The RKF45 error function should be a metric on the ρ_{q} function space. We base our metric on the L^{1}norm because for any physically valid ρ_{q} it simply gives the particle number N_{q}, which is dimensionless and necessarily exists. The L^{1} norm also has a straightforward counterpart for discretized functions, the Manhattan vector norm. Because each component ρ_{q} should fulfill the threshold criterion of the RKF45 method independently, we take the maximum of their norms:

 (28) 
The integration error e is simply the norm of the pointwise difference between two states ρ and :

e[ρ,] = ‖ρ(x) − (x)‖,  (29) 
where
ρ(
x),
(
x) refer to the 5th and 4thorder approximations, respectively. All data were produced using a tolance of
e ≤ 10
^{−6}.
4 Results
4.1 Single light spot
We start by perturbing the interface with a single pointlike spot of light. Specifically, we study how the interface responds to a spot switching on and off as well as how it approaches steady state. Our observations will form the basis for coupling the dynamics of multiple light spots by feedback control in Section 4.2.
We solve eqn (22) on a square region of area 20L × 20L with a single light spot placed at its center. Initially, only A surfactants are present. They are homogeneously distributed and in adsorption equilibrium with the bordering fluids. Then, the light spot is instantaneously turned on and switches surfactants from A to B. The resulting radial profiles for density ρ_{B}, radial velocity v_{r}, and surface tension δσ relative to their equilibrium values are displayed in Fig. 3 for various times t.

 Fig. 3 Radial profiles at various times t after a light spot at r = 0 is switched on: (a) surfactant density ρ_{B}, (b) radial velocity v_{r}, and (c) surface tension δσ relative to their equilibrium values. Inset: Blowup to show δσ at the three latest times. The solid line in each diagram shows the steady state nonequilibrium profile. Emblem in (a): the sensor senses v_{r} and ρ_{B} at a distance r from the light spot.  
As surfactants switch from A to B, type B surfactants start to diffuse radially away from the light spot and desorb (a). Simultaneously, a gradient in surface tension forms (c) and starts to drive a Marangoni current (b): at first toward the spot and later away from it. This is illustrated in Fig. 4, where we plot the radial velocity at distance r = 5L from the spot: it continuously changes from pointing toward the spot (negative) to pointing away from it (positive) at t ≈ 0.5k^{d}_{A}.

 Fig. 4 Radial velocity v_{r} at distance r = 5L from the light spot plotted versus time t for light switched on (solid, blue) and switched off in the nonequilibrium steady state (dashdotted, green). Emblem: a sensor measures v_{r} initiated by the light spot and feeds it back to the light source. Inset: Radial velocity versus time, when the light spot is switched on at t = 0 and switched off when v_{r} crosses the zero line (dashed, purple).  
Eventually, the surfactant densities approach a nonequilibrium steady state in which adsorption, diffusion, advection, and photoisomerization balance each other. In steady state, currents at the interface are nonzero because new type B surfactants continue to diffuse and drift away from the light spot until they desorb and new type A surfactants continue to replace switched surfactants by adsorption at the spot. The radial density profile ρ_{B}(r) in Fig. 3(a) decays like e^{−κr} with κ ∝ L as r → ∞.
An interesting aspect of the current system is the reversal of the direction of the radial velocity. Initially, Marangoni flow towards the light spot (v_{r} < 0) acts against the diffusive current of B surfactants, which in active emulsions can cause demixing (see ref. 31). However, here the Marangoni current ultimately aligns with diffusion of type B surfactants (v_{r} > 0), while reaching steady state. The directional change can be traced to an excess in the total surfactant density ρ = ρ_{A} + ρ_{B} > ρ^{eq}, which develops with time in an increasing region around the light spot (not shown). Because δσ ≈ −k_{B}T(ρ − ρ^{eq}) in leading order, the surface tension becomes negative [see Fig. 3(c)] and the resulting positive slope determines the outward direction of v.
In Fig. 4 we display the time evolution of the radial velocity v_{r} at a specific distance r, i.e. the response of the system to switching the light spot on. After reaching steady state, light is switched off and v_{r} relaxes towards zero. We now introduce a feedback control step in our system: when v_{r} changes sign, we switch off the light spot. The resulting time evolution of v_{r} is plotted in the inset. It gives an idea how we implement feedback control in Section 4.2.
For comparison, we include radial profiles in Fig. 5 for the case that the B surfactants desorb much faster than the A type (k^{d}_{B} = 300k^{d}_{A}) as observed in ref. 19. Now the surface tension is always positive. It has a negative slope and the Marangoni current keeps its initial direction towards the light spot in steady state. We attribute this simpler behavior to the increased desorption rate of type B, which prevents these surfactants from spreading across the interface. As a result, diffusion becomes irrelevant and steady state is dominated by Marangoni currents driving type A surfactants to the spot to replenish switched and desorbed surfactants there.

 Fig. 5 Radial profiles at various times t after a light spot at r = 0 is switched on for Mg = 6 × 10^{5} and k^{d}_{B} = 300k^{d}_{A}: (a) surface tension δσ relative to equilibrium value (inset: surfactant density ρ_{B}), and (b) radial velocity v_{r}. The solid line in each diagram shows the steadystate nonequilibrium profile.  
Adsorption and desorption rates can be influenced by changing the length of a surfactant's tail group, so that, for example, both isomers of C_{4}AzoOC_{6}E_{2} have almost identical adsorption behavior.^{18} In the following, we concentrate on the more interesting case k^{d}_{B} = k^{d}_{A}.
4.2 Two light spots coupled by feedback control
The main idea is now that one senses the flow velocities at each light spot and specifies a feedback rule that couples back to the light spots and determines for each of them if they are switched on or off [see emblem in Fig. 6(a)]. This means the system organizes its dynamic response by itself.

 Fig. 6 Example with two coupled light spots at a distance 10L from each other: (a) schematic of two light spots with their connected velocity sensors. (b) Velocities v_{1} and v_{2} sensed at the locations of the respective light spots 1 and 2 plotted over time t. (c) Dynamic fingerprints of two light spots: colorcoded time series for light spots (on: yellow/light, off: blue/dark) and velocities v_{1} and v_{2} (advection towards other spot: red, away from it: blue, darkness of shading indicates speed).  
For two light spots a simple rule is that always the spot at which the velocity is largest is on and the other spot is off. To rationalize this rule, we consider two light spots at respective positions x_{1} and x_{2} with distance Δx = x_{1} − x_{2} [see schematic in Fig. 6(a)]. In this setup the advection velocity v(x) at one spot is always directed towards or away from the other light spot due to symmetry. Accordingly, we define the advection velocity v_{l} at one spot l to be positive, when it points away from the other spot and to be negative in the opposite case:

 (30) 
Only the spot with the larger v_{l} is switched on or equivalently, the isomerization function K_{l}(t) of each spot, which we introduced in eqn (20), is

K_{l}(t) = K_{A→B}θ(v_{l}(t) −v_{j}(t)), j ≠ l.  (31) 
Here,
θ is the Heaviside step function and
K_{A→B} > 0 is the isomerization constant while a light spot is switched on. In the following, we always prepare the system in its equilibrium state and start by switching light spot 1 on at time
t = 0.
Movie M01 in the ESI† shows how the density of the type B surfactant and the Marangoni flow velocities vary in time for two light spots with distance Δx = 10L. In Fig. 6(b) we display the velocities at the two spots, which govern the feedback dynamics. As soon as the velocity v_{2} created by spot 1 at the location of spot 2 becomes positive, spot 1 is switched off and the active spot 2 generates a negative v_{1}. Over time the coupled light spots quickly establish regular oscillations in v_{1} and v_{2}, which occur in antiphase. In Fig. 6(c) we represent the on/off state of the light spots and the velocities v_{1}, v_{2} in dynamic fingerprints, which we will use further to show the dynamic response of the system.
We observed the same switching behaviour of the light spots for a large set of distances and isomerization constants (proportional to light intensity). The switching frequency f plotted in Fig. 7 hardly depends on K_{A→B} in the studied parameter range of more than three decades. Furthermore, increasing distance Δx by more than one decade reduces f by less than a factor of two. Overall, f is of the order of k^{d}_{A}, which suggests the spots are coupled through a diffusion–desorption process. An animation of the two alternating light spots is included as movie M02 in the ESI.†

 Fig. 7 Colorcoded switching frequency f of two coupled light spots plotted versus distance Δx and isomerization constant K_{A→B}.  
4.3 Coupled light spots on regular polygons
Now, we generalize the previous situation: instead of two light spots we study N spots arranged on the vertices of a regular polygon. All polygons have the same circumscribed circle with diameter 10L and the density and flow fields are computed on a square grid with a side length of 30L. The feedback control rule is now formulated as follows. Because the velocity at each spot does not have a predefined direction as in the case of 2 spots, we project it onto the outward radial direction of the circumscribed circle as shown in Fig. 8 and arrive at a scalar velocity v_{l} for each spot l: 
v_{l} = v(x_{l})·e_{l}.  (32) 

 Fig. 8 Schematic of 6 light spots positioned on the vertices of a regular hexagon, each with a radial projection vector e_{l} pointing away from the polygon's center.  
Each individual spot l is only kept switched on while v_{l} is larger than the average of all v_{l}. Equivalently, we introduce the isomerization function for spot l:

K_{l}(t) = K_{A→B}θ(v_{l}(t) − (t))  (33) 
with mean velocity
and the Heaviside step function
θ. If
N = 2 the condition simplifies to our feedback rule from Section 4.2. As before, we prepare the system in its equilibrium state but, unlike before, we linearly increase the intensity of light spot 1 from time
t = 0 to
t = 0.1
k^{d}_{A}. This initial condition may lead to some symmetry breaking in the observed dynamic patterns.
In the following we discuss the observed dynamics of the light spots using the dynamic fingerprints introduced in Fig. 6(c). Neighboring rows correspond to neighboring spots in the polygon. For all polygons animations of the switching light spots are included as Movies M03 to M10 in the ESI.†
We first concentrate on the polygons with an even number of edges. There is a strong anticorrelation between the on/off states of neighboring spots reminiscent of the oscillating state of two light spots. So after some possible transient regime two sets of nonneighboring spots alternate with each other. They create a checkerboard pattern in the dynamic fingerprint of the light spots as illustrated for N = 8 in Fig. 9. Polygons with varying N differ in how the regular pattern is reached in time. The square has a transient phase where the radial velocities at two opposing spots (1 and 3) vary much more than at the other two spots and the dynamics does not follow a clear pattern before settling in two pairs of alternating spots. The hexagon immediately enters such a regular pattern with equal variation in all velocities. The octagon (see Fig. 9) exhibits two transient phases: at first, a triple and a quintuple of spots alternate followed by an irregular transition phase into the stable oscillatory pattern of two alternating quadruples. In the triple–quintuple phase the velocities vary more at triple spots than at quintuple spots, which seems intuitive.

 Fig. 9 Dynamic fingerprints for an octagon of light spots. Top: ‘on’ and ‘off’ states of the light spots versus time indicated by yellow and blue, respectively. Each row represents one light spot numbered 1 to 8 from the bottom to the top. Bottom: Radial velocity at each light spot plotted versus time with positive and negative values indicating outward and inward direction, respectively. White and black bars in the bottom right corners represent the characteristic desorption time, i.e. 1/k^{d}_{A}.  
In polygons with an odd number of edges, always two neighboring light spots are frustrated since they cannot oscillate relative to each other. So, we expect more interesting dynamics. The equilateral triangle of light spots (not shown) immediately settles into an oscillation pattern, where a single spot (spot 1) and the pair of the other spots alternate, since by symmetry the latter switch simultaneously. The single spot is determined by the initial condition and there is a larger variation in radial velocity at spot 1, whereas the velocities at spots 2 and 3 are similar. The pentagon (Fig. 10, left) behaves like the triangle: it immediately reaches an oscillatory pattern, in which a pair and a triple of spots alternate. Here, the radial velocity at the pair spots varies more strongly than at the triple spots. The heptagon (Fig. 10, center) is the first polygon with a longlived irregular switching pattern, which does not reach a regular oscillation state. Notably, one observes some shorter on or off periods of spots, which propagate along neighborneighbor connections. Finally, the nonagon (Fig. 10, right) shows an irregular pattern as well and behaves similar to the heptagon, however, the variations in the velocities are not as pronounced. Thus, our investigations show that higher oddnumbered polygons settle into irregular patterns as suggested by the necessary frustration between neighboring spots.

 Fig. 10 Dynamic fingerprints for 5, 7, and 9sided polygons of light spots. Top: ‘on’ and ‘off’ states of the light spots versus time indicated by yellow and blue, respectively. Each row represents one light spot numbered 1 to N from the bottom to the top. Bottom: Radial velocity at each light spot plotted versus time with positive and negative values indicating outward and inward direction, respectively. White and black bars in the bottom right corners represent the characteristic desorption time 1/k^{d}_{A}.  
To quantify the irregularity in the pattern of the nonagon, we calculate the correlation matrix C_{lm} between all pairs of light spots. We define C_{lm} for each respective pair of isomerization functions K_{l} and K_{m}:

 (34) 
with

 (35) 
The resulting paircorrelation matrix is displayed in
Fig. 11. It confirms that, while distant spots are uncorrelated, neighboring spots are anticorrelated. In addition, it reveals some correlations between spots 1 and 2 and 6 and 7 which we expect to vanish in the longtime limit. The correlations between direct neighbors suggest a strong coupling between them, which is present in all studied polygons despite the fact that the projection angle in
eqn (32) approaches 90° for large
N.

 Fig. 11 Colorcoded correlation matrix C_{lm} of isomerization functions for all pairs of light spots l and m for the on–off pattern of the nonagon in Fig. 10. Green indicates inphase correlation, brown indicates antiphase correlation, and dark shading (in both cases) indicates strong correlations.  
Finally, we note for all polygons we observe that our feedback mechanism leads to net advection out of the polygons, which is indicated by the mainly blue regions in Fig. 9 and 10. So, in their centers fluid is drawn from the neighboring fluid phases towards the interface and the polygons of light spots act as pumps.
5 Conclusions
We have developed a theory for photoresponsive surfactants at a fluid–fluid interface and investigated their response to varying patterns of light. In particular, we applied feedback control to generate them. Our theory is formulated in a set of equations for the coarsegrained densities of each surfactant isomer and includes diffusive and advective Marangoni currents at the interface as well as ad/desorption and photoisomerization.
We applied our theory to surfactants with isomers of identical desorption rates. Their densities adjust to a single light spot diffusively, which reverses the direction of the radial Marangoni flow in time. Inspired by this observation, we introduced a feedback coupling between the local advection velocities and isomerization rates of two light spots and thereby generated an alternating switching of the spots. Dynamic fingerprints illustrate their temporal evolution. When we extended the feedback coupling to light spots located on the vertices of regular polygons, we observed different dynamics for even and odd numbers of spots. Evensided polygons, after some possible transient regime, settle into periodically oscillating patterns, where two sets of nextnearest neighbor spots alternate with each other reminiscent of the oscillation in the twospot system. Oddsided polygons with 3 or 5 vertices show similar behavior, however, two neighboring spots determined by initial conditions are always frustrated and switch approximately in phase. Increasing the number of edges to 7 or 9, this frustration produces irregular patterns up to the end of our observation window. The on/off states of nearest neighbor light spots are strongly anticorrelated, while distant spots show vanishing correlations.
Our findings provide a further example how feedback control can initiate complex spatiotemporal patterns in soft matter systems.^{36,37} The observed oscillations are similar to synchronised states in networks of colloidal oscillators which can be in or antiphasic depending on the number of particles involved,^{57} and networks of chemical oscillators which implement logic gates using the Belousov–Zhabotinsky reaction.^{58} While our work is motivated by surfactants which respond directly to light, similar phenomena may be observed with surfactants which respond to light indirectly via lightinduced changes in pH.^{59} We will continue in the direction outlined in this article and investigate, for example, patterns, where the light spots are driven by timedelayed feedback,^{60} which is known to either stabilize or destabilize nonlinear dynamic systems and to initiate complex patterns in space and time.^{61} In this respect, regular polygons are a narrow class of networks, which we aim to extend to more general geometries with given “connectivities” between the light spots.
Our work combines the idea of lightdriven digital microfluidics with the concept of feedback control. Overall, the implemented feedback rule generates flow fields with either regular or irregular oscillations. Thereby, a selfregulated fluid pump occurs that moves fluid out of the polygons. Since both pumping modes are produced by specific light patterns, rather than fixed geometric constraints, one can switch between both modes almost instantaneously and without mechanical adjustments. We have restricted ourselves here to pointlike spots of light. However, different light patterns produced by moving light beams, holographic technique, and light interference provide further fascinating avenues of investigation. Thus, varying light patterns in combination with photoresponsive surfactants offers the possibility of very flexible microfluidic applications.^{21}
Conflicts of interest
There are no conflicts of interest to declare.
A Surface tension
Surface tension σ is the thermodynamic conjugate force to surface area. Therefore, 
 (36) 
where F is, e.g., free energy. It is a homogeneous function, thus F = Af, where the free energy density f only depends on particle densities ρ_{A} and ρ_{B}. Therefore we can write, 
 (37) 
The particle densities are the ratios of particle numbers N_{A} and N_{B} to A,

 (38) 
and we obtain

 (39) 
or after some simplifications:

 (40) 
B Connection to Flory–Huggins theory
A modern derivation of the Flory–Huggins free energy of a binary mixture can be found in ref. 62 (Chapter 4). We start with an analogous expression for ternary mixtures, where the third component is free space with local proportion ϕ_{free}, 
ρ_{∞}^{−1}f = k_{B}T(ϕ_{A}lnϕ_{A} + ϕ_{B}lnϕ_{B} + ϕ_{free}lnϕ_{free}) + χ_{AB}ϕ_{A}ϕ_{B} + χ_{A}ϕ_{A}ϕ_{free} + χ_{B}ϕ_{B}ϕ_{free}.  (41) 
The entropic contribution from ϕ_{free} = 1 − ϕ_{A} − ϕ_{B} effectively describes a hardcore interaction between components A and B: its derivatives with respect to ϕ_{A/B} (the chemical potentials μ_{A/B}) diverge as ϕ_{A} + ϕ_{B} approach 1. Like the other interaction terms we approximate it to second order in ϕ_{A} and ϕ_{B} using (1 − ϕ)ln(1 − ϕ) ≈ −ϕ + 0.5ϕ^{2} with ϕ = ϕ_{A} + ϕ_{B}. Regrouping terms, we find 
 (42) 
In this form, f is equivalent to eqn (2), because 0th and 1st order terms do not contribute to the gradients of surface tension or chemical potentials. We can identify related quantities from both equations:

 (43) 
C Oseen tensor integration by parts
Consider Gauss's theorem applied to the following boundary integral on a compact region ⊂ ^{2} with boundary ∂ and boundary normal vector n(x): 
 (44) 
We apply the product rule to the righthand side: 
 (45) 
To show the equivalence of the second and third term, which correspond, respectively, to v from eqn (11) and (12), if → ^{2}, we need to prove that the boundary integral in eqn (45) vanishes. Consider the integral at a finite, but large distance R from the center, where σ(x) = σ^{eq}:

 (46) 
Substituting O results in the integrand:

 (47) 
Most of the resulting integrals vanish as
R → ∞ because their integrands have order 
x′
^{−2} or 
x′
^{−3}. Only one integral is nontrivial because its integrand's order is 
x′
^{−1}. We prove its convergence to zero by substituting
x′ = (
Rcos
θ,
Rsin
θ) explicitely:

 (48) 
For any
r = (
x,
y,
z) the denominator of the integrand goes to 1 as
R approaches infinity and, therefore,

 (49) 
Thus, the boundary integral in
eqn (45) indeed vanishes.
D Estimating energy constants
In their study of AzoTAB surfactants at an air–water interface Chevallier et al. provide surface tension measurements for two different cis–trans mixing ratios at various bulk concentrations (Fig. 2 in ref. 19). They measured surface tension up to 100 seconds after adding surfactants into their apparatus. For large bulk concentrations their data show a clear saturation of surface tension in time. We rely on these measured values of surface tension, which we show in the last column of Table 2 together with the bulk concentrations of both surfactants (first and second column). In addition, we use σ_{0} = 72.8 mN m^{−1} as the surface tension of a clean air–water interface at room temperature.
Table 2 Our mapping from bulk concentrations C^{eq}_{i} to surface densities ρ^{eq}_{i} and surface tension measurements σ_{eq} of data published by Chevallier et al. in ref. 19
C
^{eq}_{
trans
} (mol L^{−1}) 
C
^{eq}_{
cis
} (mol L^{−1}) 
ρ
^{eq}_{
trans
}

ρ
^{eq}_{
cis
}

σ
^{eq} (mN m^{−1}) 
11.4 
5.9 
0.96ρ_{∞} 
0.02ρ_{∞} 
37 
1.1 
5.8 
0.80ρ_{∞} 
0.14ρ_{∞} 
41 
Chevallier et al. also provide a set of equations to map bulk concentrations C^{eq}_{i} to equilibrium values of surface densities ρ^{eq}_{i},

 (50) 
The surface densities calculated from these equations are displayed in Table 2 (third and fourth column). Finally, we determine the energy constants E_{A} and E_{B} – assuming E_{AB} = (E_{A} + E_{B})/2 – by inserting both sets of data points into eqn (7) for the surface tension.
E Rescaling length and time
In Section 2.6 we reduce the dimensions of all physical quantities, using a set of characteristic parameters. Some choices, like the use of ρ_{∞} for density or k_{B}T for energy, follow naturally from the definition of f in eqn (2). However, indentifying a characteristic length and time is not as straightforward.
Therefore, we consider an idealised system, where we assume density ρ_{B} to be sufficiently large so that during photoisomerization its value remains approximately constant. Assuming also vanishing interactions between the surfactants (E_{A} = E_{B} = E_{AB} = 0), the free energy density f is

f = σ_{0} + k_{B}Tρ_{A}[ln(λ^{2}ρ_{A}) − 1],  (51) 
where the constant contribution from the B surfactant is subsumed in
σ_{0}. Using
eqn (6) for surface tension
σ and chemical potential
μ_{A} = ∂
_{ρA}f, we find

μ_{A} = k_{B}Tln(λ^{2}ρ_{A}) and σ = σ_{0} − k_{B}Tρ_{A}.  (52) 
We further introduce isomerization rates

k_{B→A} = K_{iso}ρ_{B}δ(x) and k_{A→B} = 0.  (53) 
Following the same procedure as in the main text, we arrive at a diffusion–advection–reaction equation for
ρ_{A},

 (54) 
We linearize the equation by substituting
ρ_{A}(
x) =
ρ^{eq}_{A} +
δρ_{A}(
x) and by discarding all terms of second and higher order in
δρ_{A}:

 (55) 
Taking the Fourier transform in both spatial dimensions, turns the partial differential equation into an ordinary differential equation:

 (56) 
The steadystate solution is obtained by setting
∂_{t}δ_{A} = 0. Note that
δ_{A} is a naturally dimensionless quantity and we write:

 (57) 
We identify an amplitude
C =
K_{iso}ρ_{B}/
k^{d}_{A} and two lengths, the diffusion length
l_{D} = (
D/
k^{d}_{A})
^{1/2} and the advection length
l_{A} =
k_{B}Tρ^{eq}_{A}(4
k^{d}_{A})
^{−1}. In Section 2.6 we combine
l_{D} and
k^{d}_{A} as characteristic quantities to eliminate
D in the reduced equations. The advection length
l_{A} is more situational because it depends on
ρ^{eq}_{A}. Also note that
K_{iso} only affects the amplitude of
δ_{A} (and
δρ_{A}) but none of its length scales.
F Regridding scheme
As noted in Section 3, FFTs are an efficient way to compute convolutions but they require orthogonal grids to operate on. Therefore, we need to regrid data from the hexagonal grid onto an orthogonal grid to take advantage of FFT performance.
Fig. 12 shows two different orthogonal grids adjusted to the hexagonal lattice, where only in the white rows we can directly transfer the data from the center of the hexagons to the squares. Thus we decided to use both grids for regridding our data. To preserve as much accuracy as possible we interpolate the “missing” rows in each of those orthogonal grids using a higherorder polynomial interpolation method. Note that we expect sharply peaked surface tension profiles because of the light spots and that, generally, there is a danger of oscillatory artifacts in interpolating peaked data. We can avoid these artifacts by proceeding with both grids and using the extra data to increase robustness. With the interpolated rows we have constructed two complete orthogonal grids, on which we perform the FFTbased convolution. After the convolution we recompile a hexagonal grid from both results taking only rows, which were not interpolated in the first regridding. Because the response function in eqn (14) is also peaked, this step corrects the previous interpolation and prevents oscillatory artifacts.

 Fig. 12 Regridding procedure, from left to right: (i) we divide input data into odd (grid 1) and even rows (grid 2). The remaining rows (dark gray) in each grid are interpolated from input data. (ii) Convolutions are performed on both grids to produce results 1 and 2. (iii) Finally, we compile odd rows from result 1 with even rows from result 2 into output data.  
F.1 Validation and relative error
We validate our combined regridding and convolution scheme against discrete convolution using a pointlike input function f(x) = δ(x), discretized as f_{i} analogous to eqn (25). The convolution in eqn (15) can be discretized and summed exactly on a grid with N = 140 cells of area A: 
 (58) 
Compared to this reference result, we calculate the relative error of our FFTbased method in each cell. The relative error displayed in Fig. 13 does not exceed 1.5%, independent of N. (For comparison, a simplified regridding scheme using only one orthogonal grid leads to relative errors exceeding 60% in the same benchmark.)

 Fig. 13 Relative error of our regridding procedure compared to direct convolution. The location of the input peak is indicated by a black dot and the magnitude of the error is indicated by darker shades while over and underestimation are indicated by blue and red colors, respectively.  
Acknowledgements
We thank J. Blaschke, F. Brückelmann, T. Lubensky, and O. Velev for interesting and helpful discussions. JG acknowledges funding from the German Research Foundation (DFG) via International Research Training Group 1524 and the hospitality of the University of Pennsylvania, Philadelphia, where part of the research was conducted. Support from the Collaborative Research Center 910 is also acknowledged.
References
 E. Samiei, M. Tabrizian and M. Hoorfar, Lab Chip, 2016, 16, 2376–2396 RSC.
 J. Zhang, S. Yan, D. Yuan, G. Alici, N.T. Nguyen, M. Ebrahimi Warkiani and W. Li, Lab Chip, 2016, 16, 10–34 RSC.
 L. R. Huang, E. C. Cox, R. H. Austin and J. C. Sturm, Science, 2004, 304, 987–990 CrossRef CAS PubMed.
 C. Prohm, F. Tröltzsch and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2013, 36, 118 CrossRef PubMed.
 O. Otto, P. Rosendahl, A. Mietke, S. Golfier, C. Herold, D. Klaue, S. Girardo, S. Pagliara, A. Ekpenyong, A. Jacobi, M. Wobus, N. Topfner, U. F. Keyser, J. Mansfeld, E. FischerFriedrich and J. Guck, Nat. Methods, 2015, 12, 199–202 CrossRef CAS PubMed.
 C. Schaaf and H. Stark, Soft Matter, 2017, 13, 3544–3555 RSC.
 C. V. Sternling and L. E. Scriven, AIChE J., 1959, 5, 514–523 CrossRef CAS.
 M. Mokbel, K. Schwarzenberger, K. Eckert and S. Aland, Int. J. Heat Mass Transfer, 2017, 115, 1064–1073 CrossRef CAS.
 N. Garnier, R. O. Grigoriev and M. F. Schatz, Phys. Rev. Lett., 2003, 91, 054501 CrossRef PubMed.
 M. Gugliotti, M. S. Baptista, M. J. Politi, T. P. Silverstein and C. D. Slater, J. Chem. Educ., 2004, 81, 824 CrossRef CAS.
 J. Thomson, Philos. Mag., 1855, 10, 330–333 Search PubMed.
 J. Bławzdziewicz, V. Cristini and M. Loewenberg, Phys. Fluids, 1999, 11, 251–258 CrossRef.
 M. M. Hanczyc, T. Toyota, T. Ikegami, N. Packard and T. Sugawara, J. Am. Chem. Soc., 2007, 129, 9386–9391 CrossRef CAS PubMed.
 S. Thutupalli, R. Seemann and S. Herminghaus, New J. Phys., 2011, 13, 073021 CrossRef.
 S. Shinkai, K. Matsuo, A. Harada and O. Manabe, J. Chem. Soc., Perkin Trans. 2, 1982, 1261–1265 RSC.
 H. Sakai, H. Ebana, K. Sakai, K. Tsuchiya, T. Ohkubo and M. Abe, J. Colloid Interface Sci., 2007, 316, 1027–1030 CrossRef CAS PubMed.
 J. Eastoe, M. S. Dominguez, P. Wyatt, A. Beeby and R. K. Heenan, Langmuir, 2002, 18, 7837–7844 CrossRef CAS.
 T. Shang, K. A. Smith and T. A. Hatton, Langmuir, 2003, 19, 10764–10773 CrossRef CAS.
 E. Chevallier, A. Mamane, H. A. Stone, C. Tribet, F. Lequeux and C. Monteux, Soft Matter, 2011, 7, 7866–7874 RSC.
 M. Schmitt and H. Stark, Phys. Fluids, 2016, 28, 012106 Search PubMed.
 D. Baigl, Lab Chip, 2012, 12, 3637–3653 RSC.
 N. Kavokine, M. Anyfantakis, M. Morel, S. Rudiuk, T. Bickel and D. Baigl, Angew. Chem., Int. Ed., 2016, 55, 11301 CrossRef CAS.
 E. Chevallier, C. Monteux, F. Lequeux and C. Tribet, Langmuir, 2012, 28, 2308–2312 CrossRef CAS PubMed.
 A.L. Fameau, A. Carl, A. SaintJalmes and R. von Klitzing, Phys. Chem. Chem. Phys., 2015, 16, 66–75 CrossRef CAS PubMed.
 C. T. Lee, K. A. Smith and T. A. Hatton, Macromolecules, 2004, 37, 5397–5405 CrossRef CAS.
 A. Kopyshev, N. Lomadze, D. Feldmann, J. Genzer and S. Santer, Polymer, 2015, 79, 65–72 CrossRef CAS.
 A. Diguet, R.M. Guillermic, N. Magome, A. SaintJalmes, Y. Chen, K. Yoshikawa and D. Baigl, Angew. Chem., Int. Ed., 2009, 48, 9281–9284 CrossRef CAS PubMed.
 D. Zinemanas and A. Nir, J. Fluid Mech., 1988, 193, 217 CrossRef.
 A. Nadim, Chem. Eng. Commun., 1996, 148150, 391–407 CrossRef CAS.
 C. Pozrikidis, J. Comput. Phys., 2001, 169, 250–301 CrossRef CAS.
 M. Schmitt and H. Stark, EPL, 2013, 101, 44008 CrossRef.
 G. Möller, M. Harke, H. Motschmann and D. Prescher, Langmuir, 1998, 14, 4955–4957 CrossRef.

K. J. Åström and R. M. Murray, Feedback Systems: An Introduction for Scientists and Engineers, Princeton University Press, 2008 Search PubMed.
 A. Friedman and L.S. Jiang, Commun. Partial Differ. Equ., 1988, 13, 515–550 CrossRef.
 P. Gurevich, W. Jäger and A. Skubachevskii, SIAM J. Math. Anal., 2009, 41, 733–752 CrossRef.
 C. Prohm and H. Stark, Lab Chip, 2014, 14, 2115–2123 RSC.
 M. Zeitz, P. Gurevich and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2015, 38, 22 CrossRef PubMed.
 M. L. Huggins, J. Chem. Phys., 1941, 9, 440 CrossRef CAS.
 P. J. Flory, J. Chem. Phys., 1941, 9, 660 CrossRef CAS.

J.P. Hansen and I. R. McDonald, Theory of Simple Liquids, Elsevier, 3rd edn, 2006 Search PubMed.
 U. M. B. Marconi and P. Tarazona, J. Chem. Phys., 1999, 110, 8032–8044 CrossRef CAS.
 A. J. Archer, J. Phys.: Condens. Matter, 2005, 17, 1405 CrossRef CAS.
 R. B. Jones, B. U. Felderhof and J. M. Deutch, Macromolecules, 1975, 8, 680–684 CrossRef CAS.
 A. Dominguez, P. Malgaretti, M. N. Popescu and S. Dietrich, Soft Matter, 2016, 12, 8398–8406 RSC.
 M. Riesz, Acta Math., 1949, 81, 1–222 Search PubMed.
 A. Frumkin, Z. Phys., 1926, 35, 792–802 CrossRef CAS.
 I. Langmuir, J. Am. Chem. Soc., 1918, 40, 1361–1403 CrossRef CAS.
 J. K. Ferri and K. J. Stebe, J. Colloid Interface Sci., 1999, 209, 1–9 CrossRef CAS PubMed.
 J. Bezanson, A. Edelman, S. Karpinski and V. B. Shah, SIAM Rev., 2017, 59, 65–98 CrossRef.
 M. Schmitt and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2016, 39, 80 CrossRef CAS PubMed.

T. Pang, An Introduction to Computational Physics, Cambridge University Press, 2010 Search PubMed.

Numerical Methods for AdvectionDiffusion Problems, ed. C. B. Vreugdenhil and B. Koren, Vieweg, 1993, pp. 117–138 Search PubMed.
 A. Harten, J. Comput. Phys., 1983, 49, 357–393 CrossRef.
 P. K. Sweby, SIAM J. Numer. Anal., 1984, 21, 995–1011 CrossRef.
 E. Fehlberg, Computing, 1970, 6, 61–71 CrossRef.

Digital Library of Mathematical Functions, ed. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller and B. V. Saunders, 2017, NIST, Release 1.0.16 Search PubMed.
 L. Damet, G. M. Cicuta, J. Kotar, M. C. Lagomarsino and P. Cicuta, Soft Matter, 2012, 8, 8672–8678 RSC.
 A. Wang, J. Gold, N. Tompkins, M. Heymann, K. Harrington and S. Fraden, Eur. Phys. J.Spec. Top., 2016, 225, 211–227 CrossRef CAS PubMed.
 P. Cicuta, I. Hopkinson and P. G. Petrov, J. Chem. Phys., 2001, 115, 9991–9994 CrossRef CAS.
 K. Pyragas, Phys. Lett. A, 1992, 170, 421–428 CrossRef.

Handbook of Chaos Control, ed. E. Schöll and H. G. Schuster, WileyVCH, 2nd edn, 2008 Search PubMed.

P. C. Hiemenz and T. P. Lodge, Polymer Chemistry, CRC Press, 2nd edn, 2007 Search PubMed.
Footnotes 
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7sm02101a 
‡ Note that σ_{eq} refers to an interface in equilibrium while σ_{0} refers to a clean interface; their values are usually different. 
§ Note that this is not valid for AzoTAB at an air–water interface as demonstrated by Chevallier et al.^{19} 
¶ More precisely, Koren's method is total variation diminishing in Harten's sense^{53} and robust in Sweby's sense.^{54} 
 This is similar to the steady state of a simple diffusion equation with a pointsource and constant global decay rate in two dimensions. The analytical solution for that system is the 0th order modified Bessel function of the second kind _{0}(r). For large r it scales as r^{−1/2}e^{−r} (eqn (10.40.2) in ref. 56). 

This journal is © The Royal Society of Chemistry 2018 