Interfacial dynamics in 3D binary fluid demixing: animation studies*

The late-stage phase ordering, in three dimensions, of fully symmetric binary fluid mixtures is studied via a lattice Boltzmann method. We present time-resolved maps of the fluid velocity fields and also animated visualizations of the interfacial motion. These show distinct features corresponding to regimes where viscous, crossover and inertial hydrodynamic scaling have previously been identified. Specifically, while the interface is overdamped in the viscous regime, it exhibits recoil after topological reconnection at intermediate and higher inertia; and in our most inertial runs the interface shows extensive underdamped capillary disturbances not attributable to topological reconnection events. The advantages and practicality of presenting such dynamical data in fully animated form are demonstrated and briefly discussed. This paper's animations are available from the Multimedia Enhancements page as individual files and also packed into archives (two formats).


Introduction
Spinodal decomposition occurs when a fluid mixture of two species A and B, forming a single homogeneous phase at high temperature, undergoes spontaneous demixing following a sudden temperature quench [1]. For compositions close to 50/50, one enters the 'spinodal' regime in which the initial homogeneous phase is locally unstable to small fluctuations. There then arises, after an early period of interdiffusion, a bicontinuous domain structure in which patches of A-rich and B-rich fluid are separated by sharply defined interfaces. In this late-stage structure, the local compositions of the fluid patches correspond to those of the two bulk phases in coexistence; the interfacial tension approaches σ, its equilibrium value. Local interfacial curvature causes stresses (equivalently, Laplace pressures) to arise, which drive fluid motion; for a deep quench, diffusion ceases to contribute directly to coarsening. The interface then evolves smoothly with time between isolated 'pinch-off events' or topological reconnections. For a bicontinuous structure in 3D, pinchoff events provide the coarsening mechanism whereas the reverse process of coalescence is rarely seen (this is a different situation from 2D, where coalescence and pinchoff are the same mechanism). It is usually assumed that, once started, these pinch-off events occur rapidly enough not to impede the coarsening process; and that at late times, the presence of thermal noise in the system is irrelevant. The problem is thus one of deterministic, isothermal fluid motion coupled to a moving interface. (For a critique of some of these assumptions see [2,3] and for a detailed analysis of two-fluid pinch-off, see [4]. ) The dynamics are then described by the Navier-Stokes and the advective diffusion equations [1]: where u is the fluid velocity, D t the material derivative, P the pressure, µ the chemical potential for the order parameter φ, η the viscosity, ρ the density and M an Onsager coefficient, related to the diffusivity. Here µ = δF/δφ where we choose F [φ] = (A/2)φ 2 + (B/4)φ 4 + (κ/2)(∇φ) 2 as a suitable free energy functional, selected for its simplicity. A and B are phenomenological coefficients that depend implicitly on temperature. At high temperatures A ≥ 0 and the two fluids are miscible. Demixing takes place when A ≤ 0. We consider, without loss of generality for our purposes, the particular case |A| = B, for which the equilibrium order parameter of the two equilibrium phases at low temperatures corresponds to φ = +1 and φ = −1. During the phase separation process, once the interfaces between the two phase domains are sharp, the driving term (∼ φ∇µ) in equation (1) is just the Laplace pressure of interfacial curvature [1], with interfacial tension σ = (−8κA 3 /9B 2 ) 1/2 . In this paper, these equations are solved by a lattice Boltzmann method outlined briefly in section 2 below. For simplicity, we address only the fully symmetric case of two incompressible fluids with identical physical properties (viscosity η, density ρ), and also equal volume fractions. The dynamical scaling hypothesis is that, if one defines reduced physical units of length and of time by then at late times any given characteristic structural length l(t) should evolve with time t according to where ϕ(x) is the same function for all such fluids. We generally adopt these units from now on; when we wish to refer to lengths and times measured in simulation units (i.e., the lattice constant and time step) we instead denote these by L, T .
There is now good numerical evidence for dynamical scaling in three dimensions, though not in two [5,6]; see [2] for much of the relevant data, and a survey of earlier work. Allowing for an early time diffusive offset (related to domain formation, and absorbed for each simulation run into a shift of the time origin), one finds a universal scaling curve l(t) as implied by equation (4). The curve (see figure 1) shows two limiting time regimes, separated by a broad crossover. Choosing as the length measure l ≡ 2π S(k)dk/ kS(k)dk with S(k) the density autocorrelator, the early and late time asymptotes are quantified as l = 0.07t and l = 0.9t 2/3 ; these cross at t * 10 4 , a rather large value (though formally of order unity) [2]. Note that when early diffusion and late finite size effects are carefully eliminated, any individual simulation run (these data are for 256 3 lattices) covers only about one decade of the l(t) plot. Accordingly, each run is located either wholly in one asymptotic regime or the other, or is wholly in the crossover region.
We interpret the two observed asymptotic regimes as corresponding to the viscous hydrodynamic (VH, l ∼ t) and inertial hydrodynamic (IH, l ∼ t 2/3 ) regimes as proposed for early and late times respectively by Siggia [7] and Furukawa [8]. In the former, the primary balance of terms in equations (1)-(2) is between the interfacial stress (driving) and the viscous dissipation; in the latter, it is between the driving term and acceleration. Note, however, that since the interfacial area and the kinetic energy are both decreasing, the latter balance cannot be maintained in the energy conservation equation derivable from equations (1)-(2) [9]. This conservation demands the presence of an energy cascade in which the reduction in interfacial area at large scales feeds 9.4 kinetic energy down to smaller scales where viscous dissipation finally removes energy from the system. A three-length scale analysis of this situation predicts, nonetheless, the same primary balance in the Navier-Stokes equation and the same scaling law [9].
Beyond the observed scaling asymptotes as shown in figure 1, evidence for the identification of VH and IH regimes is presented in [2] based on a number of interfacial and velocity statistics, many of them in Fourier space. (For example, there is a distinctive shoulder on the structure factor S(k) which is less marked in the inertial regime than in the viscous one.) The importance of inertia is usually measured by a Reynolds number. Such a number is usually intended to estimate the ratio of the nonlinear to the viscous term in the Navier-Stokes equation. In fluid demixing, the order-parameter Reynolds number Re φ = ll (so called because it is derivable from the l(t) curve without evaluating any velocity statistics) is normally chosen. This is of order 10 −1 for the most viscous (leftmost) runs in figure 1, of order 10 in the crossover region, and approaches 350 late in the most inertial runs shown. Despite the objections of Grant and Elder [10], it is possible in principle for Re φ to diverge with time in the IH regime; Kendon [9] has argued that, once a turbulent energy cascade develops, Re φ ceases to be an accurate measure of the ratio of convective and viscous terms in Navier-Stokes. (The latter remains finite according to [9].) It should be borne in mind, anyway, that a fully asymptotic IH regime is not yet accessible by simulation (nor likely to be soon); however, we are able to reach a regime where the inertial terms dominate over viscous ones in Navier-Stokes by at least an order of magnitude [2]. The relatively slow increase of Reynolds number (however defined) upon moving along the l(t) curve is consistent with the very wide crossover connecting VH and IH regimes.
In the present work, we present visualizations that not only confirm the identifications of the two regimes, but show with some clarity the difference in dynamics, particularly of the fluid-fluid interface, between the VH and IH regimes. Such differences are not unexpected from previous theoretical discussions [7,8]; however, their direct confirmation has had to await the simulations presented here. Despite the lack of quantitative analysis, these animations provide convincing and complementary insight to that gained through aggregate statistics of (for example) velocity autocorrelations [2]: in animations the distinct physical mechanisms for interfacial dynamics are directly observable.

Simulation methodology
All simulations were performed using Ludwig, an optimized, modular lattice Boltzmann code designed for simulating complex fluid physics on large parallel machines [11]. The binary fluid lattice Boltzmann methodology (first used in [12]) employed by Ludwig is exhaustively reported in [2], where a careful analysis of random and systematic errors is also made. In this method one simulates a linearized Boltzmann equation on a lattice. The basic dynamic quantity is the one-body distribution function; its moments are related to the appropriate conserved variables of the system. For a binary mixture, two distribution functions are introduced; one related to the density, ρ, and momentum, ρu, of the fluid, the second to the order parameter (concentration), φ. Their relaxation parameters are related to the transport coefficients of the fluid: to the viscosity for the former, to the order parameter diffusivity for the latter. We have considered parameters such that, in lattice units, ρ = 1 and φ varies between −1 and 1 (see [2,11] for further details). It has been shown that in the hydrodynamic limit the appropriate macroscopic equations (1)- (2) are recovered from the corresponding linearized lattice Boltzmann equations [2]. The advantage of keeping the diffusive order parameter field in the simulations is that topological reconnections Table 1. Lattice Boltzmann input parameters: free energy parameters A, B, κ and diffusivity M (lattice units). Fluid parameters: viscosity η, surface tension γ, characteristic length l 0 and time scales t 0 as given in equation (3) (lattice units). The density is unity and is maintained almost constant by using an ideal gas equation of state for the overall density, ρ, at very low Mach number [2]. of the interfaces are handled without encountering a singularity. In lattice Boltzmann local diffusion is the process that deals with these singularities. However, since we are interested in the limiting case when domain growth is limited by hydrodynamics, we must carefully choose the order parameter diffusivity. We show in table 1 the parameter values used in the simulation runs for which data is presented below. (Mostly these are the same parameters as used in [2] †.) A, B, and κ are chosen to to maintain a constant interfacial thickness while tension is varied. In particular, the interfacial thickness is large enough to avoid anisotropic contributions to the interfacial properties arising from the underlying lattice [2]. Finally, the order parameter mobility M is chosen to ensure that interfaces remain locally in equilibrium while the direct diffusive contribution to the domain coarsening rate is less than 2% of the total. This (necessary) compromise does not ensure that diffusion is everywhere negligible: in our animations, we frequently see how small satellite droplets, left over after a topological reconnection, disappear slowly by evaporation rather than by coalescence. (It is expected that diffusion becomes relevant when the local domain size is not large compared with the interfacial thickness.) However, we are sure that hydrodynamic flow dominates the overall coarsening dynamics, except for run DIFF, discussed separately below (figure 3).

Animation issues
The new results of this article are presented as animated graphics (section 4). The advent of epublishing allows such data to be presented within a refereed article, not merely as supplementary material, but as the main focus of discussion. We believe that real scientific insight can be conveyed by this route. Simulators, but also experimentalists, may wish to take note of our experience in preparing the animations for publication.
In contrast to the 256 3 runs used to generate figure 1, it was found that 128 3 runs was the largest for which animated datasets could be produced for Web use with reasonable resolution and download time. This is, coincidentally, about the largest system size that can be watched without visual information overload for a typical viewer; most people prefer to watch a 64 3 animation 9.6 in which the domain scale L varies from about a tenth to about a half of the linear lattice size (Λ). Since best practice [2,3] for avoidance of finite size effects is to maintain L ≤ Λ/4, such animations were created by cropping the 128 3 data.
A file size limit of about 6 Mb per animation is appropriate to maintain reasonable download characteristics. For an animation of, say, 30 seconds viewing time, this means compromising significantly on resolution. However, the results do allow appreciation of the global pattern of coarsening behaviour. A higher resolution has been adopted in the shorter animations (showing a smaller portion of the interfacial dynamics) used below to draw attention to particular physical features.
Note that the advantages of higher spatial resolution are significantly reduced if the animation then exceeds local limits and cannot run at the intended frame rate. For our coarsening data, the minimum frame rate to deceive the eye into perceiving smooth motion is about 10/s. This deception is important as motion and static features are detected by distinct parts of the human visual system [13]. Thus in animation, a dynamical mechanism can sometimes be immediately apprehended, when scrutiny of a series of still frames would never convince. (Conversely in animation the attention is drawn to parts of the image where motion is most rapid. This can allow relatively static features to go unnoticed.) Most of the discussion will be focused on 3D animations showing the interfacial motion only, where it is easier to identify and analyse the relevant dynamical features in each regime. Animated 3D velocity maps are also presented below, but are less easily assimilated, and on balance slightly less informative, for two reasons. Firstly, to avoid visual overload one can only present velocity map data for a thin slab (though the displayed velocities are fully 3D). Secondly, to make efficient use of our parallel LB algorithm, velocity data had to be output with a fairly long stride (typically every 100 steps). This leads to a slow frame rate and a large image increment between frames so that the deception of continuous motion cannot be maintained.
Animations were produced from LB datasets using the AVS package. This was configured to create order parameter isosurfaces at φ = ±0.15 coloured yellow and blue. The effect is as if the two sides of an ideal fluid-fluid interface, defined by the φ = 0 countour, were painted in the two colours (a useful aid to visualization). In addition, the finite thickness of the region bound by these surfaces gives an indication of the interfacial width. Depth perception can be improved by having two white light sources-an ambient light and a direct (bi-directional) light source.
In the velocity animations, 3D velocity maps are superimposed on the corresponding slab of the previous isosurfaces. We have coloured spectrally small pieces of 3D streamlines according to the true (not projected) 3D velocity magnitude (blue = slowest, red = fastest). All the streamlines have a fixed length; thus a short-looking red segment has a high velocity but is oriented towards or away from the viewer. Several alternative ways of imaging velocities were tried, but this one was the most successful.
The AVS images were rendered into AVI animations using DPS Video Action PVR v6.33 on a DELL Precision 610 machine. AVS image creation and their rendering into a 30 second, 9 Mb animation takes about 100 processor hours for the simulation and post-processing time, and approximately one day of trained staff time for data exploration and frame generation. AVI has superior resolution and compression characteristics to either animated GIF or MPEG; neither of these is adequate for our purposes. While AVI viewers are supplied as standard on most platforms including Unix, Web browsers are often not configured to use them by default †. † For help in accessing a viewer, see http://www.epcc.ed.ac.uk/desplat/NJP/AVI help.html. 9.7

Results
To analyse the phase separation dynamics of an immiscible binary mixture, we have considered 50/50 mixtures, starting from a homogenous φ = 0 initial state with uniform white noise added. For these initial conditions the spinodal textures consist of a bicontinuous structure, where hydrodynamics plays a relevant role. The different parameters of the simulations carried out are displayed in table 1.
Within each simulation run, there is an early-time diffusive regime which is excluded from the animations. Then, depending on where on the universal curve ( figure 1) the run lies, the main evolution will display either VH, IH or crossover regime behaviour. (The breadth of the l(t) plot ensures that, after the diffusive stage, one does not see a significant change of behaviour within a given run, i.e. only one dynamic regime can be sampled at a time.) The mechanistic differences between regimes are most easily seen when time is scaled so that the perceived overall coarsening rate is comparable between runs. Use of raw simulation units nearly achieves this (since the overall coarsening rate also sets the optimal discretization). In the animations we display the time corresponding to each frame in lattice units, T . The corresponding dimensionless time can be obtained from the data of table 1 using the fact that t = T/t 0 .
In the animations we are showing the evolution of the domain interfaces within a box of predetermined size. Hence, the interfaces are cut at the surface of the limiting box, and moving interfaces will enter and leave the box. Due to this effect, looking at the surface of the box rather than the bulk may give a misleading impression that coarsening in three dimensions proceeds by coalescence rather than pinch-off.

VH regime
This regime is the one accessed by most previous simulations, and (for reasons discussed in [2]) all experimental studies that we are aware of so far. It corresponds to reduced physical times and lengths such that t ≤ 3 × 10 2 , l ≤ 10. In this regime we have performed simulations with parameters corresponding to run VH1 of table 1, and show animations of the dynamics of the evolution of the interface in figure 2.
As expected from the fact that viscous and interfacial forces are in balance, with negligible inertia, the interface remains smooth over most of its surface (except close to pinch-off) and its motion is overdamped. To the untrained eye the interface motion superficially resembles diffusive bicontinuous domain growth (as in an A/B solid alloy, or an infinite viscosity fluid, for example; see figure 3 below), but close inspection shows clear collective motions of the interface at a scale larger than the domain size. We interpret this feature to be the result of long-range spatial correlations in the velocity field, which are expected in the VH regime [2]. The animation also reveals interface-driven squeezing of fluid out of a narrowing neck region as one approaches the topological reconnection that corresponds to breakage of the fluid neck. The collective interface motion may be related to the long range spatial correlations that characterize the fluid flow in this regime, for which a tentative explanation was given in [2] (see also section 4.3 below).
Domain pinch-off proceeds via a cylindrical tube which thins fairly symmetrically until reconnection occurs either at the midpoint (usually) or nearer one end (sometimes). Figure 2(c) shows a good example of the latter; symmetry is maintained until just before the break. The actual reconnection is achieved by rapid diffusion which sets in only when the cylinder's width becomes of order of the interfacial width. During this step there is a slight diffusive loss of  material from the cylindrical region, and any satellite fragments of the local minority phase disappear by evaporation rather than coalescence. While essential to allow reconnection in LB, the diffusivity is maintained at a low enough level that it contributes negligibly to the overall coarsening rate. In fact, diffusivity is a true physical mechanism, but it is somewhat amplified, for numerical efficiency reasons, in our simulations [2].
The post-reconnection dynamics is characterized by the retraction of the two fingers created. Initially, it may be dominated by tip evaporation, but subsequently it shows an overdamped relaxation, with each withdrawing finger collapsing smoothly back towards the surrounding interface. As expected in the absence of inertia, there is no appreciable recoil (overshoot) of the withdrawing finger, although often the surrounding interface is already moving in the direction of withdrawal and continues to do so afterwards. However the distinction between this and true recoil becomes readily apparent when inertial runs are viewed (section 4.1.3).
The presence of interface-correlated motion induced by the flow can be better appreciated by comparing the present animations with those of a fluid that combines extremely high viscosity with a relatively high mobility M . In this case diffusion is the main mechanism for domain coarsening, and the structural length grows as l ∼ t 1/3 [1]. In figure 3 we show an animation of the spinodal decomposition for a binary mixture with the parameters given for run DIFF in table 1. These are the same as for VH1 (figure 2) except with a much larger viscosity, η = 1000. This suppresses hydrodynamic coarsening and delays the early-time crossover from diffusive to viscous growth beyond the domain sizes analysed here. (We have verified that the l ∼ t 1/3 law holds for diffusive growth throughout run DIFF.) The use of such a large viscosity may induce spurious effects in the fluid dynamics at the lattice scale [2], but it does not affect the large scale interfacial motion we are interested in here. The interfaces for run DIFF look smooth and regular, as they do for run VH1. However, a careful comparison between figures 2 and 3 shows that in the diffusive case the correlated motion between interfaces, at a scale larger than the domains themselves, is lost.  t ≤ 10 6 , which corresponds to dimensionless domain sizes between 10 ≤ l ≤ 10 4 . Runs in this regime have fitted scaling exponents discernibly different from 1 or 2/3 [2]. We have analysed two simulations in this regime, one in the early crossover, the other late. As expected from the fact that this is a slow crossover not a separate scaling regime, these two runs differ significantly.

Crossover behaviour
In figure 4 we show animations for run CROSS1 in table 1, which lies in the early section of the crossover region, with t/t * close to but less than unity (see figure 1). A slight degree of interface recoil after pinch-off is detectable (if only on close inspection) after some events and takes the form of an overshoot after finger withdrawal to create a local bump of the opposite curvature to the initial finger. However, any bump spreads rapidly and is lost against the background motion of the interface, so the overall visual appearance is not very different from the VH regime; nonetheless, in some cases one sees traces of propagating capillary waves with interfacial bumps moving a distance at least of order the domain size before decaying.
All of these features are the direct result of fluid inertia; one expects them to become more pronounced as one moves further to the right on the l(t) plot. Indeed, each can be seen far more convincingly in figure 5 which includes a short high resolution (cropped) sequence from run CROSS2 of table 1. This run lies in the late crossover region, corresponding to dimensionless times 10 5 ≤ t ≤ 10 6 . This part of the l(t) plot is nearing the l ∼ t 2/3 asymptote and some would classify this as an early IH run, although the behaviour continues to change as the inertia is further increased (see below). A fairly noticeable feature, as one moves up the l(t) curve, is that a higher proportion of pinch-off events are strongly asymmetric.
(Several examples of this can be seen in figure 5.) The initial neck sometimes starts cylindrical, but an asymmetric profile, narrower at one end, is very often seen before the final disconnection. Although this may be partly a coincidence, theoretical as well as experimental observations do suggest that pinch-off ought to become strongly and universally asymmetric for low viscosity fluids, where the length l 0 is small on the length scale of observation, as certainly applies in our IH simulations [4,14]. While our pinch-off shapes do not approach any universal asymmetric profile (e.g. a narrow 'cone' contacting a nearly flat 'plate' [4,14]), the simulations possibly suggest the idea that such profiles may predominate as the 9.12 relative importance of fluid inertia is increased. Nonetheless, we cannot rule out the possibility that the asymmetric pinch-off events we see are induced by the random action of interfacial capillary waves. figure 6, we show data from run IH1 in table 1. This run is solidly in the IH regime as far as figure 1 is concerned (l ≥ 10 4 , t ≥ 10 6 ), and indeed lies close to the limit of what can be achieved in LB before encountering resolution problems with the smallest scale velocity features (the Kolmogorov dissipation length is very close to the lattice constant) [2]. The behaviour is totally dissimilar to the VH regime and quite strongly different even from the late crossover run (figure 5), although inertia is relevant in both cases.

Inertial regime In
The previously mentioned features remain clearly visible: dominantly asymmetric pinchoff events; strong recoil after pinch-off; capillary spreading of the recoil bump. On top of this, however, one detects a constant throbbing of the interface with obviously under-damped capillary wave motion. This is particularly clear at high resolution in figure 6(c). However, the onset of such behaviour is gradual: some features resulting from travelling capillary waves were already visible in figure 5(b) for the late crossover region (run CROSS2). Recently, the possible relevance of capillary wave damping as a determining factor in the domain growth rate has been pointed out and argued to lead to an asymptotic l ∼ t 4/7 behaviour at late times [15]. Our animations show that well developed capillary waves are not incompatible with the 'standard' IH domain growth law, l ∼ t 2/3 [8].
One thing we have not detected, even in our most inertial run (IH1), is the occurrence of a reverse pinch where two bumps on opposing sections of interface collide to form a new neck. In three dimensions (but not in two), this process is topologically distinct from the normal pinch-off and actually opposes coarsening by increasing the number of necks per unit volume. (The latter is related to the Euler number density of the interface.) However, if the relative amplitude of the underdamped capillary motion were to increase much beyond what we see, such events could be expected. Nonetheless, so far we have no evidence that the amplitude of the capillary waves increases faster than the domain size, as would be needed for this to happen.
The fact that the interfacial motion is distinctly different between runs CROSS2 (late crossover) and IH1 (the highest inertia we can reach) may be significant. Conventional scaling theories [8] assert that, within the IH regime, temporal and spatial correlation functions should approach a universal form. While this collapse does occur for the fluid-fluid correlator, the same does not occur for velocity autocorrelators and related quantities [2]. This has led Kendon [9] to propose a more complicated scaling analysis in which several length scales are involved. The apparently continuing evolution of the interfacial dynamics between CROSS2 and IH1 is not inconsistent with the idea that IH is not a simple scaling regime as originally proposed [8], although further runs would be needed to confirm this. The possibility of an eventual crossover to a scaling exponent for l(t) other than 2/3 remains an open question [10,15].

Static interfacial structure
As mentioned in Section 2, animations may distract viewers' attention from slowly evolving features. Therefore, in this section we show a few still frames from representative cropped datasets discussed earlier (Section 4.1).

9.15
In figure 7 one can see how, throughout the accessible part of the l(t) curve the interface remains locally quite smooth, but there are clearly more bumps on it in the IH regime; these bumps are the static signature of the capillary waves.
The regularity and smoothness of the interface in the VH regime can be understood as a result of the overdamped fluid motion; the interface reaches locally its equilibrium shape more easily in the absence of propagating perturbations. Another feature, perhaps easier to see in still frame than in animation, is a relative enhancement of the number of thin necks per unit volume as one moves up the l(t) curve. The latter is particularly visible in figure 7(b) although this is a somewhat extreme example and the quantity does fluctuate strongly in time (cf 7(c) for example). For some viewers, the gross discrepancy of static structure between figure 5 (early inertial, figure 7(b)) and figure 2(b) (viscous, figure 7(a)) is overlooked while viewing the animations but obvious in these still frames.
The experienced eye can often use the two aforementioned features (bumps and necks) in combination to locate roughly where on the l(t) plot a particular still frame comes from. This is consistent with the significant evolution of the static structure factor S(q) on moving between VH and IH regimes [2,3], where a pronounced shoulder at wavenumbers beyond the peak in S(q), visible in VH, fades away in IH. This could be linked to an increase in 'structural polydispersity' (more thin necks and/or bumpy surface features).

Velocity maps
In this section we present animated velocity maps extracted from the same datasets as presented previously (table 1).
The velocity vectors are colour indexed in simulation units which again ensures approximately (within a factor of two) that the mid-run (rms) velocities coincide across the various runs. Figure 8(a) for the VH regime shows, as expected, that the flow is fairly laminar with no sign of convection or other fluid motion not driven directly by the interface; the motion is consistent with the existence of long range spatial correlations in VH [2]. Figure 8(b) (early crossover) shows a somewhat more heterogeneous velocity field with 'hot spots' of relatively high velocity appearing from time to time. However, there is still not much visible effect of inertia on the velocity map; any incipient vorticial structure can be plausibly attributed to the swirling effects of the interface itself, although a full analysis of the velocity statistics does show a difference between the crossover regime and VH [2]. In figure 8(c) for the IH regime, the mean coarsening rate is somewhat slower than in 8(b), so the early occurrence of 'hot spots' is more significant. Also, the velocity field has, throughout the run, non-laminar structure at a scale smaller than the mean domain size. This is very visible early in the animation. Later on, there are strong heterogeneities in local fluid velocity on the domain scale and even above; these are perhaps suggestive of spatial intermittency as might be expected once a turbulent velocity cascade starts to develop. For statistical evidence of a related 'turbulence in patches' scenario, see [2].
These velocity maps provide further confirmation of the expected role of inertia in the crossover and IH regimes. However, to many viewers the interface-only animations are more convincing and show more mechanistic insight (the imagination can readily supply the whirling velocity field that the interface motion in, e.g., figure 6(b) must imply). It would be possible to take velocity animations further, to study for example animated vorticity maps-see [2] for a still frame-but this does not seem warranted at present.

Conclusions
We have presented animated datasets showing interface dynamics and velocity maps for the late stage coarsening of symmetric binary fluids. These data show clearly the increasing role of fluid inertia as one progresses along the coarsening curve l(t), and allowed us to identify a number of mechanistic features which could not be convincingly deduced from either still frames or from average statistics for structural and velocity correlations [2]. In particular, on passing through the crossover and into the inertial hydrodynamic regime we saw (i) an increasing tendency to asymmetric pinch-off during topological reconnection of fluid domains; (ii) an increasing degree of finger recoil following pinch-off, with overshoot to create a dimple of the opposite sense to the original finger; (iii) an increasing tendency to florid underdamped capillary motion on the fluid-fluid interface. The latter tendency is particularly marked in our most inertial run (IH1), where the interface is constantly alive with propagating capillary waves even far (in space and time) from any topological reconnection event. Fluid velocity maps show consistent features, with a progression from laminar to more unruly flow as one enters the inertial hydrodynamic regime; these data concur with statistical measures, reported elsewhere [2], suggesting the onset of partial fluid turbulence. In the VH regime, the structure of the (locally laminar) fluid flow appears to exhibit long-range correlations [2]; additional analysis and visualization work might expose more clearly the spatial patterns that are present (e.g. by looking at the evolution of vorticity isosurfaces). The ultimate nature of the coarsening process, beyond the range of l, t we can access, remains an open issue [2,9,10]. Experimentally this may prove academic: as discussed in [2,9], because of the intervention of gravity, it may prove extremely difficult for laboratory experiments to get further into the inertial regime than our simulations do.