Geometric Mixing, Peristalsis, and the Geometric Phase of the Stomach

Mixing fluid in a container at low Reynolds number— in an inertialess environment—is not a trivial task. Reciprocating motions merely lead to cycles of mixing and unmixing, so continuous rotation, as used in many technological applications, would appear to be necessary. However, there is another solution: movement of the walls in a cyclical fashion to introduce a geometric phase. We show using journal-bearing flow as a model that such geometric mixing is a general tool for using deformable boundaries that return to the same position to mix fluid at low Reynolds number. We then simulate a biological example: we show that mixing in the stomach functions because of the “belly phase,” peristaltic movement of the walls in a cyclical fashion introduces a geometric phase that avoids unmixing.


Introduction
How may fluid be mixed at low Reynolds number? Such mixing is normally performed with a stirrer, a rotating device within the container that produces a complex, chaotic flow. Alternatively, in the absence of a stirrer, rotation of the container walls themselves can perform the mixing, as occurs in a cement mixer. On occasions, however, mixing is attempted by a cyclic deformation of the container walls that does not allow for a net relative displacement of the corresponding surfaces, situations that often occur both in artificial devices and in living organisms. At the lowest Reynolds numbers, under what is known as creeping flow conditions, fluid inertia is negligible, fluid flow is reversible, and an inversion of the movement of the stirrer or the walls leads-up to perturbations owing to particle diffusion-to unmixing, as Taylor [1] and Heller [2] demonstrated. This would seem to preclude the use of reciprocating motion to stir fluid at low Reynolds numbers; it would appear to lead to perpetual cycles of mixing and unmixing. The question then arises of how cyclic changes in the shape of the containers could lead to efficient mixing. Consider a biological case of cavity flow: the stomach. In the stomach food and drink are mixed to form a homogeneous fluid termed chyme, which is then digested by the intestines. Gastric mixing is produced by what is called peristalsis: by the stomach walls moving in a rhythmic fashion. In mathematical terms, the shape of the stomach walls undergoes a closed cycle in the space of shapes during each peristalsis cycle. Obviously only shape cycles that do not require a cumulative net displacement between any two sections of the stomach can be considered. How then is this peristaltic movement of the stomach walls able to produce mixing, especially in animals in which the stomach dimensions are such that fluid inertia of the stomach contents is negligible?
The solution to this conundrum involves the concept of geometric phase. A geometric phase [3] is an example of anholonomy: the failure of system variables to return to their original values after a closed circuit in the parameters. In this Letter we propose what we term geometric mixing: the use of the geometric phase introduced by nonreciprocal cycling of the deformable boundaries of a container as a tool for fluid mixing at low Reynolds number. To exemplify how this process leads to efficient mixing, we use the well-known two-dimensional mixer based on the journal bearing flow but subject to a much-less-studied rotation protocol that satisfies the geometrical constraints of cyclic boundary deformations. We lastly show that peristalsis, besides its contribution to important biological functions such as fluid transport within individual tubular organs [4,5] or signaling throughout complex biological structures [6], fulfills its central role in gastric mixing and digestion [7][8][9] by operating thanks to a geometric phase in the stomach.

Journal bearing flow
Taylor [1] and Heller [2] used the Couette flow of an incompressible fluid contained between two concentric cylinders to demonstrate fluid unmixing due to the time reversibility of the Stokes regime. They showed that after rotating the cylinders through a certain angle, it is possible to arrive back at the initial state-to unmix the flow-by reversing this rotation through the same angle with the opposite sign, even when the angle is large enough that a blob of dye placed in the fluid has been apparently well mixed. Considering as parameters in this device the positions of the outer and inner cylindrical walls of the container specified respectively with the angles θ 1 and θ 2 from a given starting point, a geometric phase might arise from driving this system around a loop in the parameter space. In a fluid system in the Stokes regime, like ours, as inertia is negligible the motion is by definition always adiabatic and only induced by the change in the parameters: the positions of the cylinders. Therefore, any resulting phase after a complete cycle in parameters is a geometric phase. In the Heller-Taylor demonstration the parameter loop is very simple: θ 1 first increases a certain amount and then decreases the same amount while θ 2 remains fixed. This loop encloses no area, and reversibility ensures that the phase is zero. More complex zero-area loops can be constructed by combining in succession arbitrary pairs of reciprocal rotations of both cylinders, and they also lead to a null phase. We shall call these constructs reciprocal cycles. In order to consider less trivial loops, we may first note that the parameter space is homotopic to a 2-torus. Loops on such a space can be classified according to the number of complete turns that both parameters accumulate along the loop. Note also that a relative rotation of 2π between the walls brings the container to the original configuration except for a global rotation. Since we are interested in shape loops that can be achieved without a net cumulative displacement of the surfaces of the containers, we need to consider only the class of type-0 or contractible (to a point) loops.
All zero-area reciprocal loops are contractible, but there are many more enclosing a finite area. To obtain a finite-area non-reciprocal contractible loop we can, for instance, rotate first one cylinder, then the other, then reverse the first, and finally reverse the other. However, for concentric cylinders the streamlines are concentric circles; if we move one of the cylinders by angle θ, a tracer particle will move along a circle an angle that only depends on θ. Then it is obvious that the cumulative effect of moving one cylinder θ 1 , then the other θ 2 , then the first −θ 1 , and the second −θ 2 , is to return the particle to its original position: there is no geometric phase, and unmixing still occurs. But if we modify the Heller-Taylor setup and offset the inner cylinder, we arrive at what is known as journal-bearing flow. On introducing an eccentricity ε between the cylinders, this flow has a radial component. In the creeping-flow limit, the Navier-Stokes equations for the journal-bearing flow reduce to a linear biharmonic one, r 4 ψ = 0, for the stream function, ψ, and we may model this system utilizing an analytical solution (see [10][11][12] and the Materials and Methods section for further details). If we now perform a parameter loop by the sequence of rotations detailed above, we arrive back at our starting point from the point of view of the positions of the two cylinders, so it is, perhaps, surprising that the fluid inside does not return to its initial state. We illustrate the presence of this geometric phase in Fig 1 in which an example of the trajectory of a fluid particle is shown as the walls are driven through a nonreciprocal contractible loop. Journal-bearing flow has been much studied in the past [13][14][15][16], but never with contractible loops so that this geometric effect was never emphasized. This minor protocolary modification in a well established flow has, nonetheless, a substantial effect on the fluid dynamics as we describe below.
We note that since a flow produced by a reciprocal cycle of the boundaries induces an identity map for the positions of each fluid element at successive cycles, the problem of mixing by nonreciprocal ones is closely related to the class of dynamical systems constituted by perturbations of the identity. A fluid particle that at the beginning of the loop is in a position x ¼ ðx; zÞ, reaches, at the end of the same loop, a unique corresponding point (x 0 , z 0 ) which is a one-toone function (x 0 , z 0 ) = G[(x, z)] of the initial one. For homogeneous fluids, G must also be continuous and differentiable, whereas incompressibility implies that G preserves the area of any domain of points. In other words, incompressible flow implies Hamiltonian dynamics for the fluid particles, and the map that this dynamics induces in one loop is area preserving. For contractible zero-area loops the map is simply the identity; each particle ends in the position in which it started. Hence, a finite-area loop induces, in general, a finite deviation from the identity map and a characteristic value of the geometrical phase gives an estimate for the extent of this deviation. Since generically the geometric phase increases with the area of the loop, for small loops the map is a small perturbation away from the identity whereas loops of greater area induce larger deviations.
Let us now consider the long-term fluid dynamics elicited by a repeated realization of the same contractible non-reciprocal loop that induces a given map. The dynamics is described by the repeated iteration of this map that acts as the stroboscopic map of the time-periodic Hamiltonian system constituted by the incompressible flow periodically driven by the motion of the walls. For small loops, the map is a small perturbation of the identity that can be thought of as the implementation of the Euler algorithm for a putative continuous time dynamical system defined by this perturbation. Therefore, in 2D we expect that the iterations of the map will closely follow the trajectories of this 2D continuous system which is integrable. Therefore, fluid particles will mix very slowly in space: this is, so to speak, mixing by "quasi-static" fluids. This is nicely illustrated in Fig 2(a), where even for a square loop formed with values as large as θ = π/2 the positions of fluid particles after successive loops smoothly shift along the closed curves that are the trajectories of the continuous dynamics. The trajectories are composed of segments that nearly follow the integrable trajectories of a 2D flow (approximated as an Euler map) until it reaches the region of large phase, where chaos and heteroclinic tangles occur. There the particle jumps into another quasi-integrable trajectory, until it again reaches the region of large phase. In typical Hamiltonian chaos (the standard map, for example) the map is not a perturbation of the identity but a perturbation of a linear shear (i.e. with the canonical action-angle dynamical variables (I, φ) following I 0 = I, φ 0 = φ+I 0 ) for which reason this behavior is not normally seen [17,18]. The structure of chaos in this class of dynamics has been greatly overlooked in the literature, and the present research opens a new avenue to the understanding of this associated problem.
As the geometric phase and the corresponding perturbation from the identity map increase, the former argument begins to fail [19]. A more chaotic 2D-area preserving map emerges and with it the corresponding space-filling fully chaotic trajectories. The KAM islands typically become smaller and smaller as the characteristic values of geometric phase increase. As we see in Fig 2(b) for θ = 2π radians, and even more so in Fig 2(c) for θ = 4π radians, after 10000 cycles the fluid particle has covered most of the area available to it between the two cylinders. This is fluid mixing induced entirely by a geometric phase; we may call it geometric mixing. Geometric mixing therefore creates chaotic advection [15], as does the classical journal-bearing protocol.
In Fig 2(d)-2(f) we show the corresponding distributions of the geometric phase over the domain. The value of the geometric phase at a given initial position, obtained in terms of the final angle minus the initial angle in bipolar coordinates (see the Materials and Methods section for further details) after one iteration, F = ξ f − ξ i , is plotted on a color scale of intensities of red (positive) and blue (negative). Note that the phase goes to zero at the walls, as it must, but varies strongly across the domain. In particular, for parameters of θ = 2π radians (Fig 2(e)), we see the development of a tongue of high values of the geometric phase in one sense interpenetrating a region of high values of the phase in the opposite sense. The trajectory plotted in Fig 1 shows the origin of the tongue; fluid particles that are advected to the vicinity of the inner cylinder by the first θ 1 step are then advected to a significantly different value of r by the inner cylinder. As a result, the fluid particle is located on a completely different streamline from the first step when the outer cylinder starts rotating backwards. As may be expected, for smaller parameter values this tongue is absent (Fig 2(d)). At even higher values of θ, on the other hand, (Fig 2  (f)) the tongue wraps twice round in a highly complex fashion. In Fig 2(g)-2(i) we show by plotting the evolution of a line of initial conditions how the geometric phase is related with the dynamical structures in the flow. Fig 2(g), for θ = π/2 radians, shows that when this tongue is absent, the line segment hardly evolves; the flow is almost reversible. The line segments for Fig  2(h) and 2(i), for θ = 2π and 4π radians, on the other hand, show a great deal of stretching induced by this tongue of large geometric phase. To demonstrate this effect of the geometric phase on the flow in more detail, in Fig 3(b) we plot the length of the line segment after a single cycle against the rotation angle. A notable aspect of this plot is that it displays plateaux separated by periods of rapid growth. A comparison with Fig 2(d)-2(f) shows that it is the penetration of the tongue of large values of the geometric phase across this line segment that induces stretching. The tongue penetrates a first time before θ = 2π, and then a second time before θ = 4π, so producing two jumps; between these jumps the evolution of the line segment is much slower. For a given energy cost, which scales with the total unsigned displacement of the walls, geometric mixing is therefore more efficient for a large value of θ.
The journal-bearing flow is just one member of a class of flows that display geometric mixing. In open flows, one has instances such as the well-known Purcell swimmer that can be seen as operating through a geometric phase. Another closed flow that was studied early on in chaotic advection is the rectangular cavity flow, in which one or more of the walls of a fluid filled rectangular container can move, being set up as conveyor belts [15,20]. As in the case of the former studies of the journal-bearing, these mixing protocols imply a cumulative relative displacement of the container walls. However, in the same way as in the journal-bearing case one can introduce a geometric phase by returning all the walls to their initial relative positions after a loop in the parameters. More generally, one can conceive of flows induced by a container in which the walls do not move as rigid bodies, but instead can deform longitudinally and/or tangentially along a nonreciprocal cycle in order to produce efficient mixing. For example, one might consider the case of an elastic bag containing a fluid and subject to the action of a periodic squeezing-distention sequence around one of its sections with a compensating distention-squeezing action around another. This cycle would clearly induce a reciprocating flow unfit for efficient mixing, but again a geometric phase could be made to exist if this spatially stationary configuration were replaced by one that propagates along the bag axis.

Peristaltic mixing
The stomach is a biological instance of such a cavity flow [21,22]. The human stomach is a strong muscular receptacle between the oesophagus and the small intestine. It is not just a storage chamber for food, but also a mixer where the chyme is prepared. The human stomach has a volume L 3 of some 330 mL, while the viscosity μ of the chyme is of order 1 Pa s, its density is ρ % 10 3 kg m −3 , and the maximum flow velocities V observed are in the range 2.5-7.5 mm s −1 [21]. From these data we may estimate the Reynolds number Re = ρVL/μ to lie in the range 0.2-0.5. Thus we may conclude that in the human stomach fluid inertia has only limited importance, and in any smaller animal it will be inappreciable. We note that previous work on gastric mixing have mostly considered the case of inertial contributions [21,23,24] for which the dynamical constraints discussed herein do not apply.
The gastric mixing is brought about by peristaltic waves-transverse traveling waves of contraction-that propagate along the stomach walls at some 2.5 mm s −1 . They are initiated approximately every 20 s, and take some 60 s to pass the length of the stomach, so 2-3 waves are present at one time, while on average the stomach width as the wave passes is 0.6 times its normal width [21,22]. We thus have their velocity c = 2.5 mm s −1 , frequency ω = 0.05 Hz, and thence wavelength λ = c/ω = 5 cm, and their amplitude b = 1/2 × 0.6L % 2 cm. These waves force the stomach through a nonreciprocal loop in the space of shapes, as a result of which geometric mixing is expected. One can give a rough estimate of the size of the expected geometric phase by taking advantage of results obtained for another geometric phase problem: that of low-Reynolds-number microorganisms swimming. Many bacteria swim by deforming their bodies in the same way as the peristaltic waves of the stomach and their speed has been well estimated by modeling such deformations as plane waves [25]. Similar calculations for the stomach render the flow velocity induced by the peristaltic waves V = πc(b/λ) 2 , which comes out at approximately 1 mm s −1 , from where a displacement of about 6 cm per peristaltic cycle is expected or, considering a circular stomach of radius L, a geometric phase of the order of 2 radians.
To show the effects of this phase, we have constructed a minimal model of the stomach undergoing peristalsis, as sketched in Fig 4(a) and further detailed in the Materials and Methods section. We have intentionally reduced the geometric, dynamic, and functional complexity of the stomach and model a 2D section of a tubular stomach of uniform radius, with sealed pyloric and esophageal valves, to focus on the role peristaltic contractions may play in mixing within the enclosed inertialess cavity. Similarly, we have treated the chyme as a Newtonian fluid, leaving the complexities associated with viscoelasticity for future work, as the existence or not of geometric mixing in the stomach is independent of the viscoelastic properties of the chyme. A similar model was used in [4,26] to assess transport by peristaltic pumping in infinite slender tubes. In our model a peristaltic wave deforms the upper and lower boundaries of a symmetric cavity of aspect ratio π according to z w (x, t) = 1+bsin(kx − ωt) in the (x, z)-coordinate system. The flow within the cavity is obtained by integrating the Stokes equations for the velocity field (u, v) with the corresponding boundary conditions for the peristaltic wave at the upper boundary, u = 0 and v = @z w /@t at z = z w (x, t), and symmetry boundary conditions at z = 0. Lateral walls deform vertically to match the vertical velocity of the peristaltic wave at x = 0 and x = 2π. In Fig 4(a) black solid lines represent the streamlines of the induced fluid motion within the cavity due to the peristaltic wave. The contour plot corresponds to the timeaveraged velocity over one full peristaltic cycle. Areas of maximum average velocity are close to the axis of symmetry, whereas near the wall the average velocity is zero and no average motion is produced.
We consider the mixing of a passive scalar wð x; tÞ whose initial spatial distribution at t = 0 is given by the blurred step (wð x; t ¼ 0Þ ¼ 1 þ tanh½ðz=z w À 1=2Þ=dÞ, as represented in the contour map of Fig 4(b). The temporal evolution of this spatial concentration is obtained integrating the advection-diffusion equation for a characteristic Péclet number, Pe = cλ/D chyme representative of the mixing process within the stomach. As the characteristic diffusivity of the Geometric Mixing chyme is, at most, of order of the molecular diffusion of large macromolecules D chyme 10 −6 cm 2 /s, Pe ) 1 and advective contributions dominate the mixing process. Fig 4(c) represents the spatial concentration of the passive scalar χ after 20 peristaltic cycles (i.e. after a rescaled time T = t/T Ã = 20, where T Ã represents the cycle period) for Pe = 15 × 10 3 . The flow induced by peristalsis accumulates a finite geometric phase after each cycle, fluid elements are stretched and folded and, as a consequence, thin filaments are formed that facilitate mixing within the cavity.
We obtain the geometric phase by integrating the trajectory of passive scalars over one full cycle, with uniformly distributed initial conditions in the domain [0, 2π] × [0, z w (x,0)]. The Euclidean distance between the initial and final position after one cycle gives an estimate of the geometric phase. Contours in Fig 4(d) represent the geometric phase of the system. It can be seen that maximum displacements are observed in the central region of the cavity where filaments are created. Note that regions in Fig 4(d) with small displacements correspond to regions that remain unmixed in Fig 4(c). Thus, and despite the uniform radius of the cavity in our minimal geometric model, mixing is not spatially uniform. Regions in the central part of the cavity form thin filaments that enhance mixing, whereas regions close to the lateral and to the upper walls remain almost unmixed after 20 cycles. Even further inhomogeneities are expected for more faithful geometries, with changing average wall diameter [27], specific timing of the opening and closing of the pylorus with peristalsis [28] and interactions between the fundic/cardiac region of the stomach [29], all of which are known feature for mixing within the stomach [30].
Stomach contractions that correspond to a standing wave are akin to a zero-area reciprocal loop. As we anticipated for the journal-bearing case, reciprocal loops induce flow which does not generate any mixing. This is shown in Fig 5(b) where the concentration field after 20 cycles of the boundaries deforming as a standing wave is depicted. Since the induced geometric phase is null, mixing is only controlled by (slow) diffusion.
The importance of the geometric mixing in the stomach may be appreciated by reference to instances in which it is disrupted. The stomach is like the heart, with electrical activity from a pacemaker region stimulating oscillations; in this case being traveling waves of peristalsis. If  Fig 4(c) for the case of a stationary and random wave, respectively. this system is not functioning correctly, there can be gastroparesis or gastric fibrillation [31,32], in which the peristaltic waves become disordered. We have generated such disordered deformations by interspersing peristaltic waves whose propagation velocities c are chosen randomly from a uniform distribution of zero mean. The scalar field χ remains almost unmixed compared to the peristaltic case after an equivalent integration time, with mixing mostly controlled again by slow diffusion (Fig 5(c)). Thus, in our terms, there is poor mixing or no mixing in gastroparesis because there is not a loop around the space of forms, so no average geometric phase, and instead random peristaltic waves induce only mixing and unmixing.
To compare the degree of mixing in the three cases considered herein (peristalsis (pw), stationary (sw) and random (rw) waves), we calculate for each cycle the variance of the spatial concentration field [33,34], σ = h(χ − hχi) 2 i 1/2 , where hi denotes the spatial average. Fig 5(a) represents the evolution of σ with the number of cycles. It reveals the higher mixing efficiency realized in peristalsis by geometric mixing.

Discussion
In summary, we have introduced the concept of geometric mixing in which mixing arises as a consequence of a geometric phase induced by a contractible non-reciprocal cycle in the parameters defining the shape of the container. It turns out that the mixing efficiency estimated from the stretching of material lines is roughly proportional to the geometric phase. Mixing in the corresponding flows can be also considered as the result of chaos arising in the mapping describing the motion of fluid elements during one cycle. When the cycle is reciprocal, this map is the identity and a small departure from reciprocity corresponds to a small departure from the identity map. The chaotic properties of maps neighboring the identity have been poorly studied in the past. They also arise in the quite different context of numerical integration methods of ordinary differential equations in the limit where the step size tends to zero [19]. Our results are therefore also relevant to the characterization of chaos in this class of systems. Lastly, we have shown that such a geometric phase-the "belly phase" [35]-may be found in the stomachs of animals where Re < 1.

Supporting Information Journal bearing flow
The journal-bearing flow has been widely employed to study the process of mixing in laminar flows. Fig 6 shows a sketch of the configuration studied herein. The outer cylinder of radius R out rotates with an angular velocity O out , whereas the inner cylinder of radius R in rotates with an angular velocity O in . The eccentricity of the inner cylinder is given by ε. In the limit where viscous forces are negligible, the resulting flow is obtained by integrating the biharmonic equation for the stream function r 4 ψ = 0 with the corresponding boundary conditions at the walls of the cylinders.
Because of the linearity of the problem the solution for the stream function can be written as where ψ out is the solution for the stream function of the flow induced by the outer cylinder, whereas ψ in corresponds to the solution of the stream function of the flow induced by the inner cylinder. The angle covered by a cylinder during one cycle depends on its angular velocity according to where the subindex i denotes the outer or inner cylinder and T Ã represents the period of the cycle. Since, in the simulations considered in this paper the angular velocity of the cylinders is constant, Θ i = T Ã O i . This flow admits an exact solution [11] for the stream function when the problem is written in bipolar coordinates, (ξ, η). The cartesian coordinates (x, z) can be recovered according to Following [11] the solution for the inner and outer stream functions is given by where H = b/(c 2 +s 2 ) 1/2 , with s = sin ξ sin η and c = cosξ cos η − 1. Moreover, and with where ξ in and ξ out represent the surfaces of the inner and outer cylinders, respectively, and Δ, D, Δ Ã , h 1 , h 2 , . . ., h 8 are given by Once the flow is evaluated, the trajectories of the particles were obtained integrating The integration of Eqs (10) and (11) was carried out with a fourth order Runge-Kutta scheme. After one complete cycle, both cylinders end at their initial position, while particles departing from an initial position (ξ i , η i ) are located at (ξ f , η f ). A stroboscopic map was constructed from the iterative application of the described loops which characterizes qualitatively the dynamics of particles in the journal-bearing flow. Final positions of the particles after each full iteration were plotted to distinguish between particles that follow smooth integrable trajectories from particles following space-filling fully chaotic trajectories (see Fig 2 of the main text). We compute the geometric phase as the travelled distance measure as the difference in the bipolar angle F = ξ f −ξ i . We also evaluate the mixing efficiency by computing the length by which an initial fluid segment is stretched as a function of the number of mixing cycles. For this, we densely sample the initial fluid segment with N particles, r i (0), i = 1, . . ., N, and track the motion of each one by the method described above after T iterations of the mixing protocol, obtaining the final position of each particle, r i (n = T), i = 1, . . ., N. We then compute the total length of the final segment by adding the distances between adjacent points as P j r i ðTÞ À r iÀ1 ðTÞ j N i¼2 . The true segment length is given by the limit of this expression for N ! 1. In practice, we make sure the segment was sampled densely enough by computing the final segment length using N and 2N particles, with N sufficiently large such that both computations give the same result to within a fraction of the initial segment length.

Cavity flow
We studied the flow within a deformable symmetric cavity, sketched in Fig 7. A peristaltic wave propagates along the upper and lower walls of a two dimensional cavity of length 2π, according to z w (x, t) = 1 + b sin(kx − ωt), where b is the amplitude of the wave, k its wavenumber and ω its angular velocity. This two-dimensional configuration corresponds to a section of a cylindrically symmetric cavity around the x-axis (i.e. a symmetrically deforming tube). All lengths are normalized by the undeformed tube radius and we keep the name of the non- dimensional variables the same for simplicity. We set the cavity length-to-width aspect ratio to π in all simulations.
Numerical solution to the cavity flow. The propagation of the waves along the no-slip boundaries induces a flow within the cavity that, in the limit where viscous forces are dominant, is obtained by integrating the two-dimensional Stokes equations @u @x þ @v @z ¼ 0; ð12Þ with the corresponding boundary conditions: u = 0 and v = @z w /@t at the upper wall (i.e. z = z w ) and @u/@z = 0 and v = 0 at z = 0, to impose symmetry conditions. Lateral walls deform vertically (with a linear vertical displacement distribution) to match the vertical velocity of the peristaltic wave: u = 0 and v = z@z w /@tj x = 0 at x = 0, and u = 0 and v = z@z w /@tj x = 2π at x = 2π. Initial conditions correspond to the fluid at rest. To facilitate the numerical integration of the problem we employed the vorticity-stream function formulation [36]. Thus, the problem reduces to that of integrating where o ¼ Àð@u=@z À @v=@xÞ is the vorticity of the flow and ψ its stream function. Accordingly, the boundary conditions given above for the velocity field have to be expressed in terms of o and ψ.
We note that the computational domain changes as the moving boundaries evolve in time. However, this domain was transformed into a fixed domain by means of a mesh transformation [37] with a new computational vertical coordinate, Z = z/z w (x, t). Eqs (15) and (16) were rewritten in terms of the new variable Z and were discretized with a second-order accurate finite difference scheme. The coupled system of Eqs (15) and (16) was then solved with a lineby-line Thomas algorithm [37] at each time step. The iterative process to solve Eqs (15) and (16) starts with an initial guess for the vorticity and stream function fields which corresponded with the values of the previous time step. This system of equations was solved iteratively until convergence was achieved. The flow in a squared cavity studied in [38] was used as a check for our numerical scheme.
We integrate numerically Eqs (15) and (16) for a fixed set of parameters with non-dimensional values consistent with the available experimental data for a human stomach. We consider peristaltic waves with b = 0.4, k = 3 and ω = 3. The flow created by a standing wave was obtained by superposing two peristaltic waves with opposite velocities: the first moving from left to right with ω = 3 whereas the second propagates from right to left with ω = −3. As the resulting flows are time periodic in these two cases, the integration only needs to be carried out for a single wave period, t = T Ã = 2π/ω, saving computational time.
To generate a random wave, we choose the angular velocities, ω i , sampling them randomly at each period T from a continuous uniform distribution with zero mean in the interval [−3,3], while keeping both the wavenumber and the amplitude fixed to the same values used above.
Thus, values of ω can be either positive or negative, corresponding to waves propagating respectively from left to right or viceversa with random angular velocity. Since the resulting flow is not periodic in this case, the numerical integration was carried out for a prescribed integration time equivalent to 100 periods of the peristaltic wave, which facilitates the comparison among the different scenarios. To ensure a continuous evolution of the waveform at the end of each period (i.e. coincident with a change in ω), we shift the new waveform along the x-axis (i.e. we introduce a phase shift) to match the initial condition for the succeeding period.
Mixing within the cavity. To analyze mixing within the cavity we integrated the twodimensional advection-diffusion equation for a passive scalar wð x; tÞ: @w @t þ u @w @x þ v @w @z À 1 Pe that is initially distributed in a blurred step (wð x; t ¼ 0Þ ¼ 1 þ tanh½ðz=z w À 1=2Þ=dÞ), with δ = 0.1 its characteristic thickness (see Fig 4(b) of the main text). The advection-diffusion process is characterized by a non-dimensional Péclet number, Pe = cλ/D chyme , which in our simulations was fixed to Pe = 15 × 10 3 since in the case considered here advection largely dominates over diffusion. Eq (17) was integrated with no-flux boundary conditions uw À 1=Perw ¼ 0 at the boundaries of the computational domain, together with the symmetry boundary condition at z = 0. Since the velocity field was provided by the integration of Eqs (15)-(16), we used the same transformation of the domain and the same spatial discretization described above for the numerical integration of Eq (17) together with a Crank-Nicholson scheme for the temporal evolution.
To compare the degree of mixing in the three cases considered herein (peristalsis (pw), stationary (sw) and random (rw) waves), we calculate for each cycle the variance of the spatial concentration field [33,34], σ = h(χ − hχi) 2 i 1/2 , where hi denotes the spatial average. Fig 5(a) of the main text shows the evolution of σ with the number of cycles for the three cases we have considered.
Finally, to determine the geometric phase induced by the peristaltic wave we calculated the trajectories of uniformly distributed passive particles within the cavity for one compete period of the peristaltic wave. The trajectories were obtained integrating dx/dt = u, dz/dt = v with a fourth-order Runge-Kutta scheme. As (u, v) were evaluated at each time step on the mesh employed to integrate Eqs (15) and (16), we employed a 2D linear interpolation to obtain their values at the particles position. The distance from the position of the particle at the end of one period to its initial position was used to determine the spatial distribution of the geometric phase, F, represented in Fig 4(d) of the main text.