Coherent Particle Structures in High-Prandtl-Number Liquid Bridges

Clustering of small rigid spherical particles into particle accumulation structures (PAS) is studied numerically for a high-Prandtl-number (Pr = 68) thermocapillary liquid bridge. The one-way-coupling approach is used for calculation of the particle motion, modeling PAS as an attractor for a single particle. The attractor is created by dissipative forces acting on the particle near the boundary due to the finite size of the particle. These forces can dramatically deflect the particle trajectory from a fluid pathline and transfer it to certain tubular flow structures, called Kolmogorov–Arnold–Moser (KAM) tori, in which the particle is focused and from which it might not escape anymore. The transfer of particles can take place if a KAM torus, which is a property of the flow without particles, enters the narrow boundary layer on the flow boundaries in which the particle experiences extra forces. Since the PAS obtained in this system depends mainly on the finite particle size, it can be classified as a finite-size coherent structure (FSCS).


Introduction
Small rigid neutrally-buoyant particles are frequently employed in experimental fluid mechanics to visualize the flow or to measure its velocity, assuming they move like the fluid. Such particles are called tracers. However, generally particle motion differs from that of the surrounding fluid and particles may cluster in certain regions of the flow. This has been extensively studied for particles transported in turbulent flows (see e.g. Peacock and Dabiri 2010;Haller 2015). In this study, on the other hand, we investigate small particles dispersed in a thermocapillary liquid bridge, in which the flow is laminar. In such system particles can cluster rapidly into particle accumulation structures (PAS) of curious shapes, a phenomenon discovered by Schwabe et al. (1996). Over the last two decades, PAS in thermocapillary flows has received increasing attention both experimentally Schwabe et al. 2007;Gotoda et al. 2016;Toyama et al. 2017;Gotoda et al. 2019) and numerically Mukin and Kuhlmann 2013;Muldoon and Kuhlmann 2013;Melnikov et al. 2013;Romanò and Kuhlmann 2018;Barmak et al. 2019). It is now clear ) that the rapidly formed PAS observed in experiments depends mainly on the small but finite size of the particle, which cannot be neglected during its motion near the free surface, and on the particular structure of the flow.
PAS in a liquid bridge can only be observed if the flow arises as a three-dimensional azimuthally traveling hydrothermal wave. This wave arises due to an instability of the steady axisymmetric flow (Preisser et al. 1983;Wanschura et al. 1995;Levenstam et al. 2001). Since all azimuthal harmonics of a traveling hydrothermal wave propagate with the same rotation rate (Leypoldt et al. 2000) this flow has the distinguished property to be steady in a frame of reference rotating with the azimuthal phase velocity of the wave. This property has important implications. Bajer (1994) showed that the trajectories of fluid elements in three-dimensional steady incompressible flows are equivalent, up to singular points, to the motion in phase space of a Hamiltonian system with one and a half degrees of freedom. Therefore, according to the theory of Hamiltonian systems, the streamlines in the stationary hydrothermal wave in the rotating frame of reference must be either regular or chaotic (Ottino 1989). Accordingly, a regular streamline winds on a torus, defining a closed streamtube, which corresponds to a socalled Kolmogorov-Arnold-Moser torus (KAM torus) of the theory of Hamiltonian systems (Schuster 2005). Like KAM tori, the closed streamtubes arise as densely nested sets. Typically, the KAM tori are surrounded by chaotic streamlines. A chaotic streamline will come arbitrarily close to any point, except for points inside the KAM tori. The part of the total volume occupied by the liquid which is filled with chaotic streamlines is often called chaotic sea. This chaotic sea of streamlines is typically found along the boundaries of the flow system, here of the liquid bridge. The union of all stationary streamlines (trajectories of passive fluid elements) defines a certain geometry. This geometry is also important for chaotic mixing (Aref 2002) and has been called kinematic template by Aref (1990).
Given the kinematic template made by the KAM structure of the streamlines in the rotating frame of reference, Hofmann and Kuhlmann (2011) have first explained PAS by assuming particles behave like tracers (fluid elements) in the bulk of the liquid bridge, but are transferred from chaotic to regular streamlines when they move near the free surface. Repeated visits to the free surface can then result in a focusing of particles to a periodic orbit (limit cycle) to form PAS. This particle transfer, called streamline hopping by Mukin and Kuhlmann (2013), is accomplished by extra forces acting on a finite-size particle near the boundary . The mechanism can only be operative if the distance between a KAM torus, or a closed streamline, and the free surface is of the same order of magnitude as the spatial range over which the extra forces acting on the particle are significant Mukin and Kuhlmann 2013;Romanò and Kuhlmann 2018). This spatial range is of the order of magnitude of the particle size. The relation between the particle size and the location of the closed streamline (KAM torus) determines the existence and the character of the PAS. Accordingly, PAS is an attraction of an individual finite-size particle to an attractor and not a collective (multi-particle) effect, at least during the initial phase of PAS formation. Therefore, the dynamics of PAS formation can be studied numerically by means of one-way-coupled simulations, in which a sufficient number of non-interacting particles are initialized at random positions in a given flow. The evolution of a particle ensemble is then monitored as a function of time under a suitable bulk transport model, typically either advection or the Maxey-Riley equation (Maxey and Riley 1983). PAS has formed if eventually many or all of the particles have clustered on a particular orbit which arrives very close to the free surface. In order to account for a finite-size effect near the flow boundaries the particlesurface interaction (PSI) model proposed by Hofmann and Kuhlmann (2011) and improved by Kuhlmann (2016, 2017) can be implemented. Adopting this model, which assumes particle-boundary interaction as an inelastic collision in the direction normal to boundary, PAS observed in ground experiments for Pr = 28 liquid bridges was successfully reproduced in numerical simulations by Kuhlmann (2018, 2017). PAS-related studies in low-and moderate-Prandtl-number (Pr ≤ 28) liquid bridges have been recently reviewed in . However, ground experimental data and numerical results are lacking for Pr = 68, corresponding to 5 cSt silicone oil to be used as a working fluid in the joint Japan-European Research Experiment on Marangoni Instability (JEREMI) planned to fly on the International Space Station (ISS). A few results previously obtained on the ISS in the framework of Marangoni Experiment in Space (MEIS) demonstrated the possibility of PAS for Pr = 68 under weightlessness (Matsumoto et al. 2014). The feasibility of PAS under microgravity conditions was also demonstrated by Schwabe et al. (2006) for Pr = 15 (ndecane). In the present study we aim to numerically predict PAS that can be expected for Pr = 68 in future space experiments.

Problem Formulation
We consider an incompressible Newtonian liquid (density f , kinematic viscosity ν, thermal diffusivity κ) held between two coaxial cylindrical rods of equal radius R, forming a liquid bridge (Fig. 1). The support rods are kept at a mutual distance d and have different constant temperatures T cold = T 0 and T hot = T 0 + T , so that T is the temperature difference. The fluid motion is driven by thermocapillary stresses caused by the spatial variation of the surface tension σ [T (x, t)] along the free surface (thermocapillary effect). The governing Navier-Stokes, continuity and energy equations (for details, see Hofmann and Kuhlmann 2011) are solved subject to no-slip/no-penetration/constanttemperature boundary conditions on the rods and imposing thermocapillary stresses and adiabatic conditions on the free surface, while neglecting viscous stresses and heat fluxes in the ambient atmosphere. The equations are rendered dimensionless by using the thermal-diffusive scales: d -for length, d 2 /κ -for time, κ/d -for velocity, and ρ f κ 2 /d 2 -for pressure. θ = (T − T 0 )/ T is the reduced temperature in dimensionless form. The fluid motion is then determined by the thermocapillary Reynolds number Re = γ T d/ρ f ν 2 , where γ = −∂σ/∂T | T =T 0 > 0 is the negative surfacetension coefficient, = d/R the aspect ratio of the liquid bridge, and Pr = ν/κ the Prandtl number. The latter is fixed to Pr = 68, corresponding to 5 cSt silicone oil (Shin-Etsu 2004). The critical Reynolds number defined by a linear stability analysis of the problem is Re c = 917 (Stojanovic and Kuhlmann 2020). We find the bifurcation to be supercritical and the flow arises as a three-dimensional traveling hydrothermal wave (HTW). Since all spectral components of a fully-developed HTW have the same azimuthal phase velocity = e z , the fully-developed flow field is steady in the frame of reference rotating with the same angular velocity. Thus the flow field in the rotating frame U(x) = u(x, t 0 ) − r e ϕ is steady and a single accurate snapshot of u(x, t) at t = t 0 together with is sufficient to characterize the flow.
To explore the Lagrangian topology of the steady flow in the rotating frame, streamlines and trajectories of fluid elements are obtained by integration of the advection equation which is invariant under the transformation to a rotating frame of reference with constant rotation rate. X(t) denotes the position of an infinitesimal fluid element at time t. Equation (1) must be solved for the initial condition X(t = 0) = X 0 specifying the initial position of the fluid element. We consider particles density-matched to the liquid which are so small that inertial effects are negligible. When such particle moves far from any boundary of the flow domain it behaves like a fluid element and moves on a streamline of the flow according to (1). However, the motion of the particle is hindered by its size near the flow boundaries. Therefore, the modified particle-surface interaction (PSI) model of Hofmann and Kuhlmann (2011) and Romanò and Kuhlmann (2018) is employed. Within this model the particle centroid can approach the boundary up to some distance = a + δ, comprising of the dimensionless particle radius a = a p /d and a dimensionless lubrication gap width δ. When the particle centroid approaches the boundary from the bulk, the velocity component normal to the boundary is annihilated at the distance from the boundary, and the particle slides parallel to the boundary until it reaches a point where the normal velocity turns inward. At this point the particle is released to the bulk, and it continues to be perfectly advected by the flow.
The existence of a lubrication gap on a smooth indeformable free surface (it is nearly indeformable due to the large reference surface tension σ (T 0 ) of 5 cSt silicone oil) that a perfectly wetted particle cannot penetrate is consistent with experimental observations ). The lubrication gap width δ to be used in this model depends on the flow parameters and on the particle size and density. It has been determined by  for a similar flow driven by surface stresses by matching the periodic orbits obtained by the PSI model to the periodic orbit of a particle which was computed by two-dimensional simulations of the Navier-Stokes and Newton's equations in which all scales were fully-resolved, including the flow structure in the lubrication layer. This approach was verified in  by successful comparison with experimental results for the limit cycle to which a single particle in the steady axisymmetric flow in a liquid bridge is attracted. Moreover, using the modified PSI model with the interaction parameter determined by the functional dependence proposed in Kuhlmann (2017, 2018) were able to satisfactorily reproduce PAS observed experimentally for Pr = 28. However, for Pr = 68 there is a lack of experimental data and no fully-resolved simulations have yet been carried out. Therefore, in order to identify possible PAS we vary the value of in the range [0.001,0.05] roughly corresponding to the particle radii of interest in experiments, providing that a ≤ ≤ 2a Kuhlmann 2017, 2018).

Numerical Methods
For the numerical simulations of the fluid flow a collocated finite-volume solver developed in the opensource software package OpenFOAM is employed. The numerical solver is based on the standard pressure-based solver pisoFoam for the Navier-Stokes and energy equations in the Boussinesq approximation, extended to include the thermocapillary stresses along the free surface. It implements the PISO (Pressure-Implicit with Splitting of Operators) algorithm with two external correction steps for the pressure to accomplish the pressure-velocity coupling, while treating the non-orthogonality of the computational grid by two additional internal corrections. Second-order schemes are utilized for both the spatial discretization (leastSquares for gradient terms, Gauss QUICK for divergence terms, and Gauss linear corrected for Laplacian terms) and for the time integration (implicit backward scheme). The resulting linear algebraic equations for the pressure-projection step are solved using the preconditioned conjugate gradient (PCG) method, whereas the preconditioned biconjugate gradient (PBICG) method is employed for the momentum and energy equations. The simulation proceeds to the next time step once the relative residuals of less than 10 −8 are reached for all dependent variables, i.e. pressure, velocity, and temperature.
In the recent study of Romanò and Kuhlmann (2018) for Pr = 28, a good agreement between the numerical results obtained with this solver and the ground experiment data of Toyama et al. (2017) was demonstrated. Additional details of the numerical implementation of the solver and its comprehensive validation can be found in Kuhlmann and Lemée (2016). All simulations are carried out with a mesh consisting of 21,504,000 cells. The mesh is refined near the free surface and the solid walls with the smallest cell sizes in the radial and axial directions being ≈ 0.0011 and ≈ 0.0006, respectively. The simulations advance with a small constant time step t = 10 −8 , ensuring the maximal Courant number based on the maximal local velocity is less than 0.1, until a fully-developed traveling HTW arises. To determine the time after which the HTW can be considered developed, the velocity and temperature are monitored as functions of time at sixteen points in the liquid bridge. The peak-to-peak frequency F of these signals determines the rotation rate = (2π/m) F of the wave, where m is the fundamental wave number. The HTW is considered fully developed, once the velocity changes between successive time steps have become less than 10 −5 ReP r t in the frame of reference rotating with the angular velocity . In this case the wave has become steady in the rotating frame of reference, the frequency F has become stationary, and it can be identified as the frequency of the fully-developed HTW.
Once the flow field is obtained, streamlines are computed in the rotating frame by numerically integrating (1) using the Runge-Kutta Dormand-Prince method (Dormand and Prince 1980). The absolute and relative numerical errors estimated as a difference between the fourth-and fifth-order accurate solutions are kept less than 10 −9 by an adaptive time-stepping. An additional approximation is introduced in the computation of the streamlines by the need to interpolate the discrete velocity field U to arbitrary points in the volume. This is accomplished by an interpolation with the same order of accuracy as the finite-volume solver implemented in OpenFOAM.
Particle trajectories are computed by solving (1) in the bulk with a classical 4th-order Runge-Kutta method. It was checked that the numerical dissipation has a negligible effect on the results when using a constant time step t = 10 −6 and integrating for no longer than t = 5 in thermaldiffusive units of time. Moreover, to reduce the computational cost of the simulations, the steady-state velocity field U is interpolated on a new structured computational grid in cylindrical coordinates, which consists of 28,848,000 grid points with (N r , N ϕ , N z ) = (200, 240, 601). The viscous and thermal boundary layers are taken care of by refining the mesh near the support rods and the free-surface. Since the structured grid in cylindrical coordinates exactly fits the cylindrical shape of the liquid bridge, the PSI model for the particle motion near the free surface can be implemented without further approximation.

Flow Topology for Pr = 68
A fully-developed traveling HTW has been obtained for = 1 and Re = 1500. For these flow parameters, the HTW has the fundamental wave number m = 1 and an angular velocity of = 835.5 in the units of κ/d 2 . Owing to the importance of regular streamlines for PAS Mukin and Kuhlmann 2013;Muldoon and Kuhlmann 2013) the streamline topology of the steady flow U(x) is analyzed with focus on streamlines inside KAM tori.
Regular and chaotic streamlines are, indeed, found to coexist for the case considered. The high spatial resolution of the simulation enables to identify various sets of KAM tori. In order to uniquely define the location of the closed streamlines and KAM tori their phase relative to the phase of the HTW needs to be determined. To that end the azimuthal angle ϕ is defined relative to the angleφ θ max at which the reduced temperature θ on the free-surface r = 1/ in the mid-plane z = 0.5 reaches its global maximum θ max Fig. 2 Temperature distribution θ(r = 1/ ,φ, z = 0.5) on the free surface in the mid-plane for = 1 and Re = 1500. Hereφ is the phase angle in the rotating frame up to an arbitrary phase in the rotating frame. This condition translates to θ max = maxφ θ(r = 1/ ,φ, z = 0.5) = θ(r = 1/ , ϕ := 0, z = 0.5), whereφ is the original azimuthal coordinate which contains an arbitrary shift. This is illustrated in Fig. 2 showing the free-surface temperature along the circle at mid-plane.
The flow topology is analyzed using a Poincaré section on the plane defined by ϕ = −0.35218 and −0.35218 + π . It is shown in Fig. 3. The coordinate x denotes the horizontal coordinate in the Poincaré plane with ϕ = −0.35218 being equivalent to the positive half-axis. Two major regular regions (KAM structures) shown in color can be identified. They are organized around two central closed streamlines (L 1a and L 1b ) which wind once about the z axis, as expected for a flow with m = 1. The subscript i of L i (streamline) and T i (KAM torus) stands for the number of azimuthal periods of the structure, while an additional letter distinguishes between structures with the same azimuthal period. The system of KAM tori which encompasses L 1b is extremely squeezed and stretched near the free surface as can be seen in the upper right (hot) corner of Fig. 3. These KAM tori also partly wrap about the system of KAM tori associated with L 1a and L 2 . Axial projections as well as three-dimensional views of the closed streamlines L 1b and L 2 are shown in Fig. 4 together with the Poincaré section at ϕ = −0.35218 and −0.35218 + π . The closed streamline L 1b in Fig. 4a approaches the free surface relatively closely, up to a distance FS = 0.02794 (see Table 1). This is a hint at its potential importance for PAS when suitably-sized particles are added to the liquid (see Section "Particle Accumulation for Pr = 68"). There also exist more complex closed streamlines. An example is L 7 which closes after seven revolutions about the z axis. It is shown in Fig. 5. These examples illustrate the complexity of the streamline structures for high-Prandtl number (Pr = 68) hydrothermal waves. The intricate streamline structure also demonstrates the high numerical accuracy required, and achieved, in order to resolve it. In particular, the residual divergence error of the numerically calculated flow field must be minimized to find closed streamlines which are periodic, under integration, for a sufficiently long time.
In order to quantify the properties of the flow topology for = 1 and Re = 1500, characteristic data of the most prominent closed streamlines (L) and KAM tori (T ) are gathered in Table 1. Given are the closest approaches of the closed streamline/KAM torus to the free surface  FS and to the hot solid wall W . The cold wall is never approached by the regular streamlines as close as the hot wall. Therefore, these data are not provided. In addition, the locations of the closed streamlines are specified by providing one intersection point (r 0 , z 0 ) (fixed point) in the half-plane ϕ = −0.35218 for each streamline. In case of a torus, we specify a point which defines the torus by the streamline originating from this point. As can be seen from Table 1 Characteristic data for some of the closed streamlines L and KAM tori T near which PAS is found. The coordinates (r 0 , z 0 ) in the half-plane ϕ = −0.35218 define either a closed streamline or a streamline on the largest reconstructible KAM torus. The anglẽ ϕ θmax = −1.0098 characterizes the location of the plane ϕ = 0 with respect to the original computed flow field in the rotating frame.
= 1, Re = 1500 and Pr = 68   Table 1 the closed streamlines and KAM tori are located much closer to the free surface than to the hot wall.

Particle Accumulation for Pr = 68
In this study PAS is investigated for small spherical particles density matched to the fluid, i.e. = ρ p /ρ f = 1 and particle Stokes number St = (2/9) a 2 O(10 −5 ). For such particles suspended in the liquid with the flow velocity, inertia has a negligible effect on the PAS formation as was demonstrated in Muldoon and Kuhlmann (2016) and Romanò and Kuhlmann (2018). This is confirmed for the motion of density-matched particles considered. For instance, the trajectory X(t) of a particle with = 1 and St = 10 −5 initialized at t = 0 velocity matched to the flow at x = (0.26031, 0.77518, 0.84837) and advanced using the simplified Maxey-Riley equation (Babiano et al. 2000) deviates, at t = 5 in thermal-diffusive units of time, only by 0.001 from the trajectory of a particle advanced by the advection (1) and initiated at the same point. In fact, the PAS resulting from the two models are visually indistinguishable. Therefore, retaining the inertia term in the transport model is not required and the particle motion in the bulk of the liquid bridge can be obtained integrating (1), while near the free surface and the solid walls the particlesurface interaction (PSI) model of Hofmann and Kuhlmann (2011) is implemented. In this approach, the only parameter characterizing the particles is the thickness of the layer on the boundaries (the free surface and the solid walls) which is not accessible by the centroids of the particles. Since the PAS obtained this way is formed solely due to the finite particle size, it represents a finite-size coherent structure (FSCS) . To improve the statistics, the motion of 4,000 non-interacting particles is computed. These particles are initialized at t = 0 at random positions within the volume accessible to the particle centroids.
As an example, the axial projection of the particle distribution for = 1, Re = 1500 and = 0.037 is shown in Fig. 6 as it evolves in time. Particles which have undergone at least one collision with the free surface or the solid walls (particle-surface interaction, PSI) are colored red, while particles which have not interacted with the boundaries up to the time displayed are colored blue. The time t κ which has passed since the initialization of the particle motion is given in units of the thermal diffusion time d 2 /κ. Also given is the time t ν in units of the viscous diffusion time d 2 /ν (t ν = Pr t κ ). Note that the dots indicate the particle centroid, but not the particle size. The dynamics of clustering can be quantified by the box counting measure K(t) . This measure is defined by dividing the liquid bridge into N cells cells of equal volume. For a given total number N of particles the accumulation measure K(t) can then be defined as the normalized sum over all cells of the deviations of the number of particles k i (t) in each cell from the average number of particles in each cellN = N/N cells ( 2 ) For convenience we selected the number of cells equal to the number of particles N, thusN = N/N cells = 1. We use 10 × 40 × 10 cells in the radial, azimuthal, and axial directions, respectively. The widths of the cells in Fig. 6 Time evolution of the distribution of 4000 non-interacting particles for = 1, Re = 1500 and = 0.037 shown in axial projection. Red dot: particle has undergone at least one PSI, blue dot: particle has not undergone any PSI. The dot indicates the particle centroid, not the particle size. Time is indicated in the subcaptions the azimuthal and axial directions are ϕ = 2π/40 and z = (1 − 2 )/10, respectively. Following Muldoon and Kuhlmann (2013), in order to satisfy the requirement of equal cell volumes, the radial cell length is r i+1 − r i = (1/ − ) 2 /10 + r 2 i 1/2 − r i , i = 1, ..., 10.
The accumulation measure K(t) for the evolution of the particle ensemble shown in Fig. 6 is presented in Fig.  7 (solid blue line). The dashed red line in Fig. 7 shows the fraction N red /N of particles which, by time t, have undergone at least one collision with the free surface. The initial random distribution with K(0) = e −1 ≈ 0.37 rapidly increases within t κ < 0.05 indicating clustering. This is due to the fact that many particles from the region near the axis of the liquid bridge are directly transported to the free surface where they are removed from their streamline by the first collision. Most particles become colored red in this early phase of evolution. As a result of this process the central region of the liquid bridge becomes rapidly depleted of particles (Figs. 6b and 7). The further gradual evolution of K(t) to K(t → ∞) ≈ 0.91 is caused by the slower attraction of most of the particles to a periodic orbit by multiple PSIs. For the interaction length = 0.037 the periodic orbit practically coincides with the closed streamline L 2 of the flow (compare Figs. 6d with fig. 4b).
Some particles (≈ 7% of all particles) never approach the free surface sufficiently close to undergo a PSI, since they have been initialized in those KAM tori which are always located further from the free surface than = 0.037. These KAM tori may also house isolated regions of chaotic streamlines. Figure 8 shows examples of PAS for = 1 and Re = 1500 in axial and azimuthal projections for selected interaction parameters at time t κ = 5. This does not mean that t κ = 5 is required for PAS to form. In fact, particle accumulation occurs on a shorter time scale for many cases.  Examples of PAS in a thermocapillary liquid bridge for = 1 and Re = 1500 at time t κ = 5 (t ν = 340). Shown are axial (left) and azimuthal (right) projections of the particle configuration for different values of the interaction parameter . Red dot: particle has undergone at least one PSI, blue dot: particle has never contacted the free surface. The dot indicates the particle centroid, not the particle size This can be seen in Fig. 6, where by t κ = 0.5 most of the particles (≈ 93%) cluster around L 2 . By comparing PAS for some interaction parameters, e.g. for = 0.027 (b) and = 0.041 (c), with the closed streamlines shown in Figs. 4a and 5, respectively, it becomes evident that PAS preferentially forms near the closed streamlines (and, more generally, on KAM tori) that enter the narrow boundary layer of the thickness on the free surface. However, PAS can be formed also in the chaotic region, where the particle trajectories are closed by PSI on the free surface, even though the chaotic streamlines are open . Such a structure arises for = 0.0082 and it is shown in Fig. 8a. In this case ≈ 30% of all particles are still colored blue at t κ = 5, i.e. they have never interacted with the free surface and, therefore, have remained in the chaotic toroidal core of the flow. The larger the interaction parameter (larger particles), the more particles undergo interactions with the free surface and the more particles are depleted from the central region of the liquid bridge. In case of PAS these particles end up on the periodic attractors (compare Figs. 8a and 8c).
The finite particle size is crucial for the accumulation of particles with = 1. In a limiting case = 0 (assuming an infinitesimally small particle size and thereby disregarding PSI) small density-matched particles are simply advected by the flow and keep moving along the streamlines even in the immediate vicinity of the boundaries. The same applies to the Maxey-Riley equation (Babiano et al. 2000) which reduces to the advection equation if the densitymatched particles are also initially velocity matched to the flow. Therefore, there is no mechanism which can transfer particles from one streamline to another and the particles will remain randomly distributed.

Conclusions
Accumulation of small spherical particles density matched to the flow has been demonstrated in Pr = 68 liquid bridges. The present modeling approach uses one-way coupling of non-interacting particles, whose motion has been accurately computed in a highly-resolved periodic traveling hydrothermal wave, which is steady in the rotating frame of reference. The flow topology has been analyzed in order to identify its regular streamlines. Some of the closed streamlines which approach the free surface sufficiently close have been found to act as organizing centers for density-matched particles of particular sizes.
The KAM structure for = 1 and Re = 1500 is found to be much more intricate than the one for Pr = 4 (Mukin and Kuhlmann 2013) and Pr = 28 (Romanò and Kuhlmann 2018). It allows for particles of many different sizes to cluster in a rich variety of PAS of different shapes, which can be well understood based on the underlying KAM template.
Although the PSI model has proven successful in capturing the key effect of the finite particle size near the flow boundaries and in predicting PAS, the global parameter remains undetermined in this approach. It would be desirable to determine by additional fullyresolving numerical simulation ), depending on the particle size and, in general, also on the particle density. Another possibility is a softening of the inelastic collision underlying the PSI model by, e.g., incorporating known asymptotic solutions such as the ones of Brenner (1961) for a particle in Stokes flow moving near a boundary. A similar approach has recently been used by , who modeled the particle-boundary interaction by the leading-order lubrication approximation for a particle in Stokes flow. Furthermore, the particles will interact with each other during the final stage of PAS when many particles accumulate on a periodic orbit. Efficient models are currently being developed to take this interaction into account.
Acknowledgements Support by ESA through contract 4000121111/ 17/NL/PG/pt is gratefully acknowledged. The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC).

Funding Open Access funding provided by TU Wien (TUW)
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.