Chemomechanical simulation of soap film flow on spherical bubbles

Soap bubbles are widely appreciated for their fragile nature and their colorful appearance. The natural sciences and, in extension, computer graphics, have comprehensively studied the mechanical behavior of films and foams, as well as the optical properties of thin liquid layers. In this paper, we focus on the dynamics of material flow within the soap film, which results in fascinating, extremely detailed patterns. This flow is characterized by a complex coupling between surfactant concentration and Marangoni surface tension. We propose a novel chemomechanical simulation framework rooted in lubrication theory, which makes use of a custom semi-Lagrangian advection solver to enable the simulation of soap film dynamics on spherical bubbles both in free flow as well as under body forces such as gravity or external air flow. By comparing our simulated outcomes to videos of real-world soap bubbles recorded in a studio environment, we show that our framework, for the first time, closely recreates a wide range of dynamic effects that are also observed in experiment.

In the computer graphics community, it is now widely understood how films, bubbles and foams form, evolve and break. On the rendering side, it has become possible to recreate their characteristic iridescent appearance in physically based renderers. The main parameter governing this appearance, the thickness of the film, has so far only been driven using ad-hoc noise textures [Glassner 2000], or was assumed to be constant. The resulting renderings appear largely plausible but static, as they lack the rich dynamics known from real-world soap films.
With this paper, we aim to close this gap in order to achieve greater realism. We do so by contributing a chemomechanical framework targeted specifically at simulating the rich and detailed microscopic flow on spherical soap bubbles. Our framework employs a leading-order approximation for the soap film dynamics developed by [Chomaz 2001;Ida and Miksis 1998a]. A soap bubble is modeled as a two-dimensional flow on a static spherical surface with two associated scalar fields: the film thickness and the soap concentration. We are able to show that this state-of-the-art model, paired with a custom solver, is capable of expressing the intricate flows found on real-world soap bubbles ( Figure 1) under the mutual influence of mechanical stress, film thickness and surfactant concentration as well as body and surface forces like gravity and air friction. Our simulation is performed on a staggered grid, using finite differences in space and time. An advection scheme based on BiMocq 2 [Qu et al. 2019] minimizes numerical dissipation in order to prevent fine details from washing out over time. The resulting thickness maps are presented in real time using a custom, very efficient polarization-aware spectral rendering scheme. Besides the underlying physical model, our framework is enabled by the following key contributions: • We propose a novel advection scheme for vector and scalar quantities in spherical coordinates. Our scheme, which constructs a local coordinate frame aligned with the direction of the flow, is unconditionally stable and maintains continuous behavior even near and across the poles. • We propose an implicit update step for the soap concentration, which avoids solving a stiff indefinite system and enables the use of significantly larger time steps when applying body forces. • We introduce a novel real-time shader that is designed to reflect spectrally and polarization-dependent effects under environment lighting in a physically accurate manner. This even holds for contributions which are reflected multiple times in a spherical bubble. • We investigate the influence of material parameters, geometric scale and external forces on the flow, and compare our results to real-world examples captured under lab conditions.

RELATED WORK
With their tendency to evolve into minimal surfaces, soap films embody a fundamental mathematical and physical principle in a way that is immediately relatable and fascinating to experts and laypersons alike. Consequently, they have inspired a large body of research in mathematics, physics and materials science. Some mathematicians even went so far as to use them as analog computers to solve mathematical minimization problems [Isenberg 1978]. Within the computer graphics community, the geometric properties of minimal surfaces as well as the formation, evolution and destruction of films, bubbles and foams have inspired a large number of groundbreaking works [Da et al. 2015;Glassner 2000;Ishida et al. 2017;Iwasaki et al. 2004;Kim et al. 2015;Ďurikovič 2001;Zhu et al. 2014]. Besides their geometry, the beauty of soap bubbles also stems from their dynamic iridescent patterns. The chaotic mixing, highly nonlinear vortices and turbulence of the fluid flow are not only visually interesting, but they have also led to a body of scientific work that is just as varied and colorful as its subject of study. Examples range from the visualization and study of 2D flow [Gharib and Derango 1989] or using soap bubbles as a small-scale surrogate model for planetary atmospheres [Meuel et al. 2013;Seychelles et al. 2008] via the visualization of sound and music [Gaulon et al. 2017] to using soap films as volatile display surfaces [Ochiai et al. 2013]. In computer graphics, the simple optical effect behind the characteristic iridescent colors (thin film interference) has long been understood and used [Belcour and Barla 2017;Glassner 2000;Iwasaki et al. 2004;Jaszkowski and Rzeszut 2003;Kneiphof et al. 2019;Smits and Meyer 1992;Sun 2006]. However, the film thickness, which is the main governing parameter besides the liquid's refractive index, has rarely been driven by proper physical simulation. While some works Sethian 2013, 2016;Zhu et al. 2014] have coupled thickness in their models, they use it more or less as an intermediate variable that influences the macroscopic motion and serves as a bursting condition. Most recently, such models have been equipped with the ability to propagate turbulent flow across Plateau boundaries [Ishida et al. 2020]. In contrast, our goal is to simulate the rich and detailed dynamics of the microscopic flow that is observable through thin film interference, while staying as close as possible to a state-of-the-art physical model. We turn to the fluid mechanics and physics communities, where several comprehensive models for soap film flow have been devised. Chomaz et al. [2001] derived a highly accurate model for the dynamics of a flat soap film, which is based on the asymptotic lubrication theory, assuming the thickness of the film is small compared to its lateral extent. The main contribution of that work is on the construction of similarities between soap film flows and compressible fluid flows in a planar, two-dimensional domain. The model provided by Ida and Miksis [1998a] is in principle capable of expressing general and time-varying three-dimensional soap films. In order to simulate flow using their model [Ida and Miksis 1998b], they employ a pseudospectral Chebychev spatial collocation and restrict the solution to the one-dimensional axisymmetric case.
Since the pioneering work by Foster and Metaxas [1996], physicsbased animation of fluids has been an important topic in computer graphics due to its wide range of applications for capturing effects of smoke , free surface flow Fedkiw 2001], or fire [Nguyen et al. 2002]. We note that even lubrication theory, which forms the foundation of our framework, has at least once before been employed by other members of the computer graphics community: [Azencot et al. 2015] used it to simulate the flow of thin liquid films across solid surfaces. Regarding the geometric discretization, various possible choices exist [Ando et al. 2015;Bridson 2015;Ferstl et al. 2016;Ihmsen et al. 2014;Macklin and Müller 2013]. Grid-based Eulerian simulation of fluid remains popular and widely adopted due to its superior efficiency and versatility, despite its well-known problems in numerical diffusion. In graphics, the semi-Lagrangian advection scheme presented by Stam [1999] builds the foundation for many more advanced future developments on Eulerian fluids, including some recent impressively successful examples [Narain et al. 2019;Qu et al. 2019;Zehnder et al. 2018]. Indeed, the advection equation is such a mathematically simple, yet numerically challenging equation that consistently draws a lot of attention. The difficulty is pronounced at an even higher level when one tries to solve for fluid dynamics on a spherical geometry [Hill and Henderson 2016;Yang et al. 2019] due to the notoriously difficult "pole singularity problem" [Randall et al. 2002]. In this paper, we look into an even more challenging scenario where we need to efficiently and robustly advect multiple physical quantities related to the chemomechanical physics on a soap film.

PHYSICAL MODEL
The mechanical properties of soap solutions are characterized by the interplay of water and the soap dissolved in it. Soap molecules, which have a polar (hydrophobic) and a nonpolar (hydrophilic) end, tend to settle at the water surface, so that their hydrophobic part can avoid the contact with water. As a result, the soap concentration at surfaces is usually much higher than in the bulk fluid. Soap further acts as a surfactant, i.e., the presence of soap molecules reduces the surface tension of the fluid. When the distance between soap molecules at the surface increases, surface tension increases accordingly. By adding soap to water, it becomes possible to make bubbles that can last several seconds to minutes, since the surface tension prevents them from bursting. The resulting structure of soap films consists of three layers [Couder et al. 1989]: two water-air  Couder et al. [1989]). A thin layer of liquid (thickness 2 ) is centered around a macroscopically defined surface. The polar chemistry of surfactant molecules causes them to concentrate at the liquid-air interface. We assume the number of soap molecules in the bulk fluid (here marked in gray) to be negligible. Interference color of light reflected by a thin layer of dielectric as a function of film thickness and incidence angle. Shown are simulations with a spectral resolution of 5nm for materials of different refractive index, including the most relevant material for our purpose, water. Colors are scaled so that white corresponds to a reflectance of 100%. We refer the reader to established literature [Belcour and Barla 2017;Glassner 2000] on how to compute these colors. Somewhat counterintuitively, we note that the optical path difference decreases with increasing angle. Therefore, under oblique observation it takes a thicker film to produce the same color.
interfaces populated by soap molecules and a thin layer of bulk fluid in between ( Figure 2). The thickness of a soap film is usually around 1µm, which explains the colorful interference between light reflected at the two interfaces ( Figure 3). As a soap solution's refractive index is only weakly affected by the soap concentration, the two dominant influences on the color of the film are its large-scale geometry and hence the viewing angle, and the spatially varying thickness of the fluid layer. For the purpose of this paper, we assume the shape of the bubble to be fixed. This leaves material transport within the film manifold as the main source of texture. In order to recreate the intricate dynamics found in realworld soap films, we have to understand, model and simulate this flow.

Governing equations in 3D
We formulate the fluid flow in terms of the incompressible threedimensional Navier-Stokes equations, where ∇ is the nabla operator in three dimensions, u is the fluid velocity, is the Cauchy stress tensor [Batchelor 1967, Ch. 1.3], is the mass density, and f represents body accelerations such as gravity and air friction. At a film surface, the stress condition applies as [Couder et al. 1989;Ida and Miksis 1998a] · n = (2C − )n + ∇ s , where n is the outward normal vector at the respective surface, the air pressure, the surface tension, and 2C = −∇ · n is twice the mean surface curvature. (In general, the values of these quantities differ between one surface of the film and the other.) The 2D gradient operator ∇ s within the surface acts on a scalar field Φ as The surface tension depends on the surfactant concentration Γ, i.e., the concentration of soap molecules on the surface, where is the surface tension of pure water, and accounts for the elasticity of the film. In the small concentration range, is considered to be constant [Couder et al. 1989].
Finally, the surfactant concentration Γ is driven by the advectiondiffusion equation with being the diffusivity for surfactant molecules.

Thin-film analysis and governing equation on spheres
In the following, we will restrict ourselves to spherical bubbles with radius . For small (centimeter-sized) bubbles, this approximation is reasonable. It is thus convenient to parameterize the problem in spherical coordinates ( , , ). The fluid velocity u then reads as We further note that the thickness of a soap bubble is very small compared to its lateral extent. Using lubrication theory [Oron et al. 1997;Reynolds 1886], we reduce the three-dimensional problem to a two-dimensional one. The extent of the film along the normal direction (the thickness) is introduced as a variable rather than a third dimension of the simulation domain. Suppose the inner and outer sides of the bubble are symmetrically deformed with half thickness to either side ( Figure 2), then the kinematic condition (see supplemental material) describing the interaction between the time and spatially varying film thickness = ( , , ) and the velocity u at the interface = ± can be written as With the mean half thickness 0 and the expansion parameter = 0/ , we non-dimensionalize the variables as where is the characteristic velocity and Γ 0 is the mean surfactant concentration. We substitute these variables in Equations (1), (2), (4), (5) and (7), expand u, Γ and asymptotically with a power series, and drop all terms except those with leading order of (see [Chomaz 2001;Ida and Miksis 1998a] for more details). Note that our non-dimensionalization is adapted to spherical coordinates and thus differs from the literature examples. The nabla operator ∇ now only acts within the surface, where e and e are the respective basis vectors. The governing equations are thus reduced to where = Γ 0 / 0 2 is the Marangoni number, Re = / is the Reynolds number, and are the dynamic viscosity and the mass density of the soap solution, respectively. The thermodynamic quantity =¯combines gas constant¯and temperature [Couder et al. 1989]. ′ = / is the scaled diffusivity, and the total derivative The vector V = ( , ) ⊤ represents viscous terms including second order terms as The complete terms are provided in the supplementary document, Equation (84). For readability, we drop the primes from this point onward. Further, within the scope of this paper, we assume that the soap molecules are not diffusive ( = 0) and that viscosity can be neglected. Readers interested in viscous effects may refer to Section 7.1. In spherical coordinates, the total derivative of u can be written as and its divergence as The derivatives of a scalar field Φ (which could either be the thickness or the soap concentration Γ), are Fig. 4. Under the influence of gravity and with surface tension as opposing force, bubbles assume an equilibrium state where the film thickness gradually increases from top to bottom. This is reflected by the horizontal fringe pattern observed on real-world bubbles.
There are mainly two contributions for the evolution of soap film thickness and surfactant concentration Γ: they are passively advected by the flow field [Yang et al. 2019], but also affected by inflow or outflow as expressed by the divergence terms in Equations (10b) and (10c). Unlike the full 3D incompressible flow (Equation (1b)), the 2D flow within a thin film behaves like a compressible, elastic medium.

Surface and body forces
The most important forces governing the motion of a soap film within its manifold are surface tension, gravity and air friction. The interaction between gravity and surface tension will cause thinner films to float upwards and thicker films to drop downwards, so that a soap bubble (or film) at its equilibrium state always assumes a wedge shape ( Figure 4). Let be the gravitational acceleration scaled by 2 / . If we assume the north pole of the bubble to be pointing upwards, then the gravity vector g in the spherical coordinate system (e , e ) is g = ( sin , 0) ⊤ .
As a soap bubble is very thin, the film is easily set into motion by tangential air flow. If the surrounding air is still, it damps the flow motion. For simplicity, we assume a linear Stokes drag f air = ( Cr / )(u air − u), with Cr being the drag coefficient. Taking gravity and air friction into account, Equation (10) becomes Interestingly, from Equation (17a), we observe that surface forces as surface tension and air drag are divided by the film thickness, so that thinner films are more easily driven into motion, whereas body forces like gravity act throughout the volume of the body and thus do not depend on the thickness.

METHOD
In this section we develop novel spatial and temporal discretization schemes for the governing equations. In particular, via using a staggered spherical grid (Section 4.1), we develop an unconditionally stable advection scheme (Section 4.2) that can smoothly propagate flow across the poles, as well as a projection-like implicit solver for handling chemomechanical forces (Section 4.3).

Spatial discretization
We discretize the spherical domain with a staggered grid [Yang et al. 2019], where velocities and scalar quantities are stored at different locations (illustrated in Figures 5 and 6). This allows the accurate evaluation of the concentration gradient ∇Γ and the velocity divergence ∇ · u using central differences without the formation of checkerboard patterns [Bridson 2015, Ch. 2.4]. Also, a regular grid discretization makes it possible to build a symmetric positive definite linear system that can be solved using the conjugate gradient method.
The velocity vector is split in two components, and , which are located at the center of the cell boundaries, whereas the concentration Γ and the thickness are sampled at the cell center. Assuming the staggered grid consists of × cells, then the dimension of , Γ and is × , and the dimension of is ( − 1) × . Quantities that do not lie exactly on the respective grid points are bi-linearly interpolated between neighboring grid points (Figure 7). Special care needs to be taken when sampling near the poles, as neighborhood relations reach across the pole. At this point, both e and e experience a sign change, so velocity samples drawn from across the pole have to be negated.

Advection
Taking an operator splitting approach, we first solve the material derivative ( / ) parts of u, , and Γ in Equation (17)    velocity field through advection. Afterwards, we treat the remaining force terms in treated in a separate step (Section 4.3).
On the spherical domain, the pure advection of a time-dependent scalar field Φ(x, ), x = ( , ) ⊤ , along the velocity field u(x, ) = ( (x, ), (x, )) ⊤ can be written as the initial value problem [Pironneau et al. 1992 We seek to evaluate the advected quantity at a grid point x = ( , ) ⊤ . In keeping with standard practice in fluid simulation, we introduce a virtual particle that passes the grid point at time , and trace it backward in time. This results in the time-dependent trajectory (x , ) ( ) for the seed point (x , ) and time parameter < . Substituting this trajectory into Equation (18) [Yang et al. 2019] (e) Ours (order-1, global) (order-3, global) (order-2, global) (order-2, aligned) where the finite difference corresponds to a discretization of the time domain around with step size Δ . A particle that undergoes pure advection experiences Φ as being constant. We exploit this property to approximate the value for Φ( (x , ) ( )) using a sample taken a step of Δ backward through the flow field. A numerical integration step like Euler or Runge-Kutta can be used to obtain Φ( (x , ) ( − Δ )).
Advecting vectors on a sphere is much more challenging. As e / and e / are not equal to zero (Appendix A) as they would be on a Cartesian grid, Equation (12) contains the additional terms − 2 /tan and /tan that are not present in Equation (18). Yang et al. [2019] advect vectors in a scalar-like manner with ( , /sin ) ⊤ , and treat the two additional terms as body forces in an additional backward Euler intergration step. Their method does not work well in practice when additional force terms are present. Instead, we are able to perform unconditionally stable vector advection in a single step. We note that at the equator = /2, where tan −1 = 0, these two extra terms vanish. There, the advection of a vector field falls back to the known case of advecting a scalar field. We exploit this insight in order to obtain this desirable property anywhere on the sphere. By constructing a local coordinate system at each grid point that aligns with the velocity, we can treat the quantities there as if they were on the equator.
For each grid point x = (x , ) ( ), we draw a great circle on the sphere that passes through this point and is tangent to the velocity u at this point (see Figure 9). We denote the unit vector in u-direction to beû, and the unit vector pointing from the sphere center to x to beŵ, and the binormal unit vector to bev =ŵ ×û. Then the great circle can be parameterized by i.e., changing the arc parameter translates the current grid point x = (0) =ŵ back or forward in time with unit velocity. Following u backward for a time step Δ results in a change in the arc parameter of −Δ = −∥u∥Δ , taking us to the point This results in the following single-step advection scheme ( Figure 9): (1) Evaluate the velocity u and establish the great circle in its direction.
(2) Perform an interpolated lookup of the velocity u ′ at (−Δ ) using the technique described in Section 4.1.
(3) Decompose u ′ into tangent and binormal components ′ and ′ with respect to the great circle at this point.
The main distinction of this approach to existing work is in the choice of coordinate frame. Rather than operating in the global spherical coordinate system, it locally creates an orthonormal coordinate frame that is defined by the velocity vector and therefore data-aligned. Since the method only interpolates values from the last step, it is unconditionally stable.
When higher accuracy is desired, the same principle can also be used to implement higher-order schemes involving multi-step updates. As an example, we lay out the construction of a 2 nd -order "half-step" scheme ( Figure 21) in Appendix B, which was used to generate the results shown throughout this paper.
By design, our advection scheme does not produce artifacts when advecting quantities (scalars or vectors) around the poles, see Henderson 2016] to staggered grids, where boundary conditions are introduced at both poles to remove singularities. However, according to [Hill and Henderson 2016], their method is not free of artifacts: small disturbances will appear near the pole due to variation of grid spacing (the 1 /sin term becomes prohibitively large near the pole). By constructing a local frame in our advection scheme, the 1 /sin term always takes the value 1, so the variation in grid spacing does not affect the advection. Moreover, the so-called geometry term in [Hill and Henderson 2016;Yang et al. 2019] that is caused by coordinate orientation changes in curvilinear coordinate systems (see Appendix A) is in our case implicitly handled when transforming from global to local coordinates, and does not need to be solved separately. By performing our global-to-local coordinate transformation everywhere on the sphere, we show that all the discretization points on the sphere can be treated equally and there is no need for any special pole treatment. We expect this method to perform equally well as on a Cartesian grid.

Preserving details.
The interpolated look-up in our advection scheme causes numerical diffusion known from all semi-Lagrangian methods. This is acceptable for u and Γ, since both PDEs include diffusive terms (see Equation (10)). However, the transport equation for is non-diffusive. Therefore, it is important to prevent highfrequency details in the film thickness from blurring out over time.
To achieve this goal, we make use of BiMocq 2 [Qu et al. 2019] and extend the method to spherical coordinates. Essentially, BiMocq 2 keeps a backward mapping which maps a spatial point x( ) back to its position at the initial time 0 , as well as a forward mapping which maps a spatial point at the initial position x( 0 ) to its current position at = . Instead of repeatedly blurring the features from the last time step, we acquire the initial state directly from the backward mapping, which corresponds to the pure advection part / . The additional changes (− ∇ · u in our case) are accumulated along the forward mapping and added to the acquired value from backward mapping. At each simulation step, both mappings are updated via the advection method in Section 4.2, and the coordinates are interpolated using spherical linear interpolation. When the distortion between forward and backward mapping becomes too large, both mappings are re-initialized. This, however, introduces a small amount of numerical diffusion. We found that re-initializing both mappings when the distortion is larger than /128 provides a good trade-off between sharpness and noise. For other implementation details, such as error correction, we refer the reader to the original paper.

Concentration splitting
After solving the advection part / = 0, we now deal with force and divergence terms in the right hand sides of Equation (17). Soap film exhibits elastic properties similar to a mass-spring system. Solving such systems using an explicit time integrator would require prohibitively small step sizes to achieve a stable simulation. Instead, we construct a projection-like implicit system for Γ and u. To keep the system linear, we still solve for explicitly.
We temporally discretize the continuous equations into where Γ * , * , and u * denote the respective quantities after applying the advection step. First, we solve for Γ by rewriting Equation (24a) and applying divergence to both sides of the result, Afterwards, we combine Equations (24b) and (25) and eliminate ∇ · u, Finally, we express this linear system as a sparse matrix (Appendix C), solve it for Γ using a preconditioned conjugate gradient method [Naumov et al. 2015], and update u +1 and +1 using Equations (24a) and (24c). Note that this system is strictly symmetric positive definite unless Γ * approaches infinity, in which case we end up with a Poisson equation that is closely related to pressure projection in nearly incompressible mixed finite elements. Our implicit treatment of the concentration based on its evolution Equation (24b) allows us to take significantly larger time steps compared to what could be permitted when treating the forces explicitly. Our method also avoids solving an extremely stiff indefinite system as in standard Newton-based elasticity solvers [Gast et al. 2015].

Implementation and runtime performance
We implemented our method using CUDA and AmgX [Naumov et al. 2015], and executed it on an NVIDIA GeForce GTX 1080 graphics card. A typical resolution for our simulation grid is 1024×2048 with a step size of 0.002 s. At this setting, a single time step typically takes 16-17 conjugate gradient iterations and 1.1 s to execute. The bulk of the compute time is spent on divergence/force calculation and advection with 75% and 25%, respectively.

SOAP BUBBLE RENDERING
The iridescent effects produced by thin films are wave-optical effects that arise from constructive and destructive interference of light waves. To compute the light power that is reflected or transmitted when interacting with a thin film, we need to consider the corresponding complex amplitudes of the electromagnetic wave. This section describes a real-time renderer for spherical soap bubbles that is specifically targeted at the correct handling of polarization and spectral sampling.

Thin film model
We follow the modeling of Belcour and Barla [2017] for the reflection and transmission through a single thin film layer. In our case, the light interacts multiple times with the soap bubble surface on a given light path, so their analytical solution to the spectral integration does not directly apply here. After each reflection, we would have to retract to an RGB representation of the light, introducing an error that increases with each interaction with the soap bubble. Therefore, we numerically integrate the compounded reflectance and transmittance for each light path over the wavelengths , which we sample at 5 nm intervals, and compute the fractional light power that is carried by each light path for each wavelength.
Material model. The surface of the soap bubble is modeled locally as a thin film with parallel interfaces, sandwiched between two layers of air with refractive index = 1. The refractive index of the soap water in between is = 1.33. For a single interaction of the light with the soap bubble, the film thickness 2 is assumed to be constant. For a given light direction i and surface normal n, we define the angle of incidence in air via cos = i · n. Upon refraction into the soap film, the angle at which the light travels is found via Snell's law: sin = sin . Since the film interfaces are assumed to be parallel, the angle at which the light is refracted out of the film at both interfaces equals , and the direction of the light transmitted through the film is uninterrupted from air to air. We assume that the light enters and leaves the film at the same location on a macroscopic scale, since the film thickness is much smaller than the lateral extent of the soap bubble.
Thin film reflectance. The polarization-dependent reflectance and transmittance are the ratios of outgoing to incoming light powers at an interface between two media. In order to compute these values for two media separated by a thin film, we consider the complex-valued electro-magnetic wave amplitudes. The complexvalued reflection coefficient and transmission coefficient describe the ratios of outgoing to incoming wave amplitudes. Since the power carried by a light wave is proportional to the square of the wave amplitude, we have = | | 2 and = | | 2 . In addition to the power ratio, they also encode a phase shift of the light wave, which leads to constructive and destructive interference when two light paths interfere and their reflection or transmission coefficients are added.
The reflectance and transmittance produced by an interaction with the film is computed for each wavelength and polarization by taking all light paths inside the thin film (shown in Figure 10) into account and accumulating their wave amplitudes. The light waves are affected by the reflection coefficients , and transmission coefficients , , defined by Fresnel's equations [Born and Wolf 1970], where and act at the air-to-soap interface and and on the soap-to-air interface. The light waves of each path are also affected by a wavelength-dependent phase shift, induced by the difference in path length between successively emitted light rays at each film interface. This difference in path length is known as the optical path difference D = 4 cos , which introduces a phase shift Δ = 2 D of a light path with respect to its predecessor. The phase shifts with respect to the first light ray accumulate linearly, such that the -th ray is phase-shifted by Δ . Summing the contributions from the infinitely many emitted light rays at the top and bottom interface yields the reflectance , and transmittance of the thin-film Since we are not dealing with total internal reflections, | | < 1 holds and taking the limit of the geometric series yields a closed form solution for each wavelength.

Soap bubble ray tracing
The wave nature of light only needs to be taken into account when computing the reflectance and transmittance for a single interaction with the soap bubble. Since soap bubbles are inherently transparent, many light paths contribute to a given view ray. To compute the light transport along these paths, we multiply the reflectance and transmittance produced at each soap bubble interaction along the path. The light transport R ( ) for the -th order light path is defined as for ≥ 1, where ( ) and ( ) are the iridescent reflectance and transmittance produced at the -th interaction with the soap bubble along the path traced backwards from the observer. The reflectance ( ) and transmittance ( ) produced by different interactions along a path change only due to differences in film thickness. The angle of incidence is constant for all interactions along a light path (see Figure 11) and the index of refraction is constant as well. Polarization. Due to the spherical geometry of the soap bubble, the incoming and outgoing light directions for all interactions lie in the same plane. This implies that the light polarization basis (the decomposition into s and p component) does not change between successive soap film interactions. Therefore, we are able to correctly handle polarization effects by first evaluating Equation (29) for each polarization direction independently, and then averaging both frames.
Tracing rays. For the computation of the light paths through the soap bubble, illustrated in Figure 11, we exploit the assumption that the soap bubble is spherical. We evaluate the first = 8 light paths through the soap bubble. Since we assume distant illumination, we do not track the world-space location of the film interactions, and only consider the relevant directions: the incoming light direction i ( ) of the -th order light path to sample the environment map, the surface normal n ( ) to sample the film thickness 2 at the -th interaction, and the (virtual) viewing direction o ( ) , which is used to compute the directions for the succeeding interaction with the soap film. Given these values for the -th film interaction, the directions for the + 1-th interaction are defined as The surface normal n ( +1) is defined by a reflection of −n ( ) at the outgoing light direction o ( ) , and the new outgoing light direction o ( +1) is then the reflection of the incoming light direction i ( +1) at the surface normal n ( +1) (see Figure 11).

Spectral integration.
To produce an sRGB color image, we have to integrate the response to the respective color-matching functions for ∈ { , , } of the color space:  [Couder et al. 1989, Figure 1(a)] in the low concentration range.
c We assume advection to be the dominant transport mechanism.

RESULTS
In the following, we perform a selection of synthetic experiments and discuss the influence of the most important parameters and variables. From various sources, we attempted to gather as realistic a set of parameters as possible. The numbers used for simulations throughout this section, as well as sources for the more exotic values, are listed in Tables 1 and 2.

Mean surfactant concentration and bubble radius
According to the momentum equation (17a), a soap film under gravity but without other sources of excitation has its equilibrium state at The other two Equations (17b) and (17c) can be rewritten as i.e., Γ / remains constant. If the simulation is started with uniform surfactant concentration and thickness, we can non-dimensionalize the variables so that Γ( = 0) = ( = 0) = 1. Combining the above two equations, this yields the solution of which is From this, we can draw at least three conclusions: • Soap films under the influence of gravity are thinner at the top and thicker at the bottom, leading to colorful bands on soap bubbles ( Figure 12). As noted in Figure 3, the film color is also influenced by the viewing angle. On a bubble, the bands are bent downwards; on a flat film they appear horizontal. • A constant ratio / will lead to the same equilibrium state (albeit through a different dynamic process). From the definitions of and (Table 2), this is equivalent to keeping /Γ 0 constant. • The larger /Γ 0 , the thinner the film is at the top, and thicker at the bottom. This expected behavior is confirmed in experiment and simulation (Figures 12 and 13).

Gravity and buoyancy
The momentum equation of a soap bubble has great resemblance with the compressible Navier-Stokes equation where the surfactant concentration Γ takes the role of pressure and the variable thickness substitutes the variable density . In fact, just as smoke with smaller density flows upwards in the air, thinner soap film regions also tend to flow upwards. This is confirmed by our observation. As thinner regions on a soap bubble flow upwards (and thicker regions downwards), they form drop-shaped "islands", and leave "rivers" behind ( Figure 14).

Air friction
Soap films are highly susceptible to air flow, and assume beautiful patterns in windy environments. When a bubble is blown, a rapidly rotating wind field is produced inside and induces a shear motion on the soap film. Advection along the air flow results in thin stripes that remain stable even after the external influence has stopped (Figure 15). See Figure 16 for a false-color visualization of the external air flow and the resulting velocity field in multiple time steps of an experiment.  [Cabral and Leedom 1993]. From left to right, the correlation of the lines shows the direction of the respective vector field, while the magnitude is encoded in the color. (a) After initialization with a noise texture, the fluid sags down under the influence of gravity. (b)-(d) Over a corridor on the surface, air gradually starts flowing and slows down again. The soap fluid follows the excitation. (e) After the air has stopped flowing, the bubble remains in a rotating motion.

Evaporation
Due to evaporation, a soap bubble exposed to air becomes thinner and thinner and eventually breaks down. Since evaporation mostly depends on the exposed surface (which is constant), we model this effect by subtracting a small constant amount of at each simulation step. Once a point on the surface reaches thickness zero, the simulation is terminated. Our model in its current form does not support the simulation of bursting bubbles. Figure 17 shows a soap bubble over its whole lifetime. Starting from a random thickness distribution, thicker regions are moving downwards and horizontal color bands. Shortly before bursting, the top of the bubble becomes very thin and exhibits a gray appearance.

Real-world experiments
To capture stills and videos of real-world bubbles under laboratory conditions, we constructed simple studio environments consisting of 1200 mm×600 mm LED panels, black theater curtain and a Sony ILCE-7RM3 system camera with a ZEISS Batis 135 mm /2.8 lens.
We use Pustefix brand soap solution for all experiments.

DISCUSSION AND FUTURE WORK
We have been able to show that our model and solver, which is fast and stable, can recreate the most prominent effects found on spherical soap bubbles in the real world. An obvious next step will be to look into more general cases, like complex film shapes or groups of bubbles. Although we focus on spherical geometry in this paper, the Fig. 17. Life of a soap bubble. The bubble was initialized with a Perlin noise texture and excited by curl-noise air flow. As water evaporates and the film becomes thinner, the colorful bands move gradually downwards and the top of the bubble fades to gray. ideas underlying our scheme are not limited to spherical domains. The model in [Ida and Miksis 1998a] is valid for general manifolds; following their derivation, one arrives exactly at Equation (10), with the definition of the differential operators adjusted to the corresponding curvilinear coordinate frame. Furthermore, the idea of our advection scheme is independent of the underlying manifold shape, as long as a proper local coordinate frame is constructed. Finally, the special force and divergence terms treatment in Equation (26) holds for arbitrary shapes and can be expressed as a sparse matrix, as long as the neighborhood of each point is well-defined. However, nice-to-have properties (such as the matrix being symmetric positive definite) might be lost in other types of meshes or grids. Consequently, our scheme should not be difficult to generalize to bubbles that are diffeomorphic to a sphere. Groups of bubbles with Plateau borders, as well as bounded films, are not manifolds and hence have to be left to future consideration. The treatment of film boundaries deserves particular attention also because a complicated mechanism called marginal regeneration [Isenberg 1978] causes the film to become even thinner at the boundaries, producing regions that flow upwards erratically (see Figure 18).

Viscous film
Soap films with larger viscosity tend to be more stable and last longer, which can be achieved by adding glycerin to home-made soap solution. Some commercial soap solution also includes additional formula to make it more viscous. The soap solution we used in (a) Viscous soap film (b) Inviscid soap film Fig. 19. Two simulations with/without viscosity term Re −1 V, otherwise under the same condition and after same frame numbers. A viscous film tends to keep its texture longer in shape and has a reduced tendency to break into fractal structures.
experiments, for example, has a dynamic viscosity of 1.2 × 10 −1 Pa s, which is about 100 times greater than that of water and thus shows different dynamics. The example in Figure 19 was obtained by adding a very basic explicit step to compute viscosity; however, this is slow and unstable and thus not included in our standard solver.

Black film
As a film keeps thinning through evaporation or gravity drag, at some point it becomes so thin that destructive interference takes place, and the film appears completely black. The thickness is then about 5 nm-30 nm [Isenberg 1978;Rücker 1877]. At such a small scale, molecular forces come into play, such as Van der Waals attraction, electrostatic repulsion, and Born repulsion. These molecular forces cause black film to be surprisingly stable and form sharply defined "islands" within the colorful film ( Figure 20). In future work, it will be interesting to include such molecular forces in an extended model.

A MATERIAL DERIVATIVE IN SPHERICAL COORDINATES
The material derivative describes the rate of change of a quantity moving with a time-dependent velocity field u. We denote a dimensionless scalar quantity in spherical coordinates as Φ( , , ).
Applying the chain rule yields The dimensionless velocities in spherical coordinates are defined as Note that we included the sphere radius in our non-dimensionalization of in Section 3.2, thus we do not need to take account of the sphere radius when taking the derivatives. Substituting Equation (40) in Equation (39) gives The material derivative of a vector quantity, for example velocity, is similarly written as where To evaluate the remaining two partial derivatives, we take a step back to look at the derivatives of unit vectors in sphere coordinates. Neglecting the radial ( -dependent) component, they are Combining Equations (40), (43), (45) and (46) B 2 ND -ORDER HALF-STEP UPDATE Based on the single-step scheme detailed in Section 4.2, we construct a multi-step scheme inspired by second-order Runge Kutta that can deliver higher accuracy ( Figure 21). First, the velocity vector sampled at the grid point (red) defines a great circle (also red) which is used to construct a local coordinate frame. Following this direction backward by a half time step, a second velocity vector (blue) is sampled. After transforming this vector back to the original point, it generates a second great circle (also blue). A full backward step along this circle takes us to the look-up location ( * ) for the value that is advected to the grid point.

C LINEAR SYSTEM
We solve Equation (26) on a grid of dimension × and cell size Δ × Δ in matrix form = b, where = Γ 0,0 , . . . , Γ −1, −1 ⊤ ∈ R , = 0,0 , . . . , −1, −1 ⊤ ∈ R and ∈ R × is a sparse, block diagonal matrix with five elements in each row, corresponding to the cell itself and its four direct neighbors. In order to make symmetric positive definite, we multiply both sides of Equation (26) by sin . Then, in each row, the entries of are given by where and denote the cell center of the respective mid cell. The right hand side is given by , = sin 1 Δ − ∇ · * , u * , + CrΔ (u air ) , + Δ * , g ,