ON THE STOCHASTIC IMMERSED BOUNDARY METHOD WITH AN IMPLICIT INTERFACE FORMULATION

. In this paper, we present a consistent and rigorous derivation of some stochastic ﬂuid-structure interaction models based on an implicit inter- face formulation of the stochastic immersed boundary method. Based on the ﬂuctuation-dissipation theorem, a proper form can be derived for the noise term to be incorporated into the deterministic hydrodynamic ﬂuid-structure interaction models in either the phase ﬁeld or level-set framework. The result- ing stochastic systems not only capture the ﬂuctuation eﬀect near equilibrium but also provide an eﬀective tool to model the complex interfacial morphology in a ﬂuctuating ﬂuid.


1.
Introduction. The problem of fluid-structure interaction is of central interests in many areas of sciences and engineering ranging from the aerodynamics of flying objects, colloidal suspension in chemical reactions, to cell migration in biological environment. Such interactions may involve phenomena at many different length and time scales. While the microscopic fluctuations are generally very small on the macroscopic scale, they play essential roles in many important applications such as those in materials sciences and fluid mechanics [6,23,25,29]. An effective approach to model such small scale fluctuations has been the introduction of random noises in the macroscopic field equations. Here, we are particularly motivated by recent works on the phase field models for the vesicle-fluid interactions [12,13]. Indeed, the studies of lipid bilayer vesicles immersed in a fluid environment have attracted much attention in recent years. It is widely known that the effects of thermal fluctuations in such a system are very important for various biological processes and functions [29]. This raises the question on how to properly incorporate the fluctuation effects into the phase field-fluid interaction model or other models based on the implicit interface representations. In this work, we show how this issue can be addressed based on the stochastic immersed boundary method [2,20] using an implicit surface representation as in [12,13]. The resulting equations are systems of stochastic partial differential equations involving the fluid variables (velocity and pressure in the primitive variable formulation) and the phase field or level set variables characterizing the motion of the immersed interface.
We note that stochastic PDEs have been used widely to study the effect of random fields and a comprehensive survey of the literature is beyond the scope of this 374 QIANG DU AND MANLIN LI manuscript (we refer to [10] for some general discussions and references). On works related to the fluid-structure interactions, and in particular the stochastic immersed boundary methods, we refer to [2,20,21] and the references cited therein. There are also related works concerning particles/structures immersed in an random but prescribed fluid field, see for example, [27]. The use of implicit surface representation for an immersed boundary has also been studied by many researchers before, for instance, the level set formulations [3,7,8,9] and the phase field or diffuse interface formulation [1,4,12,16,19,32] are two popular choices. It turns out that using the principle of least actions [5], one may derive the different formulations of the immersed boundary method, including both explicit and implicit interface representations within the same framework as we briefly outline later. Then, one may mimic the derivation of the stochastic immersed boundary method (SIBM) to derive the corresponding stochastic systems of equations to account for the near equilibrium fluctuation effects when the interface is implicitly represented by either a level-set function or a phase field variable. Based on the fluctuation dissipation theorem, we then show that the phase field or level set functions are governed by a transport equation, while the random fluctuation term appears only in the momentum equation for the fluid velocity field with a suitable functional form. Moreover, such forms of stochastic equations can also be consistently derived when a mollified velocity field is used in the transport equation, with appropriate modifications made to the contributions in the momentum equation due to the interaction with the interface. The resulting stochastic system of equations, which we refer as the stochastic implicit interface model (or SIIM for short), can be used to develop effective numerical algorithms for simulating the fluid-structure interaction with complex interfacial morphologies, in particular, for handling, in a natural way, the change of topology due to stochastic fluctuations. The key models are summarized into propositions, which should not be interpreted as proven statements but as a way to highlight the equations. The rigorous analysis on the well-posedness of the resultant stochastic systems will be provided in the future works. While the framework of the SIIM established here works in general, we present several examples related to some of our past works as illustrations. We also, whenever possible, make comparisons and draw connections between the different models, either in deterministic or stochastic forms. Additional comments on future works are provided at the end of the paper.
2. The stochastic immersed boundary method. The immersed boundary method can be formulated to model the interaction between an incompressible viscous fluid and an immersed structure represented by discrete particles. The typical equation of immersed boundary method [28] is given by where X i is the position of the i th particle, f i is the force interaction between the i th particle and the fluid, and ζ r is a delta function for r = 0 or a regularized version parameterized by some positive constant r > 0. Typically, r may represent the size of a blob around the particle position in the latter case. In the derivation of the stochastic immersed boundary method (SIBM) [2,20], a potential energy E({X i }) is defined as a function of the particle position vectors {X i }, and the forces are given by For fluids in the low Reynolds number regime, the Navier-Stokes equations may be simply replaced by the Stokes equation. To model the stochastic effects, noise terms f and {g i } are first added into the equations to get: Then, the equation (2) is linearized around the equilibrium solution of the corresponding deterministic equation and the resulting linear system is subsequently expanded in the Fourier space. The correlation structure of the solution to the linearized equation in Fourier space can be explicitly computed. By imposing that the solution obeys the Boltzmann distribution as the dynamics reaches the thermal equilibrium, the proper forms of noises f and {g i } are obtained. If one views the index i as a continuous parameter, then the particles forms an immersed interface and the summation in the momentum equation can be simply replaced by an integration over the interface. For details, we refer to the original derivations given in [2,20].
3. Fluid-structure interaction with an implicit interface representation. Many interfacial problems can be effectively modeled by an implicit interface representation such as the level set methods and phase field or diffuse interface methods [1,26]. For instance, phase field has been a convenient method for modeling many interfacial problems. In general, a phase field method uses a smooth phase field function to distinguish two regions, or phases, in the physical domain. For instance, the phase field function takes value 1 in one phase and −1 in the other. The interface between the two phases is characterized by a thin transition layer around the zero level set of the phase field function, thus leading to a diffuse interface description.
An advantage of using the phase field approach is that, through an implicit interface formulation, it avoids explicitly tracking or constructing the interface between the two phases. The dynamics and energetics of the interfaces are encoded in the dynamics and variational calculus of the phase field functions. Phase field method is popular in many applications such as in the modeling of solidification dynamics, microstructure evolution, and vesicle dynamics. As a particular motivation to the current work, we recall a recent example given by the phase field vesicle-fluid interaction model for the deformation of a bio-membrane vesicle in fluid [13]. A vesicle membrane Γ, formed by lipid bilayer, has an equilibrium configuration often characterized by the minimizer of the elastic bending energy. Although more general forms of the bending energy can be considered, for illustration, we focus on the following form [30]: where Γ is the surface of vesicle membrane, H is the mean curvature of Γ. c 0 is the spontaneous curvature and κ is the bending modulus. In our earlier works, a phase field function φ is used to represent the vesicle configuration as a labeling function.
φ takes value nearly +1 inside the vesicle membrane, and −1 outside. The width of the thin transition layer is characterized by a small positive parameter . The elastic bending energy (3) can be approximated by [14,15], where Ω is the computational domain. The deterministic phase field vesicle fluid interaction model is derived through the least action principle in [13]. For the fixed time interval [0, T 0 ], the least action principle asserts the fluid flow particle trajectory x(t, α) minimizes the action functional, The fluid velocity field u(t, x) is defined, in an implicit fashion, by And the phase field function is transported by the fluid flow via which is implicitly determined by x(t, α) as well. Via suitable variations to the trajectory x = x(t, α), the deterministic phase field-fluid interaction equation [13] can then be derived as: For convenience, ∇E[φ] is used here and throughout the paper as a simple notation for the functional derivative of E[φ] with respect to the variation in φ. We allow the functional E to depend on not only φ but also its spatial derivatives, so ∇E[φ] should also account for variations to the derivatives of φ as well (see examples provided later). Note that the general form of the equation (7) is independent of the specific choice of the energy E[φ]. The derivation is thus applicable to more general settings where the interfacial energy E = E[φ] may represent other physical energies such as the surface tension or it may be given in a level-set formulation rather than a diffuse interface representation.
In this work, we present a stochastic implicit interface model which is an analog of the stochastic immersed boundary model, but with an implicit interface representation. The stochastic systems can be specialized to either the level set or the phase field representations of the immersed interface. These models, due to the implicit interface formulation, have their advantages in dealing with more complex interfacial morphologies such as those involving topological changes.

Stochastic immersed boundary with implicit surface representation.
In the immersed boundary method, the immersed structure is explicitly described by the positions of particles {X i }. And the potential energy is a function of {X i }.
In the implicit interface model, we use a continuous labeling function to characterize the interface configuration and the potential energy is in functional form. The deterministic equation (7) formally has a dissipative energy law, , and the Hessian of gives a stable steady state of the deterministic equation (7). Through the calculation of the correlation structure of the linearized equation, we can derive the proper form for the noise to be incorporated in the stochastic system based on the fluctuation-dissipation theorem. The linearization makes the derivation simpler by allowing a linear algebraic (matrix) notation, although the conclusion holds for the nonlinear system as well. For the sake of illustration, we assume a periodic boundary condition for both the flow variables and the interface marker with the computational domain Ω to be the unit cell. Let L 0 be the linearized operator of ∇E at φ 0 . As an illustrative example, if we take then we may explicitly have By linearizing the equation and introducing the noise terms f and g, we formally have the following coupled stochastic system: where ψ = φ − φ 0 , f and g are two Gaussian noises which are assumed to be delta correlated in time. The more specific forms of f and g are given below. Now by expanding (8) in Fourier space with k denoting the wave number and usingû k andψ k to denote the Fourier coefficients of u and ψ respectively, we then get where P k = I − k ⊗ k |k| 2 is the projection to ensure incompressibility of fluids [31]. In the Fourier variables, we have with {σ l,k } determined by the linearized operator L 0 , and, Moreover,f k andĝ k are the corresponding Fourier coefficients of the noise terms f and g, i.e.
f k andĝ k are finite dimensional Gaussian white noises. Their correlations are to be determined in Proposition 1 and 4 later. Let a(t) = (û,ψ) T . The linearized equation (9) can be abbreviated as an infinite system of ODEs, d dt where n(t) = (P kfk ,ĝ k ) T and the matrix operator D has the form with its entries defined by: Then D 12 = P · Π · Σ and D 21 = −Π.
It follows from the fluctuation-dissipation theorem [20,17] that, the equilibrium solution obeys the Boltzmann distribution. The probability density formally is given by where Z is a normalization constant, K b is the Boltzmann constant and T denotes the temperature.
To be more precise, we have the following observations. 1. The velocity field is real, thusû k =û −k . Moreover, due to the incompressibility of the fluid,û k · k = 0. It implies that the probability distribution of u k lives on the orthogonal complement of the vector k, rather than the whole space. 2. ψ is also a real function, henceψ k =ψ −k .
Immediately, we see thatû andψ are independent. Therefore, we can assume that, when the dynamics reaches the thermal equilibrium, the correlation matrix of the solution is given by where < · > represents the expectation, Direct calculation shows that (C 11 ) k,k = δ k,k P k and C 22 = Σ −1 . For illustrative purposes, let us consider the case of two space dimension though the derivation works in three dimension as well. Let r k and s k represent the real and imaginary parts ofû k respectively. Because of the independence implied by the Boltzmann distribution andû k =û −k , we get For k = 0, let |r k | = r and k |k| = (cos θ, sin θ), then, We can similarly calculate < s k ⊗ s k >= K b T 2 · P k . Therefore <û k ⊗û k >= K b T · δ k,k P k .
That is, (C 11 ) k,k = δ k,k P k . For C 22 , it can be computed similarly without much difficulty.
While the derivation uses a two dimensional form, one can extend the calculation to three dimension as well. We omit the details.
It follows from (10) that, where {·} * denotes the adjoint operation (transpose of the complex conjugate) and the last term N is the Ito's correction. Recall the assumption that the noise is delta correlated in time, N satisfies < n(t) ⊗ n(t ) >= N δ(t − t ). Setting t → ∞, we get from the equilibrium distribution that Taking N in block matrix form, Now, N 22 = 0 implies that the noise g = 0. To highlight the conclusion of the above derivation, we summarize the discussion into the proposition below. Proposition 1. The coupled system of equations describes the dynamics near the mechanical and thermal equilibrium. That is, φ should be purely transported via velocity field without any noise perturbation, and a proper noise f should be added to the fluid momentum equation. When expanding in Fourier space of the spatial variables, the Fourier coefficients {f k } of f satisfy for any k, k and t, t .
We note that the form of the noise term f given in the Fourier space can also be represented in more familiar terms [22], when expressed as ∇ · S with a random stree tensor S = (s i,j (t, x)) having a space time white noise correlation sturcture [11]. Moreover, in the above discussion, the energy functional can take on a very general form. We now discuss several examples here as illustration.
Example 1. We first consider a stochastic model for fluid interface interaction under the surface tension.
In the study of the surface tension based on a phase field function the energy E[φ] is defined as [6,13,18,32], which represents the surface area of the interface. A constant surface tension is assumed. The corresponding force density term is given by Substituting the above into the system (11), we get a stochastic phase field fluidinterface interaction model with the interfacial energy being governed by the surface tension. While written formally as a force term, ∇E[φ]∇φ should be best interpreted as an approximation to the extra stress due to the interface deformation. Notice that the momentum equation can also be rewritten as where the nonlinear term −1 (φ 3 − φ)∇φ is absorbed into the pressure term. Example 2. We have briefly mentioned the phase field approximation of the bending elastic energy in the studies of vesicle membranes [13]. The corresponding phase field energy is given by where the first penalty term involving M 1 is used for preserving the surface area of the interface and the second penalty term involving M 2 is for the volume conservation. For simplicity, we set c 0 = 0. The force density ∇E[φ]∇φ calculated in [12,13] is given by where Substituting into the momentum equation, we get the stochastic phase field fluidstructure interaction model with bending elastic energy. We note that with the incompressibility condition satisfied by u, we can set M 2 = 0 in the last term of the force expression (14), since the volume constraint is automatically satisfied by the transport equation for φ. Similar to the previous example, the first term in the force ∇E[φ]∇φ corresponds to the diffuse interface approximation to the extra Willmore stress due to the bending elastic energy (see [13] for a derivation of the consistency with the sharp interface model). The penalty term for the area constraint can also be formulated via introducing a Lagrange multiplier. Effectively, its contribution is then similar to that derived from the variation of the surface tension (13). Example 3. In [7,8,9], a level set model for incompressible fluid-structure interaction is discussed. Similar to the phase field approach, a level set function φ is used to represent the elastic material in fluid. The energy is defined as where E d is the energy density, is a small positive parameter, and cut-off function η(r) = 1 2 (1 + cos(πr)) for |r| ≤ 1 and η(r) = 0 otherwise. The force density is We need to point out that the force density in [9] is derived through the principle of virtual work. That is the force density, say F , is derived from The expression of ∇E[φ]∇φ in (15) differs from F up to a gradient term, which is absorbed into the pressure term. By substituting the force density (15) into the system (11), we get a stochastic level-set fluid-interface model. Depending on how the energy density E d is chosen, it can be used to model interfacial energies of various kinds, just like the phase field counterpart. All the above three examples are special cases of the SIIM given in (11). Naturally other variations can be considered in the future.
Before we end this section, a couple of remarks are in order. First, recall that we linearized the equation around an equilibrium state of the deterministic dynamics to derive the noise form. However, the only noise contribution appears in the coupled stochastic model is in the momentum equation. Moreover, it is independent of the choice of equilibrium state. This means that potentially we could adopt this model to study the case when the deterministic dynamics may have multiple equilibrium states (even though the steady state velocity field should always be zero due to the dissipative energy law).
Secondly, the noise we have derived for the stochastic implicit interface model is very much the same as the noise derived in stochastic immersed boundary method [20]. This is not surprising since the deterministic versions of the two models also share much similarities. First, both phase field or level set functions and the particles representing immersed structure are transported by fluid flow. Secondly the deterministic momentum equations of fluid can be both derived within the framework of the least action principle. In section 5, we will comment on this matter further.

5.
General stochastic implicit interface model for transport via a mollified flow. In section 4, we have derived the proper form of the noise for phase field-fluid interaction equation. It turns out the noise term only appears in the fluid momentum equation. Generally speaking, the presence of the noise term weakens the regularity of the fluid velocity field, which leads to the difficulty of interpreting transport equation of the phase field or level set function in a suitable functional setting. As a possible remedy, we consider a mollified fluid velocity field, a technique dated back to Leray [24] and used widely in the literature, including the setting of the Immersed Boundary Methods and SIBM [20].
For simplicity, we consider a symmetric kernel ζ, that is, ζ(x) = ζ(|x|), and define φ via the mollified transport equation, Here * denotes the standard convolution operator. Despite of this change, it is still possible to derive the corresponding deterministic equation through the least action principle on a periodic domain. As we will see later, the mollification does not affect the noise form derived in section 4.
Define the mollified flow map z(t, α) by Since u is divergence free, z(t, α) is a volume preserving map. The transport equation (16) implies that φ(t, x) = φ 0 (z −1 (t, x)). For a fixed time interval [0, T 0 ], define the action functional, As before, the least action principle asserts that the correct fluid particle trajectory minimizes A[x(t, α)]. Hence the variation gives with p being the pressure, that is, the Lagrange multiplier due to incompressibility of the fluid.
To be more precise about the variation, let x ε (t, α) be a family of impressible fluid particle trajectories, where ε is the index. The variation δ means Before obtaining the main result, we first prepare a lemma for technical completeness.

QIANG DU AND MANLIN LI
Converting back to Eulerian coordinate and noticing that divδx = 0, we have Since δz(0) = δx(0) = 0, we conclude that holds for any α 0 . Now we move on to the derivation of the fluid momentum equation through the least action principle. We state the result as the following proposition.
The version of the least action principle used above does not provide the term µ u for the viscous effect. The viscous term is normally added back to the equation to provide the energy dissipation. For the low Reynolds number case, we can ignore the nonlinear term u · ∇u in the momentum equation and consider simply the linear Stokes equation with the suitable modification. Taking all these into account, we summarize the results into the proposition below, again, just to highlight the equations. Proposition 3. Consider the level set or phase field function φ being transported via the mollified flow and fluid in the low Reynolds number regime, we have the fluid velocity field and φ satisfying the coupled system of equations, We remark that the deterministic Immersed Boundary Method can also be derived within the framework of the least action principle, with or without considering mollification of the fluid velocity. The corresponding action functional is This then leads to After adding the viscosity term µ u, the Immerse Boundary Method is recovered. One can find similar calculation in section 4 when we consider how to include the proper noise term to the stochastic version of model with a mollified fluid velocity. Note that formally we can derive the following dissipative energy law by multiplying u to the first equation in (18) and ∇E[φ] to the second equation, then integrating over the spatial domain: Assume that φ 0 is a non-degenerated local minimum of E[φ]. Linearizing (18) around (u ≡ 0, φ = φ 0 ) and incorporating the noises f and g to the dynamics equations, we get We then expand the linearized equations in the Fourier space, using the same notation as in section 4. Recall a(t) = (û,ψ) T , we have The only differences from the terms in the equation (10) are 1.
where I is an identity matrix, and 2. D 21 = −Π · Y . Then we may continue the calculation further to get N = −K b T (DC + CD * ) with N being the same as before except that So we reach the same conclusion that the noise g should vanish too and f remains the same as the case without the velocity mollification. As before, we summarize this discussion into the following proposition.

Proposition 4 (Stochastic implicit interface model (SIIM)).
To capture the near equilibrium fluctuation effect in the structure-fluid interaction model with an implicit interface representation, the following system of equations can be used.

The velocity field u is incompressible
and satisfies the stochastic equation where the noise f is delta correlated in time, and in Fourier expansion of the spatial variables, it satisfies (12).
As noted before, the noise f can be written as ∇ · S with a random stree tensor S [11]. We note also that if the mollified kernel ζ, in a limiting case, is a delta function, then the SIIM equation (19)(20)(21) goes back to the case discussed in section 4 with the same form of noise satisfying the same equation (12). Again, we may elect to use either the phase field or the level set formulation for the interfacial energy term E = E(φ) in the above equations of SIIM.
We make a couple of remarks. First, while the striking similarity in the SIIM with and without the mollification is not surprising as seen from the derivation, the assumption on the periodic boundary conditions simplifies much of the technical discussion. For other boundary conditions, one may consider deriving similar stochastic models, but special care near the boundary needs to be taken when a mollification is introduced. Secondly, we considered a simplified linear Stokes equation in Propositions 3 and 4. The derivation of the noise through the linearized equation is identical for both the full Navier-Stokes and linear Stokes equations. In practice though, it is often that in the low Reynolds number regime, the simplified Stokes equation can be utilized to describe the bulk fluid motion. Computationally, simulations based on the SIIM corresponding to the Stokes equation may also be more efficient to carry out than the full Navier-Stokes counterpart. We will examine this issue in future computational studies. 6. Conclusion. In this work, working with implicit interface formulations of the fluid-structure interaction models and following the ideas used for deriving the stochastic immersed boundary method, we presented in section 4 the proper form of the noise to account for near equilibrium fluctuations in the stochastic implicit interface model. In addition, in section 5 we extended the model in a consistent fashion to more general cases when the phase field or level set functions are transported by a mollified velocity field. We illustrated that both deterministic and stochastic models can be derived using the same framework as for the case without the mollification. While in some of the detailed calculations we have used the forms associated with the two dimensional space, the models remain valid in three dimensions and can be written correspondingly for the fluid momentum equation in either the Navier-Stokes or the linearized Stokes form.
The models studied here have their advantages in dealing with more complex interfacial morphologies such as those involving topological changes, thus they could be a good alternative to model the complex morphology of fluid-structure interactions subject to fluctuation effects, especially when the change of interfacial topology might take place. We thus expect applications of our model stochastic systems in many physical and biological applications.
Theoretically, a more rigorous treatment of both the deterministic and stochastic models proposed here is also possible, including the discussion of the well-posedness for mollified flows [11]. In practice, though, it is important to develop efficient numerical schemes for our model in order to apply them to simulate interesting physical and biological processes. Such tasks will be explored in our future works.