1. Introduction
The dynamics of small inclusions in incompressible flows have attracted much attention over recent years. Part of this interest originates from the relevance of this subject to a plethora of fields both in nature and technology. There is indeed a variety of problems (Lyubimov, Lyubimova & Straube Reference Lyubimov, Lyubimova and Straube2005; Ehsan et al. Reference Ehsan, Abbasi, Ohringer and Parviz2006; Balboa Usabiaga & Delgado-Buscalioni Reference Balboa Usabiaga and Delgado-Buscalioni2011; Pushkin, Melnikov & Shevtsova Reference Pushkin, Melnikov and Shevtsova2011; Melnikov & Shevtsova Reference Melnikov and Shevtsova2017; Saghir & Abdulmajeed Reference Saghir and Abdulmajeed2017; Lappa Reference Lappa2018; Vorobev & Khlebnikova Reference Vorobev and Khlebnikova2018; Aksouh et al. Reference Aksouh, Chemini, Mataoui and Poncet2020; Vorobev, Zagvozkin & Lyubimova Reference Vorobev, Zagvozkin and Lyubimova2020), differing in regard to the spatio-temporal nature of the flow (laminar or turbulent) and/or with respect to the physical effect driving it, where dispersed particles can interact with a carrier fluid and produce different phenomena of practical interest. The study of this subject stands at the crossroads of physics, chemistry, biology, materials science, mechanical engineering and the science of soft matter.
Apart from these rather technical reasons, further research in this field has been spurred by purely theoretical factors, such as the quest to identify common elegant principles underpinning these phenomena (with the key aim to form predictive links between flow properties and particle behaviour during process operations) or the desire to build mathematical tools and specific numerical techniques that can be used to get insights into them.
This interest is still alive as witnessed by the several analyses being produced on this subject nowadays. It covers an extremely vast field encompassing different ‘variants’ such as particles transported by flows produced by pressure gradients (forced convection, see, e.g. Matas, Morris & Guazzelli Reference Matas, Morris and Guazzelli2003; Lashgari, Picano & Brandt Reference Lashgari, Picano and Brandt2015), differential rotation (Subbotin & Kozlov Reference Subbotin and Kozlov2020) or induced by temperature gradients such as thermogravitational (Akbar, Rahman & Ghiaasiaan Reference Akbar, Rahman and Ghiaasiaan2009; Puragliesi et al. Reference Puragliesi, Dehbi, Leriche, Soldati and Deville2011), thermocapillary (Melnikov & Shevtsova Reference Melnikov and Shevtsova2017; Gotoda et al. Reference Gotoda, Toyama, Ishimura, Sano, Suzuki, Kaneko and Ueno2019) and thermovibrational convection (Lappa Reference Lappa2016a, Reference Lappa2019).
Transcending specific differences among all these situations, it has been recognised that, in the absence of inter-particle forces (of hydrodynamic or electrostatic nature), the dynamics are essentially driven by the peculiar ability of particles with finite size and mass to violate the incompressibility constraint that is typical of liquids. Put mathematically, this implies that a clear distinction must be introduced between the velocity of the carrier fluid and that of the particles. While the former is solenoidal (i.e. its divergence is zero), this property does not hold for the latter, which formally makes a swarm of moving particles more similar to a ‘compressible gas’ than to a liquid (Pushkin et al. Reference Pushkin, Melnikov and Shevtsova2011).
As a result, the spacing among particles can change in time leading to nonlinear patterning behaviours of great interest (such as the emergence of depleted areas or regions where particles tend to ‘concentrate’). Though some of these phenomena can be organised in well-studied universality classes (Lyubimov et al. Reference Lyubimov, Lyubimova and Straube2005; Pushkin et al. Reference Pushkin, Melnikov and Shevtsova2011; Romanò & Kuhlmann Reference Romanò and Kuhlmann2018), however, often the interpretation of experimental realizations is not as straightforward as one would imagine, due to the inherent complexity of an experimental set-up compared to the simplicity of theoretical models and many other out-of-control factors.
The present work, in particular, is concerned with the theoretical and numerical analysis of such phenomena in the so-called ‘liquid bridge’ (LB). Landmark experiments in this configuration (a drop of liquid held by surface tension between two circular rods) date back to 1996 when Schwabe, Hintz & Frank (Reference Schwabe, Hintz and Frank1996) discovered the existence of peculiar structures formed by the accumulation of dispersed solid matter. These structures or ‘attractors’ were originally observed for conditions in which the Marangoni flow is dominant with respect to the buoyancy driven (thermogravitational) convection. Experiments were conducted in normal gravity conditions with liquid bridges having axial extensions of a few millimetres. More specifically, Schwabe et al. (Reference Schwabe, Hintz and Frank1996) observed them in the presence of travelling (hydrothermal) waves. Later Tanaka et al. (Reference Tanaka, Kawamura, Ueno and Schwabe2006) and Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007) characterised more precisely this phenomenon reporting on the existence of particle aggregates perfectly aligned in space along paths with the topology of a multi-blade circuit or windmill (today generally referred to as particle accumulation structures or PAS).
An extended description of the PAS morphology can be found in various papers appearing on the subject over the last decade (Tanaka et al. Reference Tanaka, Kawamura, Ueno and Schwabe2006; Abe, Ueno & Kawamura Reference Abe, Ueno and Kawamura2007; Ueno et al. Reference Ueno, Abe, Noguchi and Kawamura2008; Niigaki & Ueno Reference Niigaki and Ueno2012; Lappa Reference Lappa2013a,Reference Lappab,Reference Lappac; Gotoda et al. Reference Gotoda, Sano, Kaneko and Ueno2015; Toyama et al. Reference Toyama, Gotoda, Kaneko and Ueno2017; Oba et al. Reference Oba, Toyama, Hori and Ueno2019; Yamaguchi, Hori & Ueno Reference Yamaguchi, Hori and Ueno2019). In particular, one of the properties that has received significant interest is the so-called number of loops (
$N$), i.e. the times that the one-dimensional structure formed by the particles wounds about the vortex core. The value of 
$N$ has been found to be an integer multiple of the azimuthal wavenumber 
$m$ (
$N=m$ or 
$N=2m$) depending on several parameters in a non-monotone (often intermittent) way. For this reason some efforts have been devoted to a reconstruction of the related basin of attraction (Gotoda et al. Reference Gotoda, Sano, Kaneko and Ueno2015; Toyama et al. Reference Toyama, Gotoda, Kaneko and Ueno2017), i.e. the set of initial conditions or influential factors leading the system to meet one or the other family of attractors.
Remarkably, a significant part of this research has also been aimed to forge a new unified concept for the explanation of these phenomena; this has been attempted through the systematic integration into related models of more physics via approaches with different (experimental, theoretical and numerical) foci (Kuhlmann et al. Reference Kuhlmann, Lappa, Melnikov, Mukin, Muldoon, Pushkin, Shevtsova and Ueno2014; Lappa Reference Lappa2016b).
All these efforts have led to the emergence of two main theories evolving in parallel in the literature, based on the concepts of ‘inertial’ or ‘finite-size’ Lagrangian coherent structures, respectively (Melnikov, Pushkin & Shevtsova Reference Melnikov, Pushkin and Shevtsova2011; Pushkin et al. Reference Pushkin, Melnikov and Shevtsova2011; Lappa Reference Lappa2013a,Reference Lappab; Melnikov, Pushkin & Shevtsova Reference Melnikov, Pushkin and Shevtsova2013; Mukin & Kuhlmann Reference Mukin and Kuhlmann2013; Muldoon & Kuhlmann Reference Muldoon and Kuhlmann2013; Lappa Reference Lappa2014; Romanò & Kuhlmann Reference Romanò and Kuhlmann2018).
As implicit in its denomination, the first model or paradigm is essentially based on the idea that inertial effects (essentially due to the different densities of the dispersed solid phase and the surrounding fluid) can cause minute changes in the trajectories of particles with respect to the streamlines of the underlying fluid flow. In turn, these changes can accumulate in time and cause a variation in the spacing among particles leading them to cluster in some specific regions. According to the model originally developed by Pushkin et al. (Reference Pushkin, Melnikov and Shevtsova2011), in particular, this process would be mediated by the time-periodic nature of the flow, able to cause a synchronization between the hydrothermal wave produced as a result of an instability of Marangoni flow and the revolution motion of particles about the Marangoni toroidal vortex. This effect would be responsible for the correspondence between some geometrical features of the PAS (the aforementioned ‘blades’) and the nodes of the hydrothermal disturbance travelling in the azimuthal direction (leading to a phenomenon similar to that of the Mexican wave ‘La Ola’, where the wave propagation is made visible by the synchronised motions of the spectators, Melnikov et al. Reference Melnikov, Pushkin and Shevtsova2013). In line with this theory, by applying particle tracking velocimetry, Gotoda et al. (Reference Gotoda, Toyama, Ishimura, Sano, Suzuki, Kaneko and Ueno2019) have recently verified that the ratio of the hydrothermal wave frequency to the turnover frequency of the particle is independent of the intensity of the Marangoni flow.
According to the alternate model proposed by Kuhlmann and co-workers (Mukin & Kuhlmann Reference Mukin and Kuhlmann2013; Muldoon & Kuhlmann Reference Muldoon and Kuhlmann2013; Romanò & Kuhlmann Reference Romanò and Kuhlmann2018), the node-blade correspondence would be of a purely geometrical nature. On the basis of this theory, the entire PAS formation process would be driven by the presence of Kolmogorov–Arnold–Moser (KAM) tori of the flow, i.e. closed stream tubes existing in the frame of reference rotating with the same angular frequency of the wave, in which the flow field is steady. With this model, the accumulation of particles inside these tori would be made possible by virtue of their interaction with the free liquid surface (forcing particles to leave their original streamlines as a result of their finite diameter).
Both models have displayed a varying degree of success in explaining experimental dynamics depending on the considered circumstances. Hypotheses that have been put forward to resolve the discrepancy between the two theories are based on the idea that particle attractors exist in the fluid regardless of the inertia of particles and that different mechanisms can drive particles to them, one being the transversal migration (with respect to the direction of the fluid flow) which is induced by inertia in the presence of shear or fluid vorticity (Melnikov et al. Reference Melnikov, Pushkin and Shevtsova2011; Lappa Reference Lappa2013a,Reference Lappab, Reference Lappa2014), this mechanism being operative throughout the bulk of the fluid, and another being related to the aforementioned interaction of particles with the free surface. Despite these developments and the considerable attention these problems have received, however, a large number of key issues remain unresolved and many aspects still elude us.
Owing to the involved nature of these studies, we do not strive to fit an adequate account of them into the framework of the present introduction; rather, here we limit ourselves to emphasizing that, while the recent resurgence in interest in the identification of the abovementioned cause-and-effect relationships has led to significant advances in our understanding of the origin and properties of these particle structures, these in turn suggest new questions, both general and system specific, which require further analysis.
In particular, here we concentrate on an aspect that apparently has not been addressed yet, i.e. the role played in such phenomena by buoyancy effects. Gravitational forces acting on particles as a result of their different density with respect to the carrier fluid may indeed represent another possible mechanism forcing them to leave their original streamlines.
The present work examines the consequence of this weakness of the existing literature, which appears to have never been questioned to date. In this regard, our analysis is specifically motivated by (and builds on) the earlier experimental studies by Sasaki, Tanaka & Kawamura (Reference Sasaki, Tanaka and Kawamura2004) and Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007) and proposes widening the range of methodologies used to examine this specific problem.
The former authors reported on a new structure that seems to escape a simple classification on the basis of existing categorizations or paradigms for this type of phenomena. In order to clarify the underlying dynamics, we concentrate on non-neutrally buoyant spheres inside liquid bridges in normal gravity conditions. Rather than seeking further verification of already existing theories or models, we try to identify new physical effects. For these specific reasons (in order to filter out mechanisms that have been already characterised by other authors or theories), we address both cases with negligible and significant density difference.
2. Mathematical and numerical model
The schematic representation of the problem under consideration is depicted in figure 1. A liquid column of aspect ratio 
$AR = L/D$ is kept between two rods of diameter 
$D$ placed a distance 
$L$ apart. The two rods are heated differentially by keeping them at different temperatures, resulting in a temperature difference 
${\rm \Delta} {T} = T_{hot} - T_{cold}$. Absence of heat fluxes through the free surface (adiabatic condition) is assumed. Moreover, a hydrostatically deformed free surface is considered.
Figure 1. Sketch of a high aspect-ratio liquid bridge deformed under the influence of gravity.
Indeed, due to gravity, the free surface between the liquid and the surrounding (gaseous) environment can undergo departure from the cylindrical configuration (this is especially true for large-aspect-ratio liquid bridges, and the final interface shape results essentially from the balance between the weight of the liquid column and the magnitude of surface tension). The entity of the deformation is typically measured through the Bond number, 
, in which 
is the acceleration due to gravity while 
$\rho _0$ and 
$\sigma _{0}$ are the density of the liquid and the interfacial tension evaluated at the reference temperature, 
$T_0$, respectively.
Assuming that the capillary number is relatively small (see, e.g. Lappa Reference Lappa2004), the shape of the static free surface can be determined by solving the Young–Laplace equation
where 
${\rm \Delta} {p}$ is the pressure difference across the interface and 
$R_{1}$, 
$R_{2}$ are the principal radii of curvature. In particular, we have conveniently reformulated (
) in non-dimensional form in the cylindrical coordinates system shown in figure 1 by replacing the principal radii of curvature with their expression in the considered frame of reference. This has led to
\begin{equation} \frac{1}{c \left[ 1 + \left( \textrm{d} c/\textrm{d} z \right)^2 \right]} - \frac{\textrm{d} c/\textrm{d} z}{c \left[ 1 + \left( \textrm{d}^2c/\textrm{d} z^2 \right)^2 \right]^{3/2}} = k_{1}+k_{2}z, \end{equation}
to be solved with the boundary conditions 
$c = c( 0 ) = c( 1 ) = D/2L$, where 
$c = c( z )$ is the non-dimensional (with respect to the LB height, 
$L$) radial position of the interface. We have solved it numerically through a shooting method: the value of the parameter 
$k_{1}={\rm \Delta} {p}/\sigma _{0}$, which essentially determines the volume ratio 
$S=4V/{\rm \pi} D^2L$ (where V is the volume of the cylinder with diameter D and height L), has been guessed and changed in an iterative way up to obtain the desired value of 
$S$ (
$S=1$ in the present work), while 
$k_{2}=Bo$ has been assigned according to the considered conditions.
2.1. Mathematical description
Flow motion is governed by the incompressible Navier–Stokes and energy equations (2.3)–(2.5) that, taking into account the classical Boussinesq approximation, can be put in non-dimensional form as (see, e.g. Romanò & Kuhlmann Reference Romanò and Kuhlmann2018)
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}} = 0, \end{gather}
\begin{gather}\frac{\textrm{D}\boldsymbol{u}}{\textrm{D} t} ={-}\boldsymbol{\nabla}{p} +\boldsymbol{\nabla}^2{\boldsymbol{u}} + GrT\boldsymbol{e}_{z}, \end{gather}
\begin{gather}\frac{\textrm{D} T}{\textrm{D} t} = \frac{1}{Pr}{\nabla}^2{T}, \end{gather}to be solved together with the tangential stress boundary condition (2.6):
\begin{equation} \left[ \boldsymbol{\nabla}{\boldsymbol{u}} + \left(\boldsymbol{\nabla}{\boldsymbol{u}}\right)^{\mathrm{T}} \right] \boldsymbol{\cdot} \boldsymbol{n} + Re \boldsymbol{\nabla}_{s}{T} = 0. \end{equation}
Here, 
$\boldsymbol {u}$ is the fluid velocity, 
$\boldsymbol {e}_{z}$ is the unit vector in the 
$z$-axis direction (we assume that the liquid bridge axis is aligned with the gravitational field), the operator 
$\textrm {D}/\textrm {D} t$ indicates material derivation, the subscript s for the gradient operator indicates derivation along the free surface and 
$\boldsymbol {n}$ is the unit normal vector at the liquid–gas interface.
The dimensionless numbers appearing in the above set of equations are the Grashof number, 
, the Prandtl number, 
$Pr = \nu / \alpha$, and the Reynolds number, 
$Re = \sigma _{T}{\rm \Delta} {T} L/ \rho _{0} \nu ^2$, where 
$\alpha$ and 
$\beta$ are the thermal diffusivity and thermal expansion coefficient of the fluid, respectively, 
$\nu$ is its kinematic viscosity and 
$\sigma _T=-\partial \sigma /\partial T|_{T=T_0}$. For consistency with a companion work (Capobianchi & Lappa Reference Capobianchi and Lappa2020), in the following we shall refer to the Marangoni number, 
$Ma=RePr$, rather than 
$Re$.
Simulations have been carried out by solving the governing equations (2.3)–(2.5) using the computational platform OpenFOAM and the well-known segregated PISO algorithm of Issa (Reference Issa1986) (with the related collocated grid arrangement of variables). With such an approach, a numerical-oscillation-free solution is ensured by relying on the Rhie–Chow correction scheme (Rhie & Chow Reference Rhie and Chow1983). In particular, parabolic equations (momentum and energy) have been integrated in time by means of a first-order accurate implicit Euler scheme and the related convective terms have been treated through ‘Gauss integration’ (i.e. using the divergence theorem to interpolate the advected field at the cell faces). All convective terms have been discretised using the second-order accurate unbounded linear central differences. A separate discussion is needed for the elliptic pressure equation typical of the PISO approach. This equation stems from the intrinsic philosophy on which all pressure-velocity ‘projection’ methods for incompressible flow are based (i.e. it is obtained by substituting the velocity field ‘corrected’ by means of the gradient of pressure in the continuity equation). With OpenFOAM it is solved numerically by means of a standard iterative method and using homogeneous Neumann boundary conditions at the boundaries (in line with the consistent pressure Poisson equation strategy originally introduced by Gresho & Sani (Reference Gresho and Sani1987)). These ‘extra’ numerical conditions (generally referred to as ‘numerical’ boundary conditions to distinguish them from the PBC, i.e. physical boundary conditions which simply follow from the ‘physics’ of the considered problem) formally result from imposing at the boundaries the effective PBC for the velocity field; by doing so there is no need to correct the fluid velocity along the periphery of the considered fluid domain (as illustrated by Gresho & Sani (Reference Gresho and Sani1987), the problem can formally be made ‘well posed’ from a mathematical point of view by setting to zero the normal component of the pressure gradient).
2.2. Particle dynamics
As witnessed by existing numerical simulations on the subject, there is consensus in the literature (see, for instance, Schwabe et al. Reference Schwabe, Mizev, Udhayasankar and Tanaka2007; Kuhlmann et al. Reference Kuhlmann, Lappa, Melnikov, Mukin, Muldoon, Pushkin, Shevtsova and Ueno2014) that inter-particle forces do not interfere with (or do not contribute to) PAS formation in liquid bridges. This explains why the existing studies on this subject rely to an unusual level on insights obtained via theoretical models and ensuing computations that do not account for particle-to-particle (mutual) effects. The ‘back influence’ of solid particles on fluid flow is also generally neglected.
Following this common way of thinking, in the present work we resort to a simple and straightforward one-way coupling approach in which flow and particle equations are solved in a segregated, sequential manner; for the sake of completeness, however, in the present section some room is dedicated to the discussion of ‘critical’ (propaedeutical) arguments about the ‘validity’ of these assumptions.
In this regard, we start from the simple remark that, though (to the best of our knowledge) no studies have appeared in the literature about the limits of applicability of the one-way coupling approach either for generic laminar convection with dispersed particles or for the specific problem of PAS, some instructive information is available for the companion case of turbulent (incompressible) flows. As an example, a general rationale to distinguish ‘dilute’ and ‘dense’ suspensions (and identify suitable computational approaches accordingly) can be found in the review by Elghobashi (Reference Elghobashi1994). Following this classification (the reader being also referred to some related geometrical considerations elaborated in the book by Bodnár, Galdi & Nečasová (Reference Bodnár, Galdi and Nečasová2017)), one should primarily rely on the so-called ‘volume fraction’ of the dispersed phase (
$\phi$) as a relevant parameter to judge on the effectively needed level of coupling (i.e. the level of mathematical sophistication required to capture the ‘physics’ of particle laden flows). The parameter 
$\phi$ is generally defined as the ratio of the volume globally occupied by the solid particles in a certain sub-region of the considered fluid domain and the volume of the considered region itself (namely 
$\phi ={N_{part}{\rm \pi} d_p^3}/{6\varOmega }$, where 
$N_{part}$ is the number of solid particles and 
$\varOmega$ is the volume of the considered portion of fluid).
According to Elghobashi (Reference Elghobashi1994), for isotropic turbulence, one-way coupling should be assumed as a relevant approximation only for 
$\phi <10^{-6}$, larger values of this parameter (
$10^{-6}\leq \phi \leq 10^{-3}$ or 
$\phi >10^{-3}$), requiring two- (back influence of solid particles on fluid flow taken into account) or four-way coupling (fluid-dynamic interactions and collisions between the particles also properly accounted for), respectively. As a comparison with available experimental evidence for PAS would immediately reveal, however, these criteria should be regarded as overly conservative for simply time-periodic flows, this being witnessed by the very good agreement between existing numerical simulations based on the one-way coupling approach and the experimental results by Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007) (where the volume fraction based on the entire volume of the liquid bridge was 
$\phi \approx 10^{-4}$). The latter would implicitly suggest that the upper limit for the range of applicability of the one-way coupling strategy for these particular problems may be located at values of the volume fraction 
$\phi \geq 10^{-4}$.
Regardless of the specific type of flow (laminar or turbulent) and the effective threshold (in terms of 
$\phi$) to be exceeded to increase the needed level of coupling, these general considerations are meaningful or useful as they show that an ideal correspondence between any numerical study based on the one-way approximation and ‘effective experimental analysis’ (aimed to verify the validity of the numerical predictions) could easily be obtained by limiting the number of experiment tracers (that is, by making 
$\phi$ sufficiently small).
Before being carried away by the simplicity of this argument, anyhow, it should also be pointed out that its generality may be invalidated on ‘some spatial scales’ (even for extremely small values of 
$\phi$) by the peculiar nature of the PAS itself. This is because of the tendency of all particles to cluster in a string having negligible transverse extension and of the intrinsically ‘local’ nature of the geometrical criteria given above. Strictly speaking, these criteria should be verified in ‘any’ given neighbourhood of the PAS. Once again, a rationale for the widespread practice to give up on the modelling of particle interactions (four-way dynamics) can be found in available experimental observations or data. As illustrated, e.g. by Capobianchi & Lappa (Reference Capobianchi and Lappa2020), experiments dealing with a variety of fluids and particle types (see the literature cited in the introduction) have shown that once particles are inside the attractor (the KAM torus), even if 
$\phi$ increases locally to a significant extent, inter-particle forces or other ‘back influence’ effects are not intense enough to force particles to leave the attractor. This may be regarded as the sought justification about the effective possibility to neglect two- and four-way effects even when all the particles are concentrated in a thin region (the PAS itself). Obviously, to make the overall approach mathematically and physically consistent, the particle Stokes number defined as 
$St = d^2_{p} / 18L^2$, where 
$d_p$ is the diameter of the particle, must be much smaller than 1, i.e. 
$St \ll 1$. This condition is needed to make the Maxey–Riley equation (used to track particle motion) valid, and to satisfy the additional requirement that particle diameter must be smaller than the transverse size of KAM tori (to make the accumulation process possible).
Yet in line with the existing literature, we neglect the Saffman lift force and the Faxén corrections (Kuhlmann et al. Reference Kuhlmann, Lappa, Melnikov, Mukin, Muldoon, Pushkin, Shevtsova and Ueno2014). The influence of the Basset force is also disregarded on the basis of the earlier analysis by Lappa (Reference Lappa2013c) (where it was clarified that this force can play a non-negligible role only for non-dimensional angular frequencies of the flow exceeding 
$10^3$; the reader being also referred to the more recent study by Capobianchi & Lappa (Reference Capobianchi and Lappa2020)). Accordingly, the motion of each individual particle can be described by the abovementioned Maxey–Riley equation cast in the following condensed form (see, e.g. Babiano et al. Reference Babiano, Cartwright, Piro and Provenzale2000):
\begin{gather} \frac{\textrm{d}\boldsymbol{v}}{\textrm{d} t} = \frac{1}{\xi + 1/2} \left[\frac{3}{2} \frac{\textrm{D}\boldsymbol{u}}{\textrm{D} t} - \frac{\left( \boldsymbol{v} - \boldsymbol{u} \right)}{St} + \frac{\left(\xi-1 \right)\boldsymbol{e}_{z}}{Fr^2} \right], \end{gather}
\begin{gather}\frac{\textrm{d}\boldsymbol{x}}{\textrm{d} t}=\boldsymbol{v}. \end{gather}
The parameters appearing in (2.7) are the particle-to-fluid density ratio, 
$\xi = \rho _{p}/\rho$, the aforementioned Stokes number and the Froude number, 
$Fr = \nu /\sqrt {g{L^3}}$.
Unlike previous numerical studies (Romanò & Kuhlmann Reference Romanò and Kuhlmann2018; Capobianchi & Lappa Reference Capobianchi and Lappa2020), we do not model particle-boundaries interactions (with the only exception of the isodense case presented in § 3.2 for which otherwise no PAS would be obtained). In other words, for all the non-isodense particle cases presented in this work the centre of the particle is allowed to reach the limits of the domain. Notably, this choice should be regarded as a natural consequence of our (intentional) choice to explore new physical effects that may contribute ‘independently’ to the formation of PAS.
With the present model, only when the particle centre reaches a physical boundary, the normal component of the velocity is annihilated. Although this crude assumption does not reflect accurately the physics of the problem, it can be used to filter out the mechanisms theorised by Kuhlmann and co-workers. According to their theory, as the fluid inside KAM tori is locked and does not mix with the fluid outside them, these fluid-dynamic structures would play the role of templates for the accumulation of particles as opposed to the streamlines of chaotic nature located outside them. The mechanism for the transfer of particles from external streamlines to these loci would be favoured by the inability of particles to reach the free surface of the liquid bridge (due to their finite size and repulsive lubrication forces). Put simply, the interaction with the free surface would be responsible for an iterative (progressive) process by which particles are continuously transferred from one streamline to another until all particles fall inside the KAM tori.
Since particle interaction with the free surface is not modelled, PAS obtained in the present work in the presence of particle buoyancy (non-isodense cases, see § 3.3) should therefore be considered independent of this mechanism (though we acknowledge that in real experiments, the process based on the interaction of particles with the free interface might be effective).
3. Results
3.1. Base flow field
As already explained to a certain extent in the earlier sections, in this work we are interested in investigating (numerically) relatively overlooked, or not yet considered particle accumulation phenomena revealed by previous experiments.
Before embarking into such analysis, it is worth recalling briefly the main properties of Marangoni flow in liquid bridges. The typical hierarchy of bifurcations displayed by the liquid bridges is as follows: initially (for relatively small values of the imposed temperature difference) the flow is steady and axisymmetric (with fluid moving along the interface from the hot rod towards the cold rod and in the opposite direction in proximity to the axis, such a motion resulting in a single toroidal roll). On increasing the magnitude of the temperature gradient, however, the initial symmetry is lost and multi-cellular flows with well-defined patterns in the azimuthal direction are established after transients have decayed away (see, e.g. Shevtsova et al. (Reference Shevtsova, Mialdun, Kawamura, Ueno, Nishino and Lappa2011) and references therein). From a temporal point of view, in general, the initially steady flow is replaced by convective states that can be spatially periodic, quasi-periodic or chaotic (Swartenbroux & Schwabe Reference Swartenbroux and Schwabe1997; Ueno, Tanaka & Kawamura Reference Ueno, Tanaka and Kawamura2003; Melnikov, Shevtsova & Legros Reference Melnikov, Shevtsova and Legros2004, Reference Melnikov, Shevtsova and Legros2005). This line of research is extremely active as demonstrated by the studies still being produced on the subject (Nishino et al. Reference Nishino, Yano, Kawamura, Matsumoto, Ueno and Ermakov2015; Wang et al. Reference Wang, Wu, Duan and Kang2017; Kang et al. Reference Kang, Wu, Duan, Hu, Wang, Zhang and Hu2019).
The transition to time dependence is generally due to the onset and growth of disturbances which travel in the fluid (waves). Due to the axisymmetry of the interface separating the liquid from the external environment, these waves always emerge in pairs. Such coupled disturbances have azimuthal components of different sign (if one wave travels in the clockwise sense, the companion one will travel in the counter-clockwise sense). Moreover, their propagation direction is inclined with respect to the axis of the liquid bridge. As a result, the two waves produced by the instability have also an axial component (directed from the hot side to the cold one, i.e. in the upstream direction). The interplay or relative amplitude of the two opposite azimuthal components is known to make possible two types of waveforms, i.e. standing waves (equal amplitude) and travelling waves (different amplitude, Lappa (Reference Lappa2009) and references therein). Several experiments and numerical studies on the subject have indeed revealed one or the other mechanism of oscillation or a precise line of evolution where initial standing waves are taken over by travelling waves (Swartenbroux & Schwabe Reference Swartenbroux and Schwabe1997; Lappa, Savino & Monti Reference Lappa, Savino and Monti2001; Schwabe & Mizev Reference Schwabe and Mizev2011).
Here we concentrate on liquid bridges of high aspect ratio in normal gravity conditions. This configuration is of special interest due to the seemingly elusive nature of PAS that it can support. In fact, previous experiments of Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007) have indicated that the region of existence of accumulation structures, i.e. the range of 
$Ma$ for which well defined PAS are observed, is extremely narrow. This can immediately be appreciated by taking a look at figure 2. This plot shows the region of existence of PAS for different aspect ratios in normal gravity conditions observed by Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007). For relatively shallow liquid bridges, i.e. for those cases where the wavenumber of the hydrothermal wave is 
$m=2$ or 
$3$, the range where PAS can be obtained for a given aspect ratio is fairly wide. By contrast, for 
$m=1$ (
$AR\approx 0.9$), the region of existence of well forming structures reduces to a very thin interval.
Figure 2. Range of existence of well developing PAS determined experimentally by Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007). Data have been rescaled according to the present definition of aspect ratio and Marangoni number.
In the remainder of this paper we consider flow conditions corresponding to this limited range. Subsequently, investigation of PAS formation is carried out by varying parametrically the other main influential factors, i.e. 
$St$, 
$\xi$ and 
$Fr$.
As anticipated in the introduction, PAS generally display the topology of a ‘closed wire’ or windmill with many blades. As seen through an axial view (i.e. by looking at the ordered structure through one of the rods supporting the liquid bridge), PAS may also resemble a ‘flower’ with the centre corresponding to the axis of the liquid bridge and a series of petals evenly spaced in the azimuthal direction (their number depending on the azimuthal wavenumber). Put simply, for a ‘classical’ PAS the string of particles forms a closed spiral, wound 
$m$ times around the toroidal vortex and passing along the free surface of the liquid bridge (Schwabe et al. Reference Schwabe, Mizev, Udhayasankar and Tanaka2007).
Variants somehow representing a departure from this classical definition have also been identified. As an example, for 
$m\geq 2$, Muldoon & Kuhlmann (Reference Muldoon and Kuhlmann2013) and Melnikov & Shevtsova (Reference Melnikov and Shevtsova2017) have reported on asymmetric versions of PAS, i.e. windmills characterised by an eccentricity with respect to the axis of the liquid bridge.
A strong asymmetry is also a feature of the PAS originally observed by Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007) for 
$m=1$ (cf. figure 3a,b).
Figure 3. Particle accumulation structures determined experimentally by Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007). Flow conditions are: 
$AR=0.9$, 
$Ma\approx 18\,000$, 
$Gr\approx 6500$ and 
$Bo\approx 4.6$. Reprinted from ‘Formation of dynamic particle accumulation structures in oscillatory thermocapillary flow in liquid bridges’, by D. Schwabe, A.I. Mizev and M. Udhayasankar, 2007, Physics of Fluids, 19, 072102-8, with the permission of AIP Publishing.
The PAS for 
$m\geq 2$ (see, e.g. Lappa Reference Lappa2013c) are essentially centro-symmetric and tend to encapsulate the region of ‘cold fluid’ which is established inside the liquid bridge (in practice, by looking at the temperature distribution in the generic cross-section, this region appears as a ‘cold spot’ with the centre coincident with the axis of the liquid bridge). As shown in figure 4(a), for 
$m=1$, the region of cold fluid (roughly corresponding to the so-called ‘return flow’) is shifted away from the axis. There are some other notable differences that also deserve some attention. Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007) described the particle structure emerging in such conditions as ‘a deformed hook that seems to rotate around the axis of the LB, which is passing through its centre of mass’. Moreover, by inspecting the related experimental results reported in figure 3(a,b), the reader may get the impression that the PAS is not closed. From now on, we will refer to this peculiar structure as the H-PAS.
Figure 4. Temperature field in supercritical conditions: (a) above view of the plane cutting the LB at its mid-height, and (b) side view on a plane cutting the LB through its axis of symmetry.
For the present work, the ‘carrier’ flow (three-dimensional oscillatory flow to be used as a basis to explore PAS formation) has been obtained by carrying out simulations in conditions close to the experimental ones (for the convenience of the reader, the input conditions for our numerical study are reported in table 1). Assuming a diameter of the rods 
$D=6$ mm as in Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007), and considering an aspect ratio 
$AR=0.9$ and a 
${\rm \Delta} T=18$ K, this is equivalent in normal gravity conditions to 
$Ma\simeq 18\,000$, 
$Gr\simeq 6430$ and 
$Bo\simeq 4.6$. These are the parameters that we have therefore used for our numerical simulations, with the only exception of the Bond number for which we decided to use a slightly smaller value (i.e. 
$Bo\simeq 4$) because of numerical issues that we found in integrating (2.2).
Table 1. Physical properties of sodium nitrate measured at the working temperature adopted for the experiments of Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007).
We wish to highlight that another source of uncertainty relates to the nature of the liquid itself used for the experiments by Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007). Indeed, these authors mixed sodium nitrate with cesium nitrate, CsNO
$_3$, in order to match the density of the different tracers adopted for the experiments. Since we could not derive precisely the flow properties of the Na
$_{1-x}$Cs
$_x$NO
$_3$ melts resulting from this adulteration (the authors did not provide precise information about the fluid properties of their mixtures), some discrepancies should be expected between our simulations and the actual experiments (owing to the unknown variations potentially affecting 
$Pr$, 
$Gr$, 
$Bo$ and 
$Fr$).
For the numerical simulations, we have initially used a mesh (
$M_1$) consisting of 
$(40\times 44\times 40)$ cells along the three coordinate directions 
$(r,\varphi ,z )$. In agreement with the experiments, we have found the resulting flow to be unstable and characterised by the presence of a hydrothermal wave travelling in the azimuthal direction with mode 
$m=1$.
Subsequently, we have carried out a mesh convergence study by implementing two further levels of refinement, 
$M_2$, 
$M_3$. In particular, mesh 
$M_2$ was obtained by increasing the number of points by a factor 
$1.5$ along each direction while, for 
$M_3$, the points were increased by the same factor only in the axial direction, 
$z$ (with respect to 
$M_2$).
The corresponding periods of oscillation, 
${T}_1$, 
${T}_2$, 
${T}_3$ extrapolated from the signals (temperature) sampled at a point next to the interface, have been found to be 
${T}_1=3.17$, 
${T}_2=3.05$ and 
${T}_3=3.02$ s, for 
$M_1$, 
$M_2$ and 
$M_3$, respectively. The relative percentage error between two consecutive runs was 
$\varepsilon _{12}\simeq 3\,\%$ and 
$\varepsilon _{23}\simeq 0.7\,\%$. Given the good rate of convergence and the absence of significant differences between the last two runs, all the subsequent simulations have been carried out using the resolution 
$M_1$.
The period extrapolated from the material available in Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007) (cf. the metadata available in Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007) accessible from figure 6a,b) is 
$T\approx 2.8$ s. The difference with respect to our computations might be ascribed to the slightly different flow conditions considered in our simulations. Nevertheless, the agreement can be judged satisfactory. Further validation of the present approach has been obtained through comparison (not shown) with the results obtained by Melnikov et al. (Reference Melnikov, Pushkin and Shevtsova2011) and Lappa (Reference Lappa2013c) in terms of angular frequency of the hydrothermal wave and PAS emerging for 
$Pr=8$, 
$AR=0.34$, 
$Ma=20\,600$, 
$m=3$ (the percentage difference being smaller than 
$1\,\%$, Capobianchi & Lappa Reference Capobianchi and Lappa2020).
3.2. Particle accumulation for isodense systems ( 
$\xi =1$)
Having built the necessary base flow for the investigation of the emergence of PAS in the same conditions originally examined by Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007), we have initially considered isodense particles like those that were used in the experiments (
$\xi =1)$. As anticipated in § 2.2, for this specific case, the model about the interaction of particles with the free interface has been enabled (i.e. the centre of particles is intentionally prevented from reaching the free surface in order to allow PAS formation in the absence of inertia).
We wish to recall at this stage that elaborating a model able to capture all the physical aspects relating to the interaction of a solid particle with a free interface is not as straightforward as one would imagine. A realistic approach should take into account the wetting properties (contact angle) and surface deflections. Actually, a particle can touch the free surface, but in that case it is adsorbed there because of the capillary energy gain (the reader being referred to, e.g. Stocco & Nobili (Reference Stocco and Nobili2017) for an exhaustive elaboration of related ‘equilibrium’ arguments and the calculation of the free energy of the system as a function of the contact angle). In the present study we refer expressly to the modus operandi implemented by Hofmann & Kuhlmann (Reference Hofmann and Kuhlmann2011), which in turn has its roots in the earlier works by Vassileva et al. (Reference Vassileva, Van den Ende, Mugele and Mellema2005) and Do-Quang & Amberg (Reference Do-Quang and Amberg2009). This may be regarded as an elastic reflection model relying on the idea that ‘a small particle cannot pierce out of a wetting liquid due to the high capillary pressure that would be associated with a surface bulging on the small scale of the particle radius’ (Hofmann & Kuhlmann Reference Hofmann and Kuhlmann2011). The velocity of the particle perpendicular to the free surface is annihilated accordingly when the distance between the centre of the particle and the surface is equal to the particle radius. As noticed by these authors, this paradigm also aligns with the hypothesis of an undeformable interface, i.e. with the assumption of an asymptotically large mean surface tension (
$Ca\rightarrow 0$). It obviously implies the perfect wettability of the particle by the liquid (Hofmann & Kuhlmann Reference Hofmann and Kuhlmann2011).
Notably, using a slightly different perspective, this interaction scheme may also be regarded as a realisation of the treatment originally elaborated by Faxén (Reference Faxén1959) using the so-called ‘method of reflections’. Using this model (see also Brenner Reference Brenner1961), the motion of a single solid sphere approaching a free surface can be made equivalent to the case of two equal spheres moving toward each other. As a result of the hydrodynamic interaction of the solid sphere with its unphysical ‘mirror image’ (located beyond the free surface), the particle is not adsorbed but slips along the direction of the interface (Happel & Brenner Reference Happel and Brenner1991).
The present section, however, should be regarded as an attempt to go beyond a mere verification of the validity of the particle-free-surface theory for PAS formation elaborated by Kuhlmann and co-workers and relying on the approach discussed above. Indeed, some effort is also provided to show that while the interaction with the free interface can promote PAS emergence, other mechanisms can proceed in the opposite direction.
Along these lines, figure 5 reports in an ordered fashion the PAS produced by our simulations for different particle sizes (
$1$, 
$20$, 
$30$ and 
$40$ mm) and 
$Fr\simeq 1\times 10^{-3}$ (the value of the Froude number being however irrelevant for these simulations given the isodense nature of the dispersed solid mass). Although the findings collected in figure 5 show similar particle patterns, it appears clear that in some circumstances the structure is not well formed (especially for the case with 
$20$ microns particles) since a non-negligible amount of particles remains freely dispersed outside the PAS circuit. Leaving aside for a while the different degree of particle accumulation, however, it can be seen that each emerging structure is constituted by a central region that coils up and then extends from the base of the LB (cooled rod) toward the opposite side (heated rod) following a path that for an external observer would appear inclined with respect to the axis of the LB. The rest of the structure is essentially a curved branch that connects two nearly opposite areas of the heated rod.
Figure 5. Particle accumulation structure seen from different lateral perspectives (
$t_{1,2,3}$) obtained by turning the view by 120 degrees for the case 
$\xi =1$ and four different tracer sizes (each row corresponding to a different particle size). From the first to last row: 
$d_p=1$ mm (
$St\approx 2\times 10^{-9}$), 
$d_p=20$ mm (
$St\approx 7.6\times 10^{-7}$), 
$d_p=30$ mm (
$St\approx 1.7\times 10^{-6}$ and 
$d_p=40$ mm (
$St\approx 3\times 10^{-6}$)).
As comparison with figure 3(a,b) would reveal, the PAS calculated numerically have recognizable features in common with the H-PAS found by Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007). The H-PAS, in fact, manifests itself as an elongated area of accumulation that extends from one rod to the other, which in the lateral views seems to connect two opposite corners of the LB (cf. figure 3b). If looked at from above (figure 3a), it would resemble a spiral with relatively limited azimuthal extent (approximately one quarter of the overall circumferential extension of the liquid bridge free surface). Similar features are also recognizable in the PAS determined numerically whose peculiar topology can be inferred from figure 6.
Figure 6. Views (lateral and from above) of the H-PAS determined numerically for the case 
$d_p=1$ mm and of the corresponding closed stream tube (centreline of the KAM torus) determined in a frame of reference rotating at the angular frequency of the travelling wave but with opposite direction (in the view from above, the larger circle represents the periphery of the closest disk, while the smaller circle is periphery of the opposite disk; the shaded region corresponds to the axial projection of the interface).
Figure 6 is particularly instructive as it also reveals the closed stream tube responsible for the emergence of the PAS, namely the KAM torus observable in the frame of reference rotating with angular speed 
$-\varOmega$, where 
$\varOmega =2{\rm \pi} /\mathrm {T}$ is the angular frequency of the hydrothermal wave.
Notably, while this torus plays the role of backbone for PAS formation, the effective tendency of isodense particles to be entrained (or repelled) by it can be put in relation with their size (which can be used to explain the variable level of particle scattering visible in figure 5 for different diameters). As illustrated by several authors (Szeri, Wiggins & Leal Reference Szeri, Wiggins and Leal1991; Babiano et al. Reference Babiano, Cartwright, Piro and Provenzale2000; Vilela, de Moura & Grebogi Reference Vilela, de Moura and Grebogi2006; Sapsis & Haller Reference Sapsis and Haller2008), in particular, if the Stokes number is large enough, mechanisms are enabled that can prevent some particles from being captured by the attractor (Gotoda et al. Reference Gotoda, Melnikov, Ueno and Shevtsova2016). Such disturbing effects can be placed in a precise mathematical context by using existing criteria.
There is a significant amount of literature in this regard. As an example, Babiano et al. (Reference Babiano, Cartwright, Piro and Provenzale2000) could successfully characterise the unstable regions of the space of parameters in which scattering of inertial particles occurs in two-dimensional time-dependent flows. Lately, Sapsis & Haller (Reference Sapsis and Haller2008) tried to eliminate some limitations affecting the study by Babiano et al. (Reference Babiano, Cartwright, Piro and Provenzale2000) and formulated another criterion able to provide the limiting conditions for which finite-size neutrally buoyant particles might be repelled by three-dimensional regions of attraction (the so-called ‘slow manifold’), i.e. the KAM torus in the present case.
Here, in particular, we take advantage of such theoretical developments (Sapsis & Haller Reference Sapsis and Haller2008) to investigate the role of the particle size on the formation of H-PAS discussed before. The criterion is briefly summarised below for the convenience of the reader.
In the case of neutrally buoyant particles, the Maxey–Riley equation (2.7) becomes
\begin{equation} \epsilon\dot{\boldsymbol{v}}-\epsilon\frac{\textrm{D}\boldsymbol{u}}{\textrm{D} t}={-}\left(\boldsymbol{v}-\boldsymbol{u}\right), \end{equation}
where 
$\epsilon =3/2St_f$ (here, 
$St_f=StRe$ is the product of the Stokes number defined above with the Reynolds number); moreover, it implies that the material derivative
\begin{equation} \frac{\textrm{D}\boldsymbol{u}}{\textrm{D} t}=\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\cdot \boldsymbol{\nabla}{\boldsymbol{u}} \end{equation}is taken along the fluid path. Alternately, the derivative taken along the particle trajectory can be defined as
\begin{equation} \frac{\textrm{d}\boldsymbol{u}}{\textrm{d} t}=\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{v}\cdot \boldsymbol{\nabla}{\boldsymbol{u}}. \end{equation}
Under the assumption that particles follow precisely the flow, i.e. by considering 
$\textrm {D}\boldsymbol {u}/\textrm {D} t=\textrm {d}\boldsymbol {u}/\textrm {d} t$, (3.1) can be greatly simplified and, for such a case, the solution of the problem would simply read 
$\boldsymbol {v}-\boldsymbol {u}=(\boldsymbol {v}_0-\boldsymbol {u}_0)\exp (-t/\epsilon )$. Put differently, this indicates that regardless of the initial conditions, 
$\boldsymbol {v}_0, \boldsymbol {u}_0$, the particle velocity, 
$\boldsymbol {v}$, would approach the velocity of the fluid, 
$\boldsymbol {u}$, exponentially (i.e. particles would tend to accumulate on the aforementioned slow manifold, which in our case is represented by the KAM torus).
Babiano et al. (Reference Babiano, Cartwright, Piro and Provenzale2000), however, argued that in the presence of finite-size particles, 
$\textrm {D}\boldsymbol {u}/\textrm {D} t\neq \textrm {d}\boldsymbol {u}/\textrm {d} t$. By substituting (3.2), (3.3) into (3.1), Sapsis & Haller (Reference Sapsis and Haller2008) derived accordingly the following alternate formulation of the Lagrangian problem for neutrally buoyant particles:
\begin{gather} \frac{\textrm{d}}{\textrm{d} t}\left[\boldsymbol{v}-\boldsymbol{u}(\boldsymbol{x},t) \right]={-}\left[\boldsymbol{\nabla}{\boldsymbol{u}(\boldsymbol{x},t)+\boldsymbol{I}/\epsilon}\right] \left[\boldsymbol{v}-\boldsymbol{u}(\boldsymbol{x},t) \right], \end{gather}
\begin{gather}\frac{\textrm{d}\boldsymbol{x}}{\textrm{d} t}=\boldsymbol{v}. \end{gather}Sapsis & Haller (Reference Sapsis and Haller2008) were then able to establish a criterion that can be mathematically expressed as follows.
Theorem For any 
$\epsilon >0$, the invariant manifold, 
$M_{\epsilon }$, repels all close enough trajectories 
$\left [\boldsymbol {x}(t),\boldsymbol {v}(t)\right ]$ as long as they satisfy
\begin{equation} \lambda_{min}\left[\boldsymbol{S}(\boldsymbol{x}(t),t)+ \boldsymbol{I}/\epsilon\right]<0, \end{equation}
where the invariant manifold 
$M_{\epsilon }$ is defined, for the case of neutrally buoyant particles (see Sapsis & Haller Reference Sapsis and Haller2008), as
\begin{equation} M_{\epsilon}={(\boldsymbol{x},\phi,\boldsymbol{v}):\boldsymbol{v}=\boldsymbol{u}(\boldsymbol{x},\phi)}, \end{equation}
$\phi \equiv t_0+\epsilon \tau$ being a dummy variable which renders the above system (3.4) and (3.5) autonomous in the variable 
$(\boldsymbol {x},\phi ,\boldsymbol {v})$. Here 
$\boldsymbol {S}(\boldsymbol {x}(t),t)=\frac {1}{2}[\boldsymbol {\nabla }{\boldsymbol {u}}+\boldsymbol {\nabla }{\boldsymbol {u}}^\textrm {T}]$ and 
$\lambda _{min}$ are the smallest eigenvalues associated with the tensor between squared brackets in (3.6).
Briefly stated, the problem relating to the application of this theorem is to determine the effective spatial distribution of 
$\lambda _{min}$.
Along these lines, in the present work the eigenvalues of the tensor 
$\boldsymbol {S}(\boldsymbol {x}(t),t)+\boldsymbol {I}/\epsilon$ have been computed for the cases mentioned before, i.e. for tracers size of 
$1$, 
$20$, 
$30$ and 
$40$ microns.
In the first case 
$\lambda _{min}$ is always positive (hence the condition defined by the criterion of Sapsis & Haller (Reference Sapsis and Haller2008) is not fulfilled). On the other hand, for the 20 microns particles case, 
$\lambda _{min}$ becomes negative in two distinct areas next to the corners between the interface and the solid rods, as evident in figure 7. For what concerns the two remaining cases with larger particles, i.e. 30 and 40 microns, the region affected by negative values of the smallest eigenvalues becomes increasingly larger. On the basis of these findings, one would expect the PAS shown in figure 5 to become less defined for subsequent increments of the particle size. Nevertheless, as clearly witnessed by the evolution shown in figure 5, while the structure initially loses definition (case for 20 microns, particles shown in the second row), for the two subsequent cases it becomes increasingly more defined. We therefore argue that the two processes described above, namely, the interaction with the free surface (where the larger the size of the particles, the more effective the role played by the interface in leading them to meet the closed streamtube) and the destabilizing effects predicted by the criterion of Sapsis & Haller (Reference Sapsis and Haller2008) for increasing 
$St$, should be seen as two competing (opposite) mechanisms. Most remarkably, these arguments (being supported by the present computations) might be considered as a justification for the non-monotone behaviour observed here. As expected, on increasing the particle Stokes number beyond a certain limit (
$St\sim 7\times 10^{-6}$ in the present conditions), PAS formation is completely prevented (the destabilizing effect becomes indeed predominant).
Figure 7. Regions where the criterion of Sapsis & Haller (Reference Sapsis and Haller2008) is verified for the case with 20 microns particles.
For the sake of completeness, we wish also to highlight that the PAS formation time (evaluated as the number of periods of the hydrothermal wave (HTW) required to make the structure recognisable) was approximately the same for all the cases addressed in the present section, namely, 4–6. The images of the H-PAS shown in figure 5, however, have been taken in correspondence of a much larger number of periods (
$>$50 in all cases).
3.3. Particle accumulation structure for the cases $\xi <1$
After completing the study of PAS occurrence for isodense particles, in this section we turn to examining situations in which the particles are not isodense. As stated in § 2.2, for these cases we intentionally disable the model about the interaction of particles with the free surface. This means that we examine circumstances where only particle inertia and buoyancy (unless 
$Fr\to \infty$) might play a relevant (influential) role during the process of particle aggregation in coherent structures.
Our decision to rely on the set-up by Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007) as a testbed for the evaluation of these effects follows from the shortage of other relevant studies to be considered for such a purpose. In other words, without a clear indication of the region of existence of PAS, we base our investigation on a trial-and-error approach and change parametrically 
$St$ and 
$\xi$ until a PAS is observed.
As a first attempt, in particular, we have started our analysis considering relatively large tracers (
$St\approx 1.7\times 10^{-4}$) and varied 
$\xi$ in a neighbourhood of 
$\xi =1$. For 
$\xi$ much larger than 1, no PAS have been detected. In fact, for these cases large particles just sediment. Similarly, for particles slightly heavier than the fluid, we did not obtain any PAS even if particles could be retained in the fluid (by virtue of the drag forces exerted on them by the moving liquid, i.e. the Marangoni flow).
For the subsequent series of numerical experiments, again, first we have used the value of 
$St$ as above. Diversely, in these cases particles less dense than the fluid (
$\xi <1$) have been considered.
By seeding relatively light particles (
$\xi \approx 0.5$) in the fluid, only particle flotation has been observed (no PAS). Most interestingly, however, something rather unexpected has been found when setting the particles density very close to that of the fluid (yet for the cases 
$\xi <1$). In some specific conditions we came across an apparently heretofore unseen structure. In particular, figure 8(a) shows the structure found numerically for 
$\xi \simeq 0.987$ and 
$St\simeq 1.7\times 10^{-4}$. The first two snapshots refer to the top view and the view from below while the other frames (side views) provide an overview of the PAS from different perspectives. Next to our findings, figure 8(b) shows the structure experimentally observed by Sasaki et al. (Reference Sasaki, Tanaka and Kawamura2004); to the best of our knowledge this is the only work which has reported such a type of PAS. We wish to emphasise that, before coming across the numerically computed structure shown in figure 8(a), we were not aware of the experimental study by these authors.
Figure 8. Particle accumulation structures determined numerically (snapshots shown in a): views from above and from below, and lateral views at different time frames. Flow conditions were 
$\xi \simeq 0.987$ and 
$St\simeq 1.2\times 10^{-4}$. Experimental images of the PAS observed by Sasaki et al. (Reference Sasaki, Tanaka and Kawamura2004) at different time frames (b) reproduced with permission of the authors.
Interestingly, if the numerically obtained structure is examined from above (equivalent to an experimental observation through a transparent disk) two main (approximately circular) turns can be distinguished. As the larger radius turn corresponds to the lower spiral turn in the side view, we infer that this PAS pertains to the same category of structures discovered by Sasaki et al. (Reference Sasaki, Tanaka and Kawamura2004) (this conclusion being supported by the direct comparison with the corresponding experimental images). It can therefore be argued that the main mark distinguishing the two PAS (H-PAS and that shown in figure 8) types is the orientation of the turns present in the spiraling structure. While in the classical case (H-PAS) such turns extend essentially in the radial and axial directions, in the second case they are perpendicular to the liquid bridge axis, i.e. they tend to develop in the azimuthal direction.
As it seems to escape a simple classification on the basis of existing categorizations, in the following we will refer to this specific structure as the S-PAS.
At this stage it is worth recalling that, for situations in which 
$m\neq 1$, some authors have reported on the existence of closed loops where the number of windings was double with respect to the azimuthal wavenumber. This specific category of PAS has led the scientific community to partition the observed structures into two main classes. These have been named, SL-I (Abe et al. Reference Abe, Ueno and Kawamura2007; Ueno et al. Reference Ueno, Abe, Noguchi and Kawamura2008) and SL-II (Tanaka et al. Reference Tanaka, Kawamura, Ueno and Schwabe2006; Niigaki & Ueno Reference Niigaki and Ueno2012), respectively, to take into account the single or double nature of the winding process. Sasaki et al. (Reference Sasaki, Tanaka and Kawamura2004) concluded that the structure shown in figure 8 could be considered a SL-II type structure, which, however, we think is not completely correct (or it requires a different interpretation, as we will illustrate in the remainder of this study). Along these lines, as evident in figure 9, an ‘independent’ KAM torus (with respect to that depicted in figure 6) can be identified in the fluid-dynamic velocity field (in the relative frame of reference rotating with the same angular velocity of the hydrothermal wave). This KAM torus clearly plays the role of template for the accumulation of particles in this case, which indicates that structures like those shown in figures 6 and 9 may even coexist in the space of phases, i.e. they may be considered as ‘multiple solutions’ for PAS originating from 
$m=1$ fluid-dynamic states.
Figure 9. Particle accumulation structures observed for the case 
$\xi \simeq 0.987$ and 
$St\simeq 1.7\times 10^{-4}$ and relative KAM torus centreline (red line).
Following up on the previous point, it should be pointed out that the existence of multiple solutions, in turn, implicitly connects to another important problem, that is, the determination of the related basin of attraction, i.e. the set of initial conditions or influential factors leading the system to develop a specific structure.
Even though no information is available in this regard in the existing literature for 
$m=1$, some useful indications, however, can be found about the cases with 
$m\geq 2$ (Mukin & Kuhlmann Reference Mukin and Kuhlmann2013; Gotoda et al. Reference Gotoda, Melnikov, Ueno and Shevtsova2016; Melnikov & Shevtsova Reference Melnikov and Shevtsova2017). Indeed, these authors tried to reconstruct the basin of attraction for different families of attractors (the interested reader being also referred to the very interesting experiments conducted by Gotoda et al. Reference Gotoda, Sano, Kaneko and Ueno2015; Toyama et al. Reference Toyama, Gotoda, Kaneko and Ueno2017). The main outcomes of all these efforts is that the basin of attraction can encompass a wide range of factors such as (but not limited to) the particle/liquid density ratio (Schwabe et al. Reference Schwabe, Mizev, Udhayasankar and Tanaka2007; Gotoda et al. Reference Gotoda, Melnikov, Ueno and Shevtsova2016), the size of particles (Romanò & Kuhlmann Reference Romanò and Kuhlmann2018) or ‘initial position’ (Melnikov & Shevtsova Reference Melnikov and Shevtsova2017), the simultaneous presence of particles with different sizes (i.e. a multi-dispersed distribution, Gotoda et al. Reference Gotoda, Melnikov, Ueno and Shevtsova2016), and the value of the Marangoni number and the aspect ratio, which, in turn, determine the structure of the carrier (Marangoni) flow (Gotoda et al. Reference Gotoda, Sano, Kaneko and Ueno2015; Toyama et al. Reference Toyama, Gotoda, Kaneko and Ueno2017). A similar study about the region of existence of the S-PAS is presented in the next section.
3.3.1. Region of existence of the S-PAS
As discussed in the preceding section, we initially found the S-PAS for one specific combination of the pair 
$(St,\xi )$. Nevertheless, in order to determine its range of existence, we have explored different conditions in the space of parameters 
$(St,\xi )$ while maintaining 
$\xi <1$.
Initially, we have maintained unaltered the density ratio and run a series of simulations adjusting the Stokes number. For a relatively small particle size, S-PAS was still visible, although below a certain threshold the structure became gradually less definite (cf. figure 10) until it eventually disappeared. On the other hand, for moderate increments of 
$St$ with respect to the value reported before (cf. § 3.3), the structure was not captured at all.
Figure 10. Snapshots of S-PAS determined numerically with flow conditions 
$\xi \simeq 0.987$ and 
$St\simeq 1.2\times 10^{-4}$ (particle diameter 250 microns).
More interesting are those cases where we adjusted both 
$St$ and 
$\xi$ at the same time. After a number of attempts, we could successfully identify a well-defined region of existence for the S-PAS (cf. data displayed in figure 11). From these data some interesting insights can be drawn.
Figure 11. Region of existence of the S-PAS in the space 
$(St,\xi )$. (a) Linear-linear scale, (b) log-linear scale. The other flow parameters are the same as those considered in the previous calculations, where the S-PAS has been observed.
By looking at figure 11(a), an asymptotic behaviour can be noticed in the limit as 
$\xi \to 1$. This suggests that as the density of the particles becomes closer to that of the fluid, the particle diameter must become increasingly larger in order to allow the formation of the S-PAS. Although from a practical point of view excessively large values of 
$St$ are of no interest, this result provides some compelling evidence that the mechanism of formation is dictated solely by the competition between Stokes drag and buoyancy (for 
$\xi \to 1$, mass effects become irrelevant). Along these lines, we wish also to recall that, for exactly isodense systems, we came across the different accumulation structure (the H-PAS) discussed in § 3.2. Therefore, the vertical asymptote might also be regarded as a boundary separating the regions of occurrence of these two different categories of PAS. Figure 11(a) also shows that, in addition to the vertical asymptote, an upper limit also exists (the reader being referred to the black symbols and the related solid line introduced as a guide for the eye). As revealed by our simulations, for larger diameters (region located over the solid line), all the particles initially uniformly dispersed in the fluid simply accumulate at the top rod due to the increased buoyancy (which eventually prevails on the Stokes drag, thereby preventing PAS formation completely).
Additional useful information can be gathered from figure 11(b), where, in order to improve further the presentation of data and reveal additional aspects, we have used a logarithmic scale for the vertical 
$St$ axis. In particular, this companion figure is instrumental in showing that all the ‘numerically determined points’ are located within a region of existence that is delimited from above by the aforementioned hard boundary (corresponding to the cutoff particle diameter beyond which PAS are no longer allowed), while no precise (hard) boundary can be identified on its bottom. In place of such a boundary, an asymptotic condition can rather be defined where the PAS formation time grows indefinitely as the particle Stokes number is progressively decreased.
In this regard, figure 11 is naturally complemented (quantitatively substantiated) by figure 12 where we have reported the formation time for the S-PAS (measured in the number of periods of the HTW) as a function of the Stokes number for some representative values of the density ratio 
$\xi$. It can clearly be seen that, for each value of 
$\xi$, the formation time becomes significantly higher when 
$St$ is decreased (similar behaviours were observed experimentally by Schwabe et al. Reference Schwabe, Mizev, Udhayasankar and Tanaka2007). In particular, this trend becomes steeper as 
$\xi$ is made smaller (this being particularly evident for 
$\xi =0.75$ and 
$2.5 \times 10^{-6} < St < 5 \times 10^{-6}$). Vice versa, the minima of the three branches collected in this figure (for each of the curves shown in figure 12 the minimum formation time is attained in correspondence of the highest value of 
$St$) simply reflect the upper threshold indicated in figure 11 by the black symbols.
Figure 12. Time formation of the S-PAS measured in the number of periods of the HTW as a function of the Stokes number for different values of the density ratio 
$\xi$.
As the reader will easily realise by inspecting figure 12, the time needed for the emergence of the S-PAS in proximity to the upper boundary of its region of existence is roughly the same time needed to obtain the H-PAS (about 4-5 cycles of the HTW). Moreover, this time does not depend appreciably on the specific density ratio considered (which leads to the remarkable conclusion that once ‘ideal’ conditions for the emergence of the S-PAS are established in terms of competition between Stokes drag and buoyancy, a sort of universal behaviour is attained, by which the formation time is essentially minimised).
The specific role played by the different forces acting on particles will be clarified in the next section on the basis of dedicated numerical experiments. As a concluding remark for this section, we wish to highlight that we carried out other simulations in conditions somehow ‘specular’ to those reported in figure 11 by considering the reciprocal of the density ratio while keeping 
$St$ unvaried. Unfortunately, for the several extra cases considered, no PAS emerged. This result however does not allow us to conclude that S-PAS does not form for systems with particles denser than the fluid. We will come back to this important concept later (see § 5). We wish also to point out that, in the light of the considerations elaborated in § 2.2, the upper boundary of the existence region shown in figure 11 might be located in a range of values of 
$St$ for which the approximation of one-way coupling is no longer applicable (this specific point would require further investigation in the framework of a more sophisticated coupling approach, a cue for a future dedicated analysis).
3.3.2. Role of different forces on the process of PAS formation
To put our results in a broader perspective and provide additional insights into the underlying physics, a series of ‘ad hoc’ numerical experiments have been performed where fluid flow conditions have been kept unchanged while addressing the role of the added mass force (by tuning 
$\xi$ and 
$Fr$ while maintaining unaltered the Stokes number). As an additional step of this strategy, for some cases, we have also removed gravity from the set of forces potentially influencing particle motion (in other words, in order to obtain a robust understanding a peculiar analysis hierarchy has been used where various processes have been enabled or disabled in order discern the influence of effects which in the real world would be inextricably intertwined).
In this regard, attention has been paid to the non-dimensional parameters directly appearing in the Maxey–Riley equation (2.7). The dimensionless group appearing in front of the buoyant term is 
$(\xi -1)Fr^{2}$. It can therefore be concluded that the Froude number and density ratio can be adjusted in such a way that this product remains constant. In this way the buoyant force acting on the generic particle would remain exactly the same (though inertia does not, as witnessed by the factor 
$\xi +1/2$ in (2.7).
In order to maintain constant also the Grashof number, while changing g we have adjusted the thermal expansion of the fluid accordingly. None of these manipulations affects 
$St$.
To fix the ideas about the application of this modus operandi, in the following we will indicate with the suffix ‘
$0$’ the parameters relating to the simulations described in § 3.2, while no suffix will be used for the new parameters. By introducing a scale factor, 
$k$, such that the new acceleration due to gravity becomes 
, in order to maintain unaltered the fluid flow conditions the new coefficient of thermal expansion has to be multiplied by the same factor, i.e. 
$\beta =k\beta _0$ (enabling flow similitude, i.e. 
$Gr=Gr_0$). With these choices the new Froude number becomes 
$Fr=\sqrt {k}Fr_0$.
Since our goal was to maintain unvaried the buoyant term in the particle tracking, we set 
$(\xi -1)Fr^{-2}=(\xi _0-1)Fr_0^{-2}$. Using this expression and taking into account the definition of 
$Fr$, the following relationship could be easily worked out, 
$\xi =k\xi _0-k+1$, allowing us to determine the value to assign to 
$k$ that corresponds to the desired value of the new density ratio.
In the framework of this procedure, we have considered two representative conditions among those reported in figure 11 for considerably different particle sizes, i.e. the cases 
$St\simeq 1.2\times 10^{-5}$, 
${\xi }_0=0.8$, and 
$St\simeq 1.7\times 10^{-4}$, 
${\xi }_0=0.987$. Then, we have set a value of 
$\xi =1\times 10^{-3}$ which on the basis of the previous considerations provided 
$k\simeq 5$ and 
$k\simeq 75$ for the two different cases, respectively. Subsequently, particles have been added to the liquid (their number was exactly the same as for the original case) and tracked until accumulation occurred. This numerical experiment has revealed that S-PAS are formed even in the presence of large density ratios (provided 
$\xi <1$) and their morphology is not significantly influenced by the ‘added mass’ effect, as qualitatively and quantitatively substantiated by figure 13.
Figure 13. Comparison between the S-PAS obtained for different density ratio and 
$St\simeq 1.2\times 10^{-5}$ according to the criteria discussed in § 3.3.2. The original flow conditions were 
${\xi }_0=0.8$ (a,c), while in the other case (b,d) 
${\xi }=1\times 10^{-3}$.
As anticipated at the beginning of this section, a further understanding of the observed dynamics has been gained by ‘switching off’ gravity in the particle tracking model (i.e. in this way we have formally considered the limiting case 
$Fr\to \infty$) while keeping all other flow conditions unvaried (also retaining the buoyancy effect in the Navier–Stokes equations). Simulations have been conducted for a subset of cases extracted from those reported in figure 11. Remarkably, none of these cases has shown any sign of structure (H-PAS or S-PAS) formation.
Some important conclusions can therefore be drawn on the basis of these numerical experiments. While particle buoyancy is clearly the key ingredient for the emergence of the S-PAS for the conditions considered in the present work, the added mass force does not contribute to the underlying mechanisms.
4. Discussion
As a natural continuation of the physical interpretations provided in § 3.3.2, the present section is dedicated to a ‘recapitulation’ of the main physical mechanisms supporting the two types of PAS examined in the present work and the related ‘differences’.
We wish to recall that the PAS reported in figures 5 and 6 were obtained for isodense particles taking into account their interaction with the free interface. For such a specific case, while the KAM torus plays the role of backbone for PAS formation, the effective tendency of isodense particles to be entrained (or repelled) by it can be put in relation with their size using existing criteria in the literature. The interaction with the free surface can force particles to change fluid streamline and be progressively transported into the KAM torus where their trajectories would become stable. As illustrated in § 3.2, however, if the Stokes number is large enough, mechanisms are enabled that can prevent some particles from being captured by the attractor. In the present manuscript we have used these criteria to explain the variable (non-monotone) level of particle scattering visible in figure 5 for different particle diameters. The varying degree of scattering is a consequence of the interplay of the interaction-with-the-free-interface effect (which supports the emergence of PAS and is strengthened for larger values of 
$St$) and the opposite role represented by the ‘repel’ condition of Sapsis & Haller (Reference Sapsis and Haller2008), i.e. (3.6), which is also enhanced for larger values of the particle Stokes number.
The cause-and-effect relationships driving particle accumulation in figure 9 are completely different. The repellence theorem by Sapsis & Haller (Reference Sapsis and Haller2008) is no longer applicable to this case as gravity is present and particles are not neutrally buoyant. Moreover, (intentionally) no particle-interface interaction has been considered. For this case, the degree of accumulation or scattering of the particles must essentially be ascribed to the particle buoyancy effect (which in turn depends on 
$St$ and the density ratio 
$\xi$). For a fixed simulation time, scattering is enhanced if, starting from the curve (solid line in figure 11) bounding the region of existence of S-PAS from above, 
$St$ is decreased. Another way to think about this trend is to consider (see figure 12) that, for a fixed 
$\xi$, the structure formation time increases as 
$St$ becomes smaller and can even become virtually infinite (too long to be measured by means of the numerical approach employed in the present work) if the particle Stokes number is further decreased.
To see how these dynamics work from a purely mathematical point of view, it is sufficient to consider that since for the S-PAS no particle-surface interaction model has been implemented, the phenomenon is governed by the Maxey–Riley equation only. In particular, the main mechanism allowing particles to change streamline in the fluid is represented mathematically by the third term on the right-hand side of (2.7), namely, 
$(\xi -1)/(Fr^{2}(\xi +1/2))$. By exerting a disturbing role on the particles (which would otherwise follow the streamlines of the fluid-dynamic field) this term allows particles to jump from one streamline to another.
In other words, for the conditions considered in § 3.3, the dynamics are essentially governed by the interplay between this buoyancy term and the drag (second term on the right-hand side of (2.7)); only if 
$\xi$ and 
$St$ pertain to the existence region shown in figure 11, an equilibrium state can be attained where these two terms balance each other. In such conditions, once particles are captured by the KAM torus, they remain inside it. If such equilibrium is not attained, the disturbing role played by the particle buoyant term is so strong that it can make the accumulation process unstable (forcing all particles to cluster in proximity to the top rod).
5. Conclusions
Starting from the realization that some relevant models and theories have been appearing over the last years to interpret the PAS phenomena, the first objective of this work was to pursue the investigation of some still poorly addressed or overlooked effects. For these reasons, a peculiar modelling and analysis hierarchy has been implemented in order to ‘filter out’ already known mechanisms and concentrate selectively on other phenomena or processes. This specific approach has been made possible by the numerical framework used to address the problem (potentially allowing the user to enable or disable specific forces or effects by acting on the governing equations and related boundary conditions).
Attention has been expressly dedicated to the problem of particle accumulation structures in a high aspect-ratio (
$AR=0.9$) liquid bridge under the influence of gravity.
The fluid flow conditions have been assumed to be similar to those considered by Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007) in previous on ground experiments for which a particle accumulations structure emerged in a nearly isodense system (
$\xi \simeq 1$).
The execution of the numerical simulations, based on an Eulerian–Lagrangian one-way coupling approach, has been articulated in two steps. First, the base (carrier) flow field has been determined transcending the presence of particles. The solution has been found to be in good agreement with the conditions determined experimentally. Ensuing simulations conducted considering isodense particles have revealed the ability of these particles to demix from the fluid and accumulate forming a pattern resembling that found by Schwabe et al. (Reference Schwabe, Mizev, Udhayasankar and Tanaka2007).
In particular, numerical experiments have been conducted for a set of different particle sizes. For the smaller particles case (1 micron in size), a well-defined structure has been obtained, with features in accordance with those observed experimentally. For a larger particle case (particle diameter 20 microns), a less defined (but still recognizable) structure has been found with several particles being scattered within the whole LB domain. An explanation for such a trend has been elaborated on the basis of the criterion by Sapsis & Haller (Reference Sapsis and Haller2008). Simulations executed for larger values of the particle Stokes number, however, have revealed that the competition between the beneficial interaction-with-the-free-interface effect (enhanced for larger values of 
$St$) and the disturbing role mathematically represented by the condition of Sapsis & Haller (Reference Sapsis and Haller2008) can give rise to non-monotone behaviours.
Perhaps most important of all has been the identification of new routes to as-yet unknown phenomena, which are enabled when non-isodense particles are considered. A series of numerical experiments carried out for particles of various dimensions and densities have revealed that a sub-region of the space of parameters exists supporting conditions favourable for the formation of another pattern of accumulation.
The structures emerging within such a region of existence have been found to be consistent with a completely different (almost unknown to the scientific community) type of PAS, previously observed only by Sasaki et al. (Reference Sasaki, Tanaka and Kawamura2004) (for which we have coined the definition of S-PAS). This structure differs from those reported in earlier work on the subject due to its peculiar morphology that can hardly be categorised on the basis of other geometrical (topological) paradigms for PAS.
Dedicated numerical experiments have been conducted accordingly with the intention to elucidate the role played in the formation of this alternate structure by various forces. Such analysis has revealed that, for the conditions considered in the present work, the formation of the S-PAS must primarily be ascribed to the influence of Stokes drag and buoyancy when particles lighter than the carrier liquid are considered.
The main outcome of the present study is that particle buoyancy must be considered as a third possible root cause (in addition to those already known in the literature relying on purely inertial or particle-interface interaction effects) for the emergence of PAS. This additional mechanism can make different types of PAS (e.g. H-PAS and S-PAS) multiple solutions to the same problem, their basin of attraction differing with regard to the influence exerted by gravity on particles (to elucidate further the significance of this observation, one should keep in mind that these structures might manifest at the same time in the same physical space if isodense and non-isodense particles were initially seeded in the fluid).
Despite the tendency of the system to develop solutions that may coexist in space, surprisingly, no structure of the S-PAS type have been found numerically for conditions corresponding to particles heavier than the surrounding fluid (
$\xi >1$), which is a relatively counter-intuitive finding if one considers that many experiments have been conducted using particles of such a kind. In light of the mathematical and physical arguments reported in the present work, however, a natural justification for this apparent idiosyncracy would follow naturally simply considering that (1) the number of influential parameters governing the emergence of PAS is relatively large (and, as discussed in § 4, the present study also ‘adds’ gravity and buoyancy to them), and (2) what sets the specific case 
$m=1$ apart from the others is the very elusive nature of particle structures for this specific wavenumber. In line with these considerations, we argue that the present findings should not be misread as implying that the S-PAS does not exist for 
$\xi >1$; rather its existence should be sought in a different sub-region of the space of parameters (presumably corresponding to different values of the Grashof number, which in the present work was fixed to 
$Gr\simeq 6430$). The Grashof number can indeed have an effect on the PAS by modifying the value of the critical Marangoni number and, therefore, the interval of 
$Ma$ in which PAS can emerge. Future efforts shall therefore be expressly devoted to clarify the role of this parameter on PAS formation through adequate parametric exploration of the five-dimension (
$\xi$, 
$St$, 
$Ma$, 
$Gr$, 
$Fr$) space.
Acknowledgements
This work has been supported by the UK Engineering and Physical Sciences Research Council (EPSRC grant EP/R043167/1) and the UK Space Agency in the framework of the JEREMI (Japanese European Research Experiments on Marangoni Instabilities) project sponsored by ESA and JAXA. The authors would like to express their gratitude to Professor D. Schwabe for his valuable suggestions and indications. The authors wish also to thank Dr A. Passalacqua from Iowa State University for providing them with the subroutines required for eigenvalue computation in OpenFOAM.
Declaration of interests
The authors report no conflict of interest.
