Feedback control of photoresponsive fluid interfaces

Abstract

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 reactors¹ and precise sorting for micron-sized objects.² 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,¹⁵ spiropyran,¹⁶ or stilbene¹⁷ are photoresponsive: they reversibly change shape under illumination with specific wavelengths of light and thereby interact differently with bordering fluids and neighboring surfactants.¹⁸ 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,²⁵ and photoresponsive polymer brushes.²⁶ Beyond tunable materials, they allow to move and position individual surfactant-covered droplets,²⁷ to turn droplets into self-propelled particles,²⁰ and to pump fluid near a planar water surface.¹⁹ 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 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.

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,³³ 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.²¹

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 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,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 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 . 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 σ₀ of a clean interface without surfactants:Entropic terms are proportional to temperature k_BT 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⁴¹ for colloids and generalized to binary mixtures by Archer.⁴² 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 j_q, which depend on functional derivatives of F:

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_BT/γ, we arrive atLinear 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

Like f, σ is based on the usual adiabatic approximation of local equilibrium. We substitute f from eqn (2) and arrive atIntuitively, 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)∇_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:

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

The tensor depends only on the average viscosity [small eta, Greek, macron] = (η₊ + η₋)/2. The velocity field v initiated by an arbitrary surface tension profile σ isWe 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 fieldwhere the velocity response function g(r) is a three-dimensional vector field, which decays like |r|⁻²:

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

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:

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,⁴⁵ a form of the inverse half Laplacian (−Δ_x)^−1/2 operating on δσ. The operator is defined by its Fourier transform, such that

where [scr F, script letter F] denotes the two-dimensional Fourier transform, and [scr F, script letter F]⁻¹ 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.⁴⁶ Because it takes into account surface coverage and local interactions on the interface, it is an extension of the Langmuir adsorption model.⁴⁷ 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) readwhere 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σ.⁴⁸ 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 solvingfor 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 point-like light spots at positions x_l is described by a superposition of Dirac-δ functions with isomerization functions K_l(t):

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.5kBT
E AB	1.3kBT
E B	1.1kBT
ρ ^eqA	0.5ρ∞
ρ ^eqB	0.0ρ∞
Mg	1.2 × 10^5
k ^dB	1k^dA
k ^aA	2.12k^dA
k ^aB	1.73k^dA
K A→B	10D

---

2.6 Closed and reduced equations

We nondimensionalize all quantities using k_BT 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²/D, to the time for advection by Marangoni flow, t_Mg = L/v_Mg = LΔσ/[small eta, Greek, macron]. Thus,

To arrive at the last expression, we set Δσ = k_BTρ_∞ and used the definition for L.

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

withand 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 ρ^eq_q at the boundary but also outside [scr B, script letter B], so that

δσ(x ∉ [scr B, script letter B]) = σ(ρ^eq_A,ρ^eq_B) − σ_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 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¹⁹ assuming E_AB = (E_A + E_B)/2 (for details see Appendix D). Our choice of Mg = 1.2 × 10⁵ is also based on their data, assuming D = 900 μm² s⁻¹, ρ_∞ = 5 nm⁻², k^d_A = 6.62 s⁻¹, and T = 300 K, and a larger value of the average viscosity [small eta, Greek, macron] = 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.⁴⁹ 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 000 to 127 000 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 [capital script C]_i with area A. For example, the photoisomerization field k(x) = Kδ(x − x̃) for a single light spot is discretized by

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 ϕ:

where r is the upwind gradient ratio between cells i and i + 1 of the advected field ρ_q:

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|⁻¹. Rather than performing the direct convolution, which would need an [scr O, script letter O](N²)-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 log 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.⁵⁵ 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 L¹-norm because for any physically valid ρ_q it simply gives the particle number N_q, which is dimensionless and necessarily exists. The L¹ 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:

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⁻⁶.

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 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 non-equilibrium steady state (dash-dotted, 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_BT(ρ − ρ^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⁵ 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 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 C₄AzoOC₆E₂ have almost identical adsorption behavior.¹⁸ 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₁ and v₂ 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 v₁ and v₂ (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₁ and x₂ with distance Δx = |x₁ − x₂| [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:

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₂ 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₁. Over time the coupled light spots quickly establish regular oscillations in v₁ and v₂, which occur in antiphase. In Fig. 6(c) we represent the on/off state of the light spots and the velocities v₁, v₂ 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 Color-coded 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) − v̄(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.1k^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 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.

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 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.

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/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:

withThe 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.

Fig. 11 Color-coded 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 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,⁵⁷ and networks of chemical oscillators which implement logic gates using the Belousov–Zhabotinsky reaction.⁵⁸ 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.⁵⁹ 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,⁶⁰ which is known to either stabilize or destabilize nonlinear dynamic systems and to initiate complex patterns in space and time.⁶¹ 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.²¹

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,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,

The particle densities are the ratios of particle numbers N_A and N_B to A,

and we obtainor after some simplifications:

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,

ρ_∞⁻¹f = k_BT(ϕ_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 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ϕ² with ϕ = ϕ_A + ϕ_B. Regrouping terms, we find

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:

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]² with boundary ∂[scr A, script letter A] and boundary normal vector n(x):We apply the product rule to the right-hand side:

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]², 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:

Substituting O results in the integrand:

Most of the resulting integrals vanish as R → ∞ because their integrands have order |x′|⁻² or |x′|⁻³. Only one integral is nontrivial because its integrand's order is |x′|⁻¹. We prove its convergence to zero by substituting x′ = (R cos θ,R sin θ) explicitely:For any r = (x,y,z) the denominator of the integrand goes to 1 as R approaches infinity and, therefore,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 σ₀ = 72.8 mN m⁻¹ 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,

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_BT 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 = σ₀ + k_BTρ_A[ln(λ²ρ_A) − 1], (51)

where the constant contribution from the B surfactant is subsumed in σ₀. Using eqn (6) for surface tension σ and chemical potential μ_A = ∂_ρAf, we find

μ_A = k_BT ln(λ²ρ_A) and σ = σ₀ − k_BTρ_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,We linearize the equation by substituting ρ_A(x) = ρ^eq_A + δρ_A(x) and by discarding all terms of second and higher order in δρ_A:Taking the Fourier transform in both spatial dimensions, turns the partial differential equation into an ordinary differential equation: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: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_BTρ^eq_A(4k^d_A[small eta, Greek, macron])⁻¹. 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 δ[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.

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 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:

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.)

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

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

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7sm02101a

‡ Note that σ_eq refers to an interface in equilibrium while σ₀ 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.¹⁹

¶ More precisely, Koren's method is total variation diminishing in Harten's sense⁵³ and robust in Sweby's sense.⁵⁴

|| 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]₀(r). For large r it scales as r^−1/2e^−r (eqn (10.40.2) in ref. 56).

---