Open Access Article
This Open Access Article is licensed under a Creative Commons Attribution-Non Commercial 3.0 Unported Licence

Feedback control of photoresponsive fluid interfaces

Josua Grawitter *ab and Holger Stark a
aInstitut für Theoretische Physik, Technische Universität Berlin, Hardenbergstr. 36, 10623 Berlin, Germany. E-mail:;
bDepartment of Physics and Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104-6396, USA

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 even-sided polygons, where two sets of next-nearest 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 reactors1 and precise sorting for micron-sized 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 curvature7,8 or temperature9,10 gradients, or from gradients in the concentration11–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 stilbene17 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 surfactant-covered droplets,27 to turn droplets into self-propelled 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 descriptions28–31 have not considered photoisomerization, or in the case of ref. 20 considered it only in the high-intensity 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.

image file: c7sm02101a-f1.tif
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 point-like 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 self-generated oscillations of two neighboring spots. We generalize our approach to regular polygons. While we mostly observe regular oscillatory patterns, seven and nine-sided 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 patterns34,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 equations28–31 governing the coupled dynamics of surfactant densities ρq,
image file: c7sm02101a-t1.tif(1)
Here, the diffusive currents jDq result from thermal motion and local interactions of the surfactants, the advection velocity vA from drift due to the Marangoni effect, and the reaction terms Sadsq describe ad- and desorption to and from a fluid–fluid interface, respectively. In addition, we introduce reaction terms Sisoq 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 coarse-grained Helmholtz free energy F of a compressible mixture of two components, A and B, with mean-field interactions. We define F as the integral over a free energy density f on the interface, such that image file: c7sm02101a-t2.tif. The first term of f contains entropic contributions of each component (see e.g. Chapter 3.1 in ref. 40), the second term describes mean-field 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:
image file: c7sm02101a-t3.tif(2)
Entropic terms are proportional to temperature kBT and contain the thermal de Broglie wavelength λ, which does not affect dynamics in our case. Interaction terms are characterized by three energy constants EA, EAB and EB 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 Tarazona41 for colloids and generalized to binary mixtures by Archer.42 The original derivation builds on a variant of the N-particle 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 jq, which depend on functional derivatives of F:

image file: c7sm02101a-t4.tif(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 = kBT/γ, we arrive at
image file: c7sm02101a-t5.tif(4)
image file: c7sm02101a-t6.tif(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 image file: c7sm02101a-t7.tif

image file: c7sm02101a-t8.tif(6)
Like f, σ is based on the usual adiabatic approximation of local equilibrium. We substitute f from eqn (2) and arrive at
image file: c7sm02101a-t9.tif(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 (Ei > 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 = (vx,vy,vz):

r[η(z)∇rv] − ∇rp = −δ(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:
vz(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

image file: c7sm02101a-t10.tif(10)
The tensor depends only on the average viscosity [small eta, Greek, macron] = (η+ + η)/2. The velocity field v initiated by an arbitrary surface tension profile σ is
image file: c7sm02101a-t11.tif(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
image file: c7sm02101a-t12.tif(12)
where the velocity response function g(r) is a three-dimensional vector field, which decays like |r|−2:
image file: c7sm02101a-t13.tif(13)

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

image file: c7sm02101a-f2.tif
Fig. 2 Fluid velocity field initiated by a positive point-like 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:

image file: c7sm02101a-t14.tif(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

image file: c7sm02101a-t15.tif(15)
where [scr F, script letter F] denotes the two-dimensional Fourier transform, and [scr F, script letter F]−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
image file: c7sm02101a-t16.tif(16)
where we also introduced adsorption and desorption rates kaq and kdq 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 kaq 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 ρeqq by solving
image file: c7sm02101a-t17.tif(17)
for given rate ratios kaA/kdA and kaB/kdB.

2.5 Photoisomerization

We describe photoisomerization using rate equations with local isomerization rates kij in each direction:
SisoA(ρA,ρB) = kB→AρBkA→BρA(18)
SisoB(ρA,ρB) = kA→BρAkB→AρB(19)
For our specific examples in Section 4 we only allow for photoinduced isomerization from A to B surfactants and therefore set kB→A = 0. Note that total surfactant density is conserved during photoisomerization, which is quantified by SisoA + SisoB = 0.

The isomerization rates in our examples are spatially inhomogeneous. For example, a set of N point-like light spots at positions xl is described by a superposition of Dirac-δ functions with isomerization functions Kl(t):

image file: c7sm02101a-t18.tif(20)
Each Kl is proportional to the intensity of its light spot l. Thus, when the light spot is turned on, Kl switches from 0 to a constant value KA→B specified in Table 1.

Table 1 Our estimates for various material and photoisomerization parameters. The values for kaA/B are calculated using eqn (17)
Parameter Value
E A 1.5kBT
E AB 1.3kBT
E B 1.1kBT
ρ eqA 0.5ρ
ρ eqB 0.0ρ
Mg 1.2 × 105
k dB 1kdA
k aA 2.12kdA
k aB 1.73kdA
K A→B 10D

2.6 Closed and reduced equations

We nondimensionalize all quantities using kBT for energies, desorption rate kdA for inverse time, saturation density ρ for densities, and image file: c7sm02101a-t19.tif 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, tdiff = L2/D, to the time for advection by Marangoni flow, tMg = L/vMg = LΔσ/[small eta, Greek, macron]. Thus,

image file: c7sm02101a-t20.tif(21)

To arrive at the last expression, we set Δσ = kB and used the definition for L.

Non-dimensionalizing eqn (1), (3), (14), (16) and (19) leads to reduced equations for each component q:

image file: c7sm02101a-t21.tif(22)
image file: c7sm02101a-t22.tif(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 [scr B, script letter B] of the interface, where we solve eqn (22) with open, Dirichlet-like boundary conditions: we set the densities ρq to their equilibrium values ρeqq at the boundary but also outside [scr B, script letter B], so that

δσ(x[scr B, script letter B]) = σ(ρeqA,ρeqB) − σeq = 0.(24)
Thereby, the region of integration of the integral in eqn (14) is restricted to [scr B, script letter B].

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 EA, EB, and EAB are based on data for surface tensions of surfactant mixtures collected by Chevallier et al. for AzoTAB at a water–air interface19 assuming EAB = (EA + EB)/2 (for details see Appendix D). Our choice of Mg = 1.2 × 105 is also based on their data, assuming D = 900 μm2 s−1, ρ = 5 nm−2, kdA = 6.62 s−1, and T = 300 K, and a larger value of the average viscosity [small eta, Greek, macron] = 2 mPa s.

The equilibrium surface densities ρeqA and ρeqB are assumed to be fully adjustable by changing the partial pressure of surfactants in the fluids, and for simplicity, we set ρeqB = 0 and ρeqA = 0.5ρ. For the same reason we assume identical desorption rates kdB = kdA 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 KA→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 56[thin space (1/6-em)]000 to 127[thin space (1/6-em)]000 hexagons (depending on the problem) and approximate all functions f by their values fi at the center xi of each cell i (similar to ref. 50). If we know f exactly, we approximate it by the spatial average over each cell [capital script C]i with area A. For example, the photoisomerization field k(x) = (x[x with combining tilde]) for a single light spot is discretized by

image file: c7sm02101a-t23.tif(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 upwind-difference 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 flux-limiting function ϕ:

image file: c7sm02101a-t24.tif(26)
where r is the upwind gradient ratio between cells i and i + 1 of the advected field ρq:
image file: c7sm02101a-t25.tif(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 [scr O, script letter O](N2)-operation on a mesh with N cells, we compute the integral in Fourier space. The fast Fourier transform (FFT) is an [scr O, script letter O](N[thin space (1/6-em)]log[thin space (1/6-em)]N)-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 5th-order 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 5th-order 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 L1-norm because for any physically valid ρq it simply gives the particle number Nq, which is dimensionless and necessarily exists. The L1 norm also has a straight-forward 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:

image file: c7sm02101a-t26.tif(28)

The integration error e is simply the norm of the point-wise difference between two states ρ and [small rho, Greek, tilde]:

e[ρ,[small rho, Greek, tilde]] = ‖ρ(x) − [small rho, Greek, tilde](x)‖,(29)
where ρ(x), [small rho, Greek, tilde](x) refer to the 5th- and 4th-order 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 point-like 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 vr, and surface tension δσ relative to their equilibrium values are displayed in Fig. 3 for various times t.

image file: c7sm02101a-f3.tif
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 vr, 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 vr 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.5kdA.

image file: c7sm02101a-f4.tif
Fig. 4 Radial velocity vr at distance r = 5L from the light spot plotted versus time t for light switched on (solid, blue) and switched off in the non-equilibrium steady state (dash-dotted, green). Emblem: a sensor measures vr 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 vr 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 (vr < 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 (vr > 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 δσ ≈ −kBT(ρρ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 vr 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 vr relaxes towards zero. We now introduce a feedback control step in our system: when vr changes sign, we switch off the light spot. The resulting time evolution of vr 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 (kdB = 300kdA) 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.

image file: c7sm02101a-f5.tif
Fig. 5 Radial profiles at various times t after a light spot at r = 0 is switched on for Mg = 6 × 105 and kdB = 300kdA: (a) surface tension δσ relative to equilibrium value (inset: surfactant density ρB), and (b) radial velocity vr. The solid line in each diagram shows the steady-state 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 C4AzoOC6E2 have almost identical adsorption behavior.18 In the following, we concentrate on the more interesting case kdB = kdA.

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.
image file: c7sm02101a-f6.tif
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 v1 and v2 sensed at the locations of the respective light spots 1 and 2 plotted over time t. (c) Dynamic fingerprints of two light spots: color-coded time series for light spots (on: yellow/light, off: blue/dark) and velocities v1 and v2 (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 x1 and x2 with distance Δx = |x1x2| [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 vl at one spot l to be positive, when it points away from the other spot and to be negative in the opposite case:

image file: c7sm02101a-t27.tif(30)

Only the spot with the larger vl is switched on or equivalently, the isomerization function Kl(t) of each spot, which we introduced in eqn (20), is

Kl(t) = KA→Bθ(vl(t) −vj(t)), jl.(31)
Here, θ is the Heaviside step function and KA→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 v2 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 v1. Over time the coupled light spots quickly establish regular oscillations in v1 and v2, which occur in antiphase. In Fig. 6(c) we represent the on/off state of the light spots and the velocities v1, v2 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 KA→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 kdA, 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.

image file: c7sm02101a-f7.tif
Fig. 7 Color-coded switching frequency f of two coupled light spots plotted versus distance Δx and isomerization constant KA→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 vl for each spot l:
vl = v(xlel.(32)

image file: c7sm02101a-f8.tif
Fig. 8 Schematic of 6 light spots positioned on the vertices of a regular hexagon, each with a radial projection vector el pointing away from the polygon's center.

Each individual spot l is only kept switched on while vl is larger than the average of all vl. Equivalently, we introduce the isomerization function for spot l:

Kl(t) = KA→Bθ(vl(t) − [v with combining macron](t))(33)
with mean velocity image file: c7sm02101a-t28.tif 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.1kdA. 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 non-neighboring 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.

image file: c7sm02101a-f9.tif
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/kdA.

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 long-lived 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 neighbor-neighbor 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 odd-numbered polygons settle into irregular patterns as suggested by the necessary frustration between neighboring spots.

image file: c7sm02101a-f10.tif
Fig. 10 Dynamic fingerprints for 5, 7, and 9-sided 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/kdA.

To quantify the irregularity in the pattern of the nonagon, we calculate the correlation matrix Clm between all pairs of light spots. We define Clm for each respective pair of isomerization functions Kl and Km:

image file: c7sm02101a-t29.tif(34)
image file: c7sm02101a-t30.tif(35)
The resulting pair-correlation 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 long-time 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.

image file: c7sm02101a-f11.tif
Fig. 11 Color-coded correlation matrix Clm 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 in-phase correlation, brown indicates anti-phase 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 coarse-grained 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. Even-sided polygons, after some possible transient regime, settle into periodically oscillating patterns, where two sets of next-nearest neighbor spots alternate with each other reminiscent of the oscillation in the two-spot system. Odd-sided 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 anti-correlated, while distant spots show vanishing correlations.

Our findings provide a further example how feedback control can initiate complex spatio-temporal 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 anti-phasic 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 light-induced 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 time-delayed 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 light-driven digital microfluidics with the concept of feedback control. Overall, the implemented feedback rule generates flow fields with either regular or irregular oscillations. Thereby, a self-regulated 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 point-like 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,
image file: c7sm02101a-t31.tif(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,
image file: c7sm02101a-t32.tif(37)

The particle densities are the ratios of particle numbers NA and NB to A,

image file: c7sm02101a-t33.tif(38)
and we obtain
image file: c7sm02101a-t34.tif(39)
or after some simplifications:
image file: c7sm02101a-t35.tif(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,
ρ−1f = kBT(ϕA[thin space (1/6-em)]ln[thin space (1/6-em)]ϕA + ϕB[thin space (1/6-em)]ln[thin space (1/6-em)]ϕB + ϕfree[thin space (1/6-em)]ln[thin space (1/6-em)]ϕfree) + χABϕAϕB + χAϕAϕfree + χBϕBϕfree.(41)
The entropic contribution from ϕfree = 1 − ϕAϕB effectively describes a hard-core 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
image file: c7sm02101a-t36.tif(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:

image file: c7sm02101a-t37.tif(43)

C Oseen tensor integration by parts

Consider Gauss's theorem applied to the following boundary integral on a compact region [scr A, script letter A][Doublestruck R]2 with boundary ∂[scr A, script letter A] and boundary normal vector n(x):
image file: c7sm02101a-t38.tif(44)
We apply the product rule to the right-hand side:
image file: c7sm02101a-t39.tif(45)

To show the equivalence of the second and third term, which correspond, respectively, to v from eqn (11) and (12), if [scr A, script letter A][Doublestruck R]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:

image file: c7sm02101a-t40.tif(46)

Substituting O results in the integrand:

image file: c7sm02101a-t41.tif(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′ = (R[thin space (1/6-em)]cos[thin space (1/6-em)]θ,R[thin space (1/6-em)]sin[thin space (1/6-em)]θ) explicitely:
image file: c7sm02101a-t42.tif(48)
For any r = (x,y,z) the denominator of the integrand goes to 1 as R approaches infinity and, therefore,
image file: c7sm02101a-t43.tif(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 cistrans 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 Ceqi to surface densities ρeqi 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 Ceqi to equilibrium values of surface densities ρeqi,

image file: c7sm02101a-t44.tif(50)

The surface densities calculated from these equations are displayed in Table 2 (third and fourth column). Finally, we determine the energy constants EA and EB – assuming EAB = (EA + EB)/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 kBT 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 (EA = EB = EAB = 0), the free energy density f is

f = σ0 + kBA[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 = ∂ρAf, we find
μA = kBT[thin space (1/6-em)]ln(λ2ρA) and σ = σ0kBA.(52)
We further introduce isomerization rates
kB→A = KisoρBδ(x) and kA→B = 0.(53)
Following the same procedure as in the main text, we arrive at a diffusion–advection–reaction equation for ρA,
image file: c7sm02101a-t45.tif(54)
We linearize the equation by substituting ρA(x) = ρeqA + δρA(x) and by discarding all terms of second and higher order in δρA:
image file: c7sm02101a-t46.tif(55)
Taking the Fourier transform in both spatial dimensions, turns the partial differential equation into an ordinary differential equation:
image file: c7sm02101a-t47.tif(56)
The steady-state solution is obtained by setting tδ[small rho, Greek, circumflex]A = 0. Note that δ[small rho, Greek, circumflex]A is a naturally dimensionless quantity and we write:
image file: c7sm02101a-t48.tif(57)
We identify an amplitude C = KisoρB/kdA and two lengths, the diffusion length lD = (D/kdA)1/2 and the advection length lA = kBeqA(4kdA[small eta, Greek, macron])−1. In Section 2.6 we combine lD and kdA as characteristic quantities to eliminate D in the reduced equations. The advection length lA is more situational because it depends on ρeqA. Also note that Kiso only affects the amplitude of δ[small rho, Greek, circumflex]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 higher-order 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 FFT-based 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.

image file: c7sm02101a-f12.tif
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 point-like input function f(x) = δ(x), discretized as fi 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:
image file: c7sm02101a-t49.tif(58)

Compared to this reference result, we calculate the relative error of our FFT-based 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.)

image file: c7sm02101a-f13.tif
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.


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.


  1. E. Samiei, M. Tabrizian and M. Hoorfar, Lab Chip, 2016, 16, 2376–2396 RSC.
  2. J. Zhang, S. Yan, D. Yuan, G. Alici, N.-T. Nguyen, M. Ebrahimi Warkiani and W. Li, Lab Chip, 2016, 16, 10–34 RSC.
  3. L. R. Huang, E. C. Cox, R. H. Austin and J. C. Sturm, Science, 2004, 304, 987–990 CrossRef CAS PubMed.
  4. C. Prohm, F. Tröltzsch and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2013, 36, 118 CrossRef PubMed.
  5. 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. Fischer-Friedrich and J. Guck, Nat. Methods, 2015, 12, 199–202 CrossRef CAS PubMed.
  6. C. Schaaf and H. Stark, Soft Matter, 2017, 13, 3544–3555 RSC.
  7. C. V. Sternling and L. E. Scriven, AIChE J., 1959, 5, 514–523 CrossRef CAS.
  8. M. Mokbel, K. Schwarzenberger, K. Eckert and S. Aland, Int. J. Heat Mass Transfer, 2017, 115, 1064–1073 CrossRef CAS.
  9. N. Garnier, R. O. Grigoriev and M. F. Schatz, Phys. Rev. Lett., 2003, 91, 054501 CrossRef PubMed.
  10. M. Gugliotti, M. S. Baptista, M. J. Politi, T. P. Silverstein and C. D. Slater, J. Chem. Educ., 2004, 81, 824 CrossRef CAS.
  11. J. Thomson, Philos. Mag., 1855, 10, 330–333 Search PubMed.
  12. J. Bławzdziewicz, V. Cristini and M. Loewenberg, Phys. Fluids, 1999, 11, 251–258 CrossRef.
  13. M. M. Hanczyc, T. Toyota, T. Ikegami, N. Packard and T. Sugawara, J. Am. Chem. Soc., 2007, 129, 9386–9391 CrossRef CAS PubMed.
  14. S. Thutupalli, R. Seemann and S. Herminghaus, New J. Phys., 2011, 13, 073021 CrossRef.
  15. S. Shinkai, K. Matsuo, A. Harada and O. Manabe, J. Chem. Soc., Perkin Trans. 2, 1982, 1261–1265 RSC.
  16. H. Sakai, H. Ebana, K. Sakai, K. Tsuchiya, T. Ohkubo and M. Abe, J. Colloid Interface Sci., 2007, 316, 1027–1030 CrossRef CAS PubMed.
  17. J. Eastoe, M. S. Dominguez, P. Wyatt, A. Beeby and R. K. Heenan, Langmuir, 2002, 18, 7837–7844 CrossRef CAS.
  18. T. Shang, K. A. Smith and T. A. Hatton, Langmuir, 2003, 19, 10764–10773 CrossRef CAS.
  19. E. Chevallier, A. Mamane, H. A. Stone, C. Tribet, F. Lequeux and C. Monteux, Soft Matter, 2011, 7, 7866–7874 RSC.
  20. M. Schmitt and H. Stark, Phys. Fluids, 2016, 28, 012106 Search PubMed.
  21. D. Baigl, Lab Chip, 2012, 12, 3637–3653 RSC.
  22. N. Kavokine, M. Anyfantakis, M. Morel, S. Rudiuk, T. Bickel and D. Baigl, Angew. Chem., Int. Ed., 2016, 55, 11301 CrossRef CAS.
  23. E. Chevallier, C. Monteux, F. Lequeux and C. Tribet, Langmuir, 2012, 28, 2308–2312 CrossRef CAS PubMed.
  24. A.-L. Fameau, A. Carl, A. Saint-Jalmes and R. von Klitzing, Phys. Chem. Chem. Phys., 2015, 16, 66–75 CrossRef CAS PubMed.
  25. C. T. Lee, K. A. Smith and T. A. Hatton, Macromolecules, 2004, 37, 5397–5405 CrossRef CAS.
  26. A. Kopyshev, N. Lomadze, D. Feldmann, J. Genzer and S. Santer, Polymer, 2015, 79, 65–72 CrossRef CAS.
  27. A. Diguet, R.-M. Guillermic, N. Magome, A. Saint-Jalmes, Y. Chen, K. Yoshikawa and D. Baigl, Angew. Chem., Int. Ed., 2009, 48, 9281–9284 CrossRef CAS PubMed.
  28. D. Zinemanas and A. Nir, J. Fluid Mech., 1988, 193, 217 CrossRef.
  29. A. Nadim, Chem. Eng. Commun., 1996, 148-150, 391–407 CrossRef CAS.
  30. C. Pozrikidis, J. Comput. Phys., 2001, 169, 250–301 CrossRef CAS.
  31. M. Schmitt and H. Stark, EPL, 2013, 101, 44008 CrossRef.
  32. G. Möller, M. Harke, H. Motschmann and D. Prescher, Langmuir, 1998, 14, 4955–4957 CrossRef.
  33. K. J. Åström and R. M. Murray, Feedback Systems: An Introduction for Scientists and Engineers, Princeton University Press, 2008 Search PubMed.
  34. A. Friedman and L.-S. Jiang, Commun. Partial Differ. Equ., 1988, 13, 515–550 CrossRef.
  35. P. Gurevich, W. Jäger and A. Skubachevskii, SIAM J. Math. Anal., 2009, 41, 733–752 CrossRef.
  36. C. Prohm and H. Stark, Lab Chip, 2014, 14, 2115–2123 RSC.
  37. M. Zeitz, P. Gurevich and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2015, 38, 22 CrossRef PubMed.
  38. M. L. Huggins, J. Chem. Phys., 1941, 9, 440 CrossRef CAS.
  39. P. J. Flory, J. Chem. Phys., 1941, 9, 660 CrossRef CAS.
  40. J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, Elsevier, 3rd edn, 2006 Search PubMed.
  41. U. M. B. Marconi and P. Tarazona, J. Chem. Phys., 1999, 110, 8032–8044 CrossRef CAS.
  42. A. J. Archer, J. Phys.: Condens. Matter, 2005, 17, 1405 CrossRef CAS.
  43. R. B. Jones, B. U. Felderhof and J. M. Deutch, Macromolecules, 1975, 8, 680–684 CrossRef CAS.
  44. A. Dominguez, P. Malgaretti, M. N. Popescu and S. Dietrich, Soft Matter, 2016, 12, 8398–8406 RSC.
  45. M. Riesz, Acta Math., 1949, 81, 1–222 Search PubMed.
  46. A. Frumkin, Z. Phys., 1926, 35, 792–802 CrossRef CAS.
  47. I. Langmuir, J. Am. Chem. Soc., 1918, 40, 1361–1403 CrossRef CAS.
  48. J. K. Ferri and K. J. Stebe, J. Colloid Interface Sci., 1999, 209, 1–9 CrossRef CAS PubMed.
  49. J. Bezanson, A. Edelman, S. Karpinski and V. B. Shah, SIAM Rev., 2017, 59, 65–98 CrossRef.
  50. M. Schmitt and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2016, 39, 80 CrossRef CAS PubMed.
  51. T. Pang, An Introduction to Computational Physics, Cambridge University Press, 2010 Search PubMed.
  52. Numerical Methods for Advection-Diffusion Problems, ed. C. B. Vreugdenhil and B. Koren, Vieweg, 1993, pp. 117–138 Search PubMed.
  53. A. Harten, J. Comput. Phys., 1983, 49, 357–393 CrossRef.
  54. P. K. Sweby, SIAM J. Numer. Anal., 1984, 21, 995–1011 CrossRef.
  55. E. Fehlberg, Computing, 1970, 6, 61–71 CrossRef.
  56. 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.
  57. L. Damet, G. M. Cicuta, J. Kotar, M. C. Lagomarsino and P. Cicuta, Soft Matter, 2012, 8, 8672–8678 RSC.
  58. 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.
  59. P. Cicuta, I. Hopkinson and P. G. Petrov, J. Chem. Phys., 2001, 115, 9991–9994 CrossRef CAS.
  60. K. Pyragas, Phys. Lett. A, 1992, 170, 421–428 CrossRef.
  61. Handbook of Chaos Control, ed. E. Schöll and H. G. Schuster, Wiley-VCH, 2nd edn, 2008 Search PubMed.
  62. P. C. Hiemenz and T. P. Lodge, Polymer Chemistry, CRC Press, 2nd edn, 2007 Search PubMed.


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 sense53 and robust in Sweby's sense.54
|| This is similar to the steady state of a simple diffusion equation with a point-source 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 [scr K, script letter K]0(r). For large r it scales as r−1/2er (eqn (10.40.2) in ref. 56).

This journal is © The Royal Society of Chemistry 2018