Numerical investigation of the dynamics of a rigid spherical particle in a vortical cross-slot flow at moderate inertia

The study of flow and particle dynamics in microfluidic cross-slot channels is of high relevance for lab-on-a-chip applications. In this work, we investigate the dynamics of a rigid spherical particle in a cross-slot junction for a channel height-to-width ratio of 0.6 and at a Reynolds number of 120 for which a steady vortex exists in the junction area. Using an in-house immersed-boundary-lattice-Boltzmann code, we analyse the effect of the entry position of the particle in the junction and the particle size on the dynamics and trajectory shape of the particle. We find that the dynamics of the particle depend strongly on its lateral entry position in the junction and weakly on its vertical entry position; particles that enter close to the centre show trajectory oscillations. Larger particles have longer residence times in the junction and tend to oscillate less due to their confinement. Our work contributes to the understanding of particle dynamics in intersecting flows and enables the design of optimised geometries for cytometry and particle manipulation.


Introduction
The deformability of blood cells is a valuable label-free biomarker to indicate changes in the cell state 1 . Changes in the nucleus or cytoskeleton of cells are associated with various diseases, including malaria 2,3 , sickle cell disease 4 , cancer 5 and sepsis 6 . For example, oral cancer cells have been found to be more deformable than healthy oral epithelial cells 7 , thus facilitating metastasis 8 . Leucocytes have been shown to become stiffer or slower, passing through pores upon stimulation with chemokines associated with infection, such as sepsis 9 . Conventional methods for single-cell analysis, for instance, atomic force microscopy 10 , optical stretching 11 and micropipette aspiration 12 , are useful for fundamental research in mechanobiology, but their low throughput limits the potential for clinical applications.
Fast and accurate analysis of blood samples is important for identifying diseased cells, thus accelerating clinical decision-making 13 . Inertial microfluidics (IMF) has been widely used for cell manipulation, separation, and characterisation in biomedical applications 14 . IMF operates in a flow regime where inertial forces influence the fluid and particle dynamics 15 . In the 1960s, Segré and Silberberg 16,17 observed radial migration of rigid spherical particles towards an annulus at around 60% of the radius of a tube when the Reynolds number is of the order 10-100. Since then, several studies investigated the nature of the inertial lift forces in straight channels [18][19][20] . In the 2000s, it was recognised that the Segré-Silberberg effect could be exploited in a microfluidic setting 21 and in more complex channel geometries, such as curved 22,23 , serpentine 24,25 and spiral channels [26][27][28] .
Since its inception, IMF has played an important role in image-based deformability cytometry (DC). DC uses different channel geometries, such as linear channels 29 , constrictions 30 and cross-slot channels 31 . While traditional cytometry techniques, such as atomic force microscopy, micropipette aspiration and optical stretchers, are limited to low cell processing rates (10-100 cells/h) 29 , DC yields significantly higher throughput (around 10-1000 cells/s) 32 , making it more suitable for diagnostic applications where large quantities of cells should be measured to have reliable measurements.
Cross-slot channels ( Fig. 1) are particularly appealing for DC since the flow field in their junction features both high stresses, therefore deforming cells, and a relatively low flow speed, thus facilitating optical imaging. The fluid flow in cross-slot junctions is complex and has been extensively investigated. For specific ranges of channel aspect ratio and Reynolds number, a spiral vortex can appear in the junction 33 , which has been explained by a symmetry-breaking bifurcation above a critical Reynolds number for a given channel aspect ratio 34 . Furthermore, Burshtein et al. 35 determined the characteristics of a steady-state vortex in cross-slot flows at higher Reynolds numbers, such as the vortex core structure and time-periodic fluctuations around the stagnation point in the junction, by tuning the aspect ratio of the device. The design of the junction has been optimised to achieve more homogeneous strain rates 36,37 , while Zhang et al. 38 identified different flow regimes in the cross slot by changing the flow conditions. Recently, a crossslot device has aided the early detection of sepsis by using the deformability of leucocytes as a biomarker 39 .
Despite earlier works focused on fluid dynamics and particle trapping [40][41][42] in cross-slot junctions, only a limited amount of studies have considered moderate inertia 43,44 but focused mainly on cell deformation and not on the transient behaviour. The interaction of a particle and the vortex at moderate inertia has not been thoroughly investigated. It is currently unclear how flow field characteristics and particle

Definition of low-velocity zones A and B
x z c d Fig. 1 a 3D illustration of the cross-slot geometry with two inlets and two outlets. b Definition of the initial position (x 0 , z 0 ) of the particle's centre of mass on the inlet cross-section. The origin (red dot) is located on the inlet's centre line, and the particle initially moves along the y-axis (pointing into the plane of view). c Top view of streamlines showing the spiral vortex at Re = 120 and channel aspect ratio α = 0.60, as defined in Section "Physical model, parameters and non-dimensional groups". The colour bar indicates the velocity magnitude (normalised by the average inlet velocity). d Definition of low-velocity zones A and B in the junction, which are important for the particle dynamics. View along the outlet axis properties correlate with observable metrics, such as the shape of the particle trajectory and the residence time in the junction. In this paper, we use immersed-boundary-lattice-Boltzmann (LB) simulations to explore the dynamics of a spherical, rigid particle in an incompressible Newtonian liquid in the steady-vortex regime (Section "Physical model and numerical methods"). We pay particular attention to the characterisation of the geometry, observables and simulation initialisation (Section "Cross-slot geometry and simulation set-up"). We investigate the effect of initial particle position and size on the trajectory and residence time of the particle in the cross-slot junction at a Reynolds number of 120 and a channel aspect ratio of 0.6 (Section "Results and discussion"). Our results show that particle trajectory and residence time strongly depend on the initial lateral position. Under the investigated flow conditions and channel geometry, the vortex affects larger particles less compared to smaller ones. In Section "Conclusions", we provide a concluding discussion. We hope that our results stimulate further work, for example, the investigation of deformable particles in the cross-slot junction.

Cross-slot geometry and simulation set-up
The geometry consists of a cross-slot channel, as shown in Fig. 1a. The origin of the Cartesian coordinate system is at the centre of the cross-slot geometry. The flow in the inlets is along the y-axis (inlet axis), the flow in the outlets is along the x-axis (outlet axis), and the z-axis denotes the height direction (height axis). Depending on the combination of Reynolds number, height-to-width aspect ratio and the ratio of the outlet flow rates, different flow regimes inside the cross-slot junction can be observed, for example, a vortex which might be steady or unsteady 35,38,45 . In this work, we are interested in the steady-vortex regime since the steady vortex leads to reproducible and non-trivial particle behaviour, which is relevant for flow cytometry applications. Therefore, we set the Reynolds number to Re = 120, the channel aspect ratio to α = 0.6, and the ratio of both outlet flow rates to 1, which leads to the formation of a steady vortex 34 .
Since both inlet flow rates and both outlet flow rates are identical, the resulting flow field, including the vortex, is symmetric with respect to the inlet-height plane (yz-plane) and under a 180°-rotation about the outlet axis (x-axis) in the absence of the particle. Therefore, it is sufficient to enter the particle only on one side of the inlet-height plane, and all initial particle positions are x 0 > 0. Furthermore, we only inject particles at one inlet (y 0 < 0). Due to the presence of the vortex spiralling about the outlet axis, we need to explore the entire height range of initial particle positions. This procedure allows us to scan all physically unique initial particle positions without running multiple simulations that are physically equivalent.
The height H is constant within the entire domain, and both inlets and both outlets have the same width W. We use 36 grid points for H and 60 grid points for W, giving 45 grid points for the hydraulic diameter D h . The length of the inlets and outlets is L in = L out = 4D h , which we found to be sufficient to give results independent of the inlet and outlet lengths. The relaxation time is chosen as τ = 0.55Δt in all simulations, resulting in a viscosity of ν = 1/60Δx 2 /Δt. The peak and average values of the inlet velocity profile are u max ¼ 0:068 Δx=Δt and u ¼ 0:033 Δx=Δt, respectively. We consider particle sizes resulting in confinement values χ in the range [0.35, 0.65]. The particle is neutrally buoyant in all cases, ρ p /ρ f = 1.
To reduce the required time for the vortex to form, we introduce a force field that accelerates the fluid in the junction's upper and lower halves in opposite directions along the inlet axis for a finite amount of time. Once the force has been switched off, the vortex flow is allowed to converge to a steady state with a convergence criterion of 10 −6 for the relative difference of the flow field between consecutive time steps. It is reported that the right or left-handed orientation of the spiral vortex occurs with the same probability 34 . For consistency, we initialise the simulations in such a way that the vortex is always left-handed, as shown in Fig. 1d. Eventually, the particle is injected at the desired initial position next to the inlet at y 0 < 0.
We always initialise particles with a finite value of x 0 since a particle initially on the mid-plane between the outlets would be caught in the vortex indefinitely in the absence of distortions. Since numerical noise in the simulation would eventually let the particle escape in one of the outlets, the residence time would be determined by numerical noise and not deterministic physics. In real-world experiments, channel and particle imperfections always lead to finite residence times in the geometry considered here.
Due to inertial forces, particles initialised near the inlet would migrate away from their prescribed lateral coordinates (x 0 , z 0 ) while moving along the inlet, therefore losing control over the location where the particle enters the actual junction. To avoid this premature particle migration in the inlet, we lock the x-and z-components of the particle position while allowing the particle to rotate and move along the inlet axis freely. The 'guidance' of the particle is switched off when the particle reaches a point located D h upstream of the junction; the particle is allowed to move freely afterwards. Our simulations show that, at an upstream distance of about one hydraulic diameter from the junction, the velocity profile in the inlet starts to be affected by the junction. Therefore, the initial cross-sectional position (x 0 , z 0 ) defines at which point the particle enters the region of influence of the junction.
In the interest of using a concise notation, we definê Therefore, the accessible range of initial positions for a point particle in the inlet channel isx 0 2 ½À1; þ1 and z 0 2 ½À1; þ1.

Results and discussion
We present and discuss our simulation results in two main sections. Section "Effect of initial particle position" covers the effect of initial particle position (x 0 ,ẑ 0 ) on the trajectory and residence time of the particle. Section "Effect of particle confinement" focuses on the effect of particle size, including a neutrally buoyant point particle.

Effect of initial particle position
To study the effect of the initial particle position, we investigated a range of initial positions (x 0 ;ẑ 0 ), as summarised in Table 1. A first simulation study revealed that the richest particle dynamics could be found for small values of x 0 . Hence, we ran more subsequent simulations for smaller values ofx 0 . For the values ofẑ 0 , we chose a starting position on the mid-plane (ẑ 0 ¼ 0) and the equilibrium positions of differently sized particles (χ between 0.35 and 0.75) in a straight duct with the same aspect ratio (α = 0.60) as the inlet of the cross-slot geometry at Re = 120, see App. Equilibrium positions in a straight channel are equivalent to the inlet of the cross-slot. The same confinement values are used in Section "Effect of particle confinement" to investigate the effect of particle size. The particle is neutrally buoyant and has a fixed confinement of χ = 0.4. The channel aspect ratio is α = 0.6, and the Reynolds number is 120. Figure 2a, b shows some representative particle trajectories inside the cross-slot junction for two different values ofx 0 and a range of values ofẑ 0 . The trajectories of rigid particles found here are qualitatively similar to those presented by Gosset et al. 31 in their supporting videos and Zhang et al. 46 .

Particle trajectories
It is inconvenient to show the particle trajectories for all the simulations. Instead, we extract key characteristics from the particle trajectories that are also experimentally relevant. Figure 2c illustrates the chosen 'macroscopic' metrics of the inlet-outlet projection (xy-projection, 'top view') of the particle trajectories: the number N of peaks observed in the central area of the cross-slot junction, the amplitude A of each peak, and the lateral distance L between consecutive peaks. A peak is defined as a minimum or maximum of the inlet-outlet projection of the particle trajectory. The residence time T res is the length of time the centre of mass of the particle spends in the central junction area,x 2 ½À1; þ1 andŷ 2 ½À1; þ1.
Particles initialised very close to the mid-plane of the inlet, x 0 ¼ 0:0033 in Fig. 2a, give rise to damped oscillations in the Table 1 Non-dimensional initial particle positions (x 0 ,ẑ 0 ) used in Section "Effect of initial particle position" Definition of a peak (N ), the peak amplitude (A ) and the axial inter-peak distance (L).

Fig. 2
Particle trajectories projected onto the inlet-outlet plane (top view). The channel aspect ratio is α = 0.6, the Reynolds number is Re = 120, and the particle confinement is χ = 0.4. Different curves show different values ofẑ 0 for ax 0 ¼ 0:0033 and bx 0 ¼ 0:15. Solid and dashed lines refer to negative and positive values ofẑ 0 , respectively. c Illustration of a typical particle trajectory. Peaks (blue dots) are extrema of the curve,Â i denotes the ith amplitude (normalised by W/2) of a peak with respect to the centre-line of the outlet (x-axis), andL i is the axial distance (normalised by W/2) between the (i − 1)st and ith peaks (withL 1 playing a special role, denoting the distance between the first peak and the centre-line of the inlet). Note that the aspect ratio of the plots has been changed to improve readability trajectories. For moderate distances from the mid-plane, Fig. 2b, the oscillations are either weak or completely absent. Figure 2 also shows thatẑ 0 has a stronger effect on the shape of the trajectories whenx 0 is small. We found that particles only show appreciable oscillations when the initial positionx 0 is smaller than 0.1. We will discuss these findings in the remainder of this section by providing further 'microscopic' information: instantaneous particle velocities (translational and angular) and forces (lift and drag) which we use to rationalise the behaviour of the macroscopic observables.
In the following, we will focus on an analysis of the macroscopic metrics defined in Fig. 2c (number of peaks, peak amplitude, inter-peak distance) and the residence time to characterise the particle trajectories and their dependence on initial positions. Figure 3a shows the number of peaks, N, for all investigated initial positions (x 0 ;ẑ 0 ). Only for the narrow rangex 0 <0:05, some trajectories have two or more peaks. Increasingx 0 from 0.0033 to 0.01 reduces N by one for all the investigated values ofẑ 0 . Upon a further increase ofx 0 , the number of peaks reduces to one or zero, depending on z 0 . No peaks are observed for any trajectory withx 0 ! 0:25.

Number of peaks of particle trajectories
The initial height of the particle,ẑ 0 , has a less pronounced effect on N. Trajectories with an initial position within the range À0:22 ẑ 0 þ0:22 have the same number of peaks for identical values of jẑ 0 j although the trajectories look different (cf. Fig. 2). For jẑ 0 j>0:22, trajectories withẑ 0 >0 tend to have fewer peaks than those withẑ 0 <0 for the same magnitude ofẑ 0 . For the entire range ofx 0 , particles initially located atẑ 0 ¼ 0:5 follow trajectories without any peaks.

The peak amplitude of particle trajectories
The normalised peak amplitude,Â ¼ A=ðW =2Þ, gives finer information about the shape of the trajectories. The data in Fig. 3b show the general trend that, for trajectories with more than a single peak, the amplitude of consecutive peaks is smaller than the amplitude of the  the same trend has also been found for other values ofx 0 (data not shown). c Amplitudes of the first peak, A 1 , for all investigated initial positions, as long as the peak exists. d, e Normalised inter-peak distance,L, of trajectories as defined in Fig. 2c for all investigatedẑ 0 for dx 0 ¼ 0:0033 and ex 0 ¼ 0:01. f Normalised inter-peak distance between the first and second peak,L 2 , forx 0 ¼ 0:0033 and 0.01. Lines are guides for the eyes previous peak, indicating that the spiralling motion of the particle is damped towards the outlet. The peak amplitudes are larger when particles are released atẑ 0 ! 0, compared to trajectories withẑ 0 <0 highlighting the influence of the vortex handedness on the particle's motion in the junction under the investigated conditions. Figure 3c shows the behaviour of the normalised amplitude of the first peak,Â 1 , for all studied initial positions. For a given value ofẑ 0 , there is little change of A 1 withx 0 , evidenced by the collapse of data for different values ofx 0 , as long as the first peak exists. Changingx 0 merely affects the number of peaks, but not the amplitude of the first peak. The amplitudeÂ 1 increases withẑ 0 for negative and positive values ofẑ 0 , respectively, while the amplitude forẑ 0 ¼ 0 escapes the trend and defines a local maximum.

Axial inter-peak distance
Like the amplitude of the peaks,Â, the axial distance between consecutive peaks,L ¼ L=ðW =2Þ, serves as a quantitative measure of the particle trajectories. We defineL i as the normalised axial distance (along the outlet axis) between the peaks i − 1 and i, as long as both peaks exist. The axial distanceL 1 denotes the distance of the first peak from the inlet centre-line as illustrated in Fig. 2c. Figure 3d, e shows results forx 0 ¼ 0:0033 and 0.01, respectively. In both panels, only those curves are shown for which at least two peaks exist such thatL 2 is defined. Since trajectories withẑ 0 >0 tend to have fewer peaks than those withẑ 0 <0, there are fewer data points for the former. In all observed trajectories, the inter-peak distance grows between consecutive pairs of peaks, i.e.L iþ1 >L i . Furthermore,L i for a fixed i increases with jẑ 0 j: the farther away from the horizontal mid-plane a particle is located initially, the longer its spiralling trajectory is stretched toward the outlet. The data also show thatL is larger for positive z 0 than for negative z 0 with the same magnitude jẑ 0 j. Comparing data from Fig. 3d, e, we see thatL increases withx 0 : particles initially farther away from the vertical mid-plane have trajectories with larger inter-peak distances. Specifically, Fig. 3(f) shows the increase ofL 2 withx 0 .

Residence time
The residence time of a particle in the junction (cf. Section "Particle trajectories") is an important parameter to consider. Due to the symmetry of the problem, a particle withx 0 ¼ 0 would in principle remain in the junction forever, T res ! 1. For any finite value ofx 0 , the residence time should be finite since the particle is located closer to one outlet and can leave the junction area more easily. In fact, Fig. 4a shows that T res sharply decreases withx 0 for smallx 0 and drops further whenx 0 is increased.
Although T res is largely determined byx 0 , there is also a clear dependency onẑ 0 (Fig. 4b). For all studied values ofx 0 , there is a mild increase of T res withẑ 0 : the residence time is shorter when the particle is initialised in the lower channel half (ẑ 0 <0) compared to initial positions in the upper channel half.
In order to understand the behaviour of the residence time, we first analyse the velocity and acceleration of the particle along the outlet axis (x-axis) inside the junction (Fig.  4c, d). Particles that are initially close to the vertical midplane (smallx 0 ) accelerate more slowly than a particle that enters the junction farther away (largerx 0 ). The slow acceleration of a particle with smallx 0 is caused by the particle covering nearly identical portions of volume on either side of the mid-plane; hence the net force accelerating the particle towards one of the outlets is nearly zero. Once the particle is slowly moving towards one outlet, it gets increasingly exposed to the drag force exerted by the flow in that outlet channel, and the particle acceleration, and therefore velocity, increases. A particle that enters the junction far away from the mid-plane, on the one hand, is closer to the outlet and, on the other hand, experiences a larger drag force early on and leaves the junction sooner; hence, the residence time is shorter. The residence time also depends onẑ 0 (Fig. 4b), but to a lesser extent than it depends onx 0 . This dependency stems from the particle acceleration _ u x being smaller for particles entering the junction above the horizontal mid-plane (ẑ 0 >0). The mechanism behind this effect will become obvious at the end of this section. Interestingly, the acceleration increases exponentially with time, as revealed in Fig. 4d.
The behaviour of the number of peaks, N, and the interpeak distance,L, is determined by an interplay of the residence time and the ability of the particle to follow the vortex and revolve about the outlet axis (x-axis). The slower the particle moves towards the outlet and the faster it revolves about the outlet axis, the larger the number of peaks and the smaller the inter-peak distance. Figure 5a depicts the distribution of the vorticity magnitude of the unperturbed flow field in the junction on the horizontal mid-plane (ẑ ¼ 0). The vortex can be recognised by a region of high vorticity stretched along the outlet axis. The vortex causes the particle to revolve about the outlet axis with an angular velocity Ω x , thus creating the spiralling trajectories illustrated in Fig. 2a. Figure 5b shows the angular velocity, averaged over the time the particle is located inside the junction area, Ω x . The number of peaks N, residence time T res and the average angular velocity Ω x are related according to remembering that a full revolution of 2π radians leads to two peaks as defined in Fig. 2c. Particles initially located closer to the vertical mid-plane (smallx 0 ) experience a higher angular velocity than particles entering the junction farther away from the mid-plane (largex 0 ). Combining the faster revolution Ω x with the smaller linear acceleration _ u x (and therefore larger residence time) for particles withx 0 ! 0 explains why these particles have the largest number of peaks, N, according to Eq. (2). Since peaks are only counted in the junction area, which has a fixed size, a trajectory with more peaks has smaller inter-peak distancesL, which explains the trend ofL becoming smaller forx 0 ! 0. Looking at particles with smallx 0 , Fig. 5b reveals that Ω x is significantly smaller for particles withẑ 0 >0, while the residence time in Fig. 4(b) shows only a moderate increase for these particles, compared to particles withẑ 0 <0. Consequently, particles withẑ 0 >0 show fewer peaks than particles withẑ 0 <0 as long asx 0 is small. The smaller revolution frequency of particles withẑ 0 >0 can be understood as follows: Fig. 5c shows that, upon entering the junction, these particles move over the vortex, rather than getting trapped in the vortex, then reach a region with relatively small vorticity before finally entering the vortex. Averaged over the time they spend in the junction, these particles, therefore, do not revolve as much as the particles that are able to enter the vortex immediately.
In order to understand the variation of the peak amplitude with the initial positionsx 0 andẑ 0 , we first inspect the particle trajectories as seen from the outlet (projection on the inlet-height plane, Fig. 5c). A particle that is released atẑ 0 <0 experiences an immediate upward lift force and is drawn into the vortex where the particle revolves with small amplitude. Contrarily, a particle with largerẑ 0 features a flatter trajectory  Fig. 4 (a,b) Residence time, T res (as defined in Section "Particle trajectories"), as a function of ax 0 and bẑ 0 . The residence time is normalised by the advection time, Eq. (7). Time evolution of c particle velocity and d logarithm of particle acceleration inside the junction along the outlet axis (x-axis) tends to avoid the vortex and approaches low-velocity zone B. Only after the encounter with zone B the particle is pulled into the vortex, resulting in a larger peak amplitudeÂ 1 . None of the particles simulated here is able to reach zone A; this would only be possible for particles with very small confinement χ released atẑ 0 ! À1. The different shapes of particle trajectories can be understood by inspecting the streamlines in Fig. 1d: streamlines with slightly negativeẑ 0 turn upwards and reach the vortex immediately, while streamlines with slightly positiveẑ 0 tend to avoid the vortex and approach zone B; only streamlines near the bottom of the channel are drawn towards zone A. The independence of A 1 fromx 0 can be explained by the shape of the highvorticity region in Fig. 5a. Since the vorticity profile is largely constant along the outlet axis inside the entire junction area, the initial interaction between the particle and the vortex is not strongly sensitive tox 0 . Finally, Fig. 5c explains the trend of the residence time to increase withẑ 0 (Fig. 4b). Particles with largerẑ 0 spend some time in the low-velocity zone B, which delays the particle's acceleration towards the outlet, therefore increasing the residence time in the junction.

Effect of particle confinement
In this section, we analyse the effect of the size of the particle on its dynamics in the cross-slot junction. We investigate five different confinement values, χ ∈ [0.35, 0.65] and the limiting case of a mass-less point particle, χ → 0. In order to keep the number of free parameters manageable, we focus onx 0 ¼ 0:0033, which leads to the largest number of peaks andx 0 ¼ 0:15, which is the lowest value ofx 0 for which nearly no peaks have been observed in Fig. 3(a). Additionally, we fix the initial positions atẑ 0 ¼ ± 0:31 which is the lateral equilibrium position of a particle with χ = 0.65 in a straight channel with the same dimensions and Reynolds number as the inlet channel of the cross-slot geometry. The hydrodynamic conditions and channel aspect ratio are kept constant at Re = 120 and α = 0.6, respectively. Initial position atx 0 ¼ 0:0033 In Section "Effect of initial particle position", we found that the vortex has a more profound influence on the particle dynamics whenx 0 ¼ 0:0033 than for the other investigated initial positions. Figure 6 shows the particle trajectories in the junction for the investigated confinement values forẑ 0 ¼ À0:31 in (a, e) andẑ 0 ¼ 0:31 in (b, f). As the confinement χ increases, the characteristic spiral trajectories seen for particles with low confinement are less pronounced or no longer observed at all. Forẑ 0 ¼ À0:31, the particle with χ = 0.65 leaves the junction without any oscillation. Forẑ 0 ¼ 0:31, no particle with χ ≥0.5 oscillates in the junction. A mass-less point particle, which follows the fluid streamlines, shows the most oscillatory behaviour of all particles forẑ 0 ¼ À0:31; contrarily, when initialised atẑ 0 ¼ 0:31, the point particle avoids the vortex and hardly oscillates.
The revolution frequency of the particle about the outlet axis, Ω x , provides a complementary view. When particles enter the junction below the horizontal mid-plane,ẑ 0 ¼ À0:31 (Fig. 7a), Ω x is highest for particles with χ < 0.5, including the point particle, and Ω x decays quickly for larger confinement. Forẑ 0 ¼ 0:31 (Fig. 7b), particles with 0.35 ≤ χ < 0.5 show a behaviour similar to that forẑ 0 ¼ À0:31. In contrast, both the point particle and the particles with χ ≥ 0.5 hardly oscillate when released above the horizontal mid-plane,ẑ 0 ¼ 0:31. While the oscillation of the point particle is suppressed by the particle reaching the lowvelocity zone B, larger particles are strongly confined when they approach the top wall of the junction and, hence, are unable to revolve freely. Figure 7c displays the residence time, T res , for all investigated cases. For particles that are released below the horizontal mid-plane,ẑ 0 ¼ À0:31, confinement has a  Fig. 1d) negligible effect on the residence time: for all particles, including the point particle, T res % 2t ad . The insensitivity of the residence time to the confinement can be explained by the observation that all particles remain in the vortical zone throughout their passage (cf. Fig. 6e), where particles are efficiently advected towards the outlet. The situation is different forẑ 0 ¼ 0:31 where the residence time strongly depends on particle size. While particles with 0.35 ≤ χ ≤ 0.5 leave the junction nearly as quickly as their respective counterparts atẑ 0 ¼ À0:31, both the point particle and the larger particles (χ > 0.5) require more time to leave the junction area. The increased residence time of the point particle is caused by the point particle reaching and remaining in the low-velocity zone B, as seen in Fig. 6f. The large particles with χ > 0.5 do not reach zone B; instead, they are nearly in contact with the top wall, which reduces their ability to follow the flow around the vortex or towards the outlet.
Initial position atx 0 ¼ 0:15 In Section "Effect of initial particle position", we demonstrated that the vortex has little effect on the motion of the particle when the particle enters the junction farther away from the vertical mid-plane, for example Here, we show that, even if particle confinement is varied, the vortex continues to have little influence on the motion of the particle atx 0 ¼ 0:15. Figure 6 depicts the particle trajectories in the junction for the investigated confinement values forẑ 0 ¼ À0:31 in (c, g) andẑ 0 ¼ 0:31 in (d, h). For bothẑ 0 ¼ À0:31 and 0.31, the particle experiences a strong acceleration towards the outlet once it enters the junction area since most of the particle volume is located in the outlet-facing half of the junction. During its passage, the particle does not reach the vortex core and, therefore, does not undergo a strong revolution about the outlet axis (Ω x not shown forx 0 ¼ 0:15). Consequently, in nearly all situations, the particle leaves the junction without displaying a peak in its trajectory. The only exception is the point particle for z 0 ¼ À0:31 with a single peak, although the point particle essentially follows a similar trajectory as the finite-size particles.
The residence time T res is consistently lower forx 0 ¼ 0:15 compared tox 0 ¼ 0:0033 (Fig. 7c). The shorter residence time is caused by the particle withx 0 ¼ 0:15 being accelerated towards the outlet earlier and, at least to some extent, the shorter distance to travel to the outlet. The residence time is essentially independent of confinement χ, and it is nearly the same forẑ 0 ¼ À0:31 andẑ 0 ¼ 0:31: all particles follow roughly the same trajectory and therefore experience a similar acceleration along the outlet axis (Fig. 6).

Physical model and numerical methods
We briefly outline the physical model (Section "Physical model, parameters and non-dimensional groups") and the numerical methods (Section "Numerical methods") used in this work.

Physical model, parameters and non-dimensional groups
We consider a single rigid, spherical particle immersed in an incompressible, isothermal, viscous Newtonian liquid. The incompressible Navier-Stokes equations and Newton's laws of motion determine the dynamics of the liquid and the particle, respectively. The no-slip condition is assumed at all solid-liquid boundaries.
The relevant parameters characterising the fluid flow are the liquid density ρ f , the kinematic viscosity ν, and the mean velocity u in the inlet of the cross-slot geometry. The particle is fully characterised by its radius a and density ρ p . As will be described in more detail in Section "Cross-slot geometry and simulation set-up", the geometry is fully determined by the width W and height H of the inlet and outlet channels of the crossslot geometry.
Following 34 , we define the channel Reynolds number as The non-dimensional confinement of the particle in the cross-slot geometry is The aspect ratio of the cross-section of the inlet and outlet channels is given by Note that we chose α < 1 in our study, justifying the use of H in the definition of the confinement, Eq. (4). The hydraulic diameter is important for the characterisation of the inflow behaviour of the particle. We further define the advection time scale to non-dimensionalise time ast ¼ t=t ad .

Numerical methods
We employ the lattice-Boltzmann (LB) method with the D3Q19 lattice 47 , the BGK collision operator 48 and the Guo forcing scheme 49 to solve the Navier-Stokes equations. The relaxation time τ of the BGK collision operator is linked to the kinematic viscosity according to where c s ¼ ffiffiffiffiffiffiffi ffi 1=3 p Δx=Δt is the lattice speed of sound, Δx is the lattice spacing and Δt is the time step.
The multi-direct-forcing immersed-boundary method with a three-point stencil is used to satisfy the no-slip boundary condition on the surface of the moving rigid particle 50,51 . The half-way bounce-back scheme 52 is employed to satisfy the no-slip condition on the stationary walls of the cross-slot geometry. Additionally, due to nonnegligible particle inertia, a coefficient of restitution is included to model perfectly elastic collisions to avoid any overlap between the particle and the wall.
We impose a fully developed rectangular Poiseuille flow with the desired mean velocity u at both inlets 53 . The boundary condition at both outlets is based on an extrapolation approach 54 modified in such a way that the average pressure on the outlet plane is constant in time. Therefore, the inlet velocity and the outlet pressure are imposed, while the outlets permit the swirling flow caused by the vortex at the centre of the cross-slot geometry to leave the domain without undesired reflections at the outlet planes.
To ensure that our simulation code accurately captures the fluid dynamics and the fluid-particle interaction, we compared the flow structure in the junction with available data 34 (App. Flow structure in the cross-slot junction) and the lateral equilibrium positions in a square cross-section with published results 55 (App. Test case: effect of initial position and confinement on lateral migration in the square duct).

Conclusions
We used the immersed-boundary-lattice-Boltzmann method to numerically investigate the dynamics of a rigid neutrally-buoyant spherical particle in a cross-slot junction for a channel height-to-width ratio of 0.6 and a Reynolds number of 120 for which a steady vortex exists in the junction. We characterised the effects of the initial position and the particle confinement χ = 2a/H (where a is the particle radius and H is the height of the channel) on the particle trajectories in the junction. Particles generally show a swirling motion within the cross-slot. Thus, observed trajectories were described in terms of macroscopic observables, such as the number of peaks of each trajectory (N), the peak amplitudes (A), the axial interpeak distances (L), and the residence time within the junction area (T res ). The simulations also provide access to microscopic observables, such as instantaneous linear and angular velocity and acceleration, which were used to rationalise the particle's behaviour.
There are several key observations when varying the location at which the particle with fixed size enters the junction:

•
The residence time T res strongly depends on the initial position of the particle. A particle entering the junction closer to the vertical mid-plane is accelerated more slowly towards its closest outlet, thus leading to a larger residence time. There is a mild effect of the initial height of the particle on the residence time. The behaviour of the residence time can be explained by the interaction of the particle with low-velocity areas in the flow field.
• Two competing mechanisms determine the number of peaks N and the axial inter-peak distance L seen in the particle trajectories: particles that are able to follow the vortex and revolve about the outlet axis more easily show a larger number of peaks and shorter inter-peak distances. Contrarily, particles with a shorter residence time have a trajectory with fewer peaks and larger interpeak distances. As a net effect, the number of peaks and inter-peak distance depend strongly on the lateral position and weakly on the height position of the particle when entering the junction.

•
The amplitude A of consecutive peaks of a given trajectory generally decrease as the particle follows the dissipating vortex along the outlet channels. The initial height of the particles has a strong effect on the amplitude: depending on the initial position, the particle either moves along the vortex, resulting in a larger amplitude, or it moves against the vortex, resulting in fast migration towards the centre of the junction and a smaller amplitude. The properties of the particle trajectories also depend on the particle confinement: • When entering the junction close to the vertical mid-plane, larger particles have a slightly higher residence time T res since a large fraction of the particle is initially located in the region drawn into the second outlet.

•
The trajectories of particles with small confinement show more peaks N with smaller inter-peak distance L because smaller particles are able to revolve faster about the outlet axis.
• Particles with larger confinement show a higher amplitude A of the first peak since these particles tend to pass the vortex first before being drawn into the vortex. The cross-slot geometry has been successfully employed in flow cytometry to deform cells 30,31,43 , where cell deformation is used as a bio-marker for the cell state 39 . In deformability cytometry using a junction, cells or particles are usually inertially focused before entering the junction 31 . Therefore, cells or particles are expected to enter and traverse the junction with typical trajectories, which can be quantified using the metrics used in the present study. The videos of Gossett et al. 31 display polystyrene beads, droplets and cells performing an oscillating trajectory in the cross-slot junction as seen from above (along the z-axis), similar to the ones we observed in Section "Effect of initial particle position". Zhang et al. 46 reported oscillations in the trajectory of glass beads in a cross-slot channel with a square cross-section at Re = 80. The trajectories of the rigid particle in our study show a good qualitative agreement with those trajectories seen under steadyvortex conditions published previously 31,46 .
Our study provides a deeper understanding of the dynamics of a rigid particle interacting with the vortex in a cross-slot junction at moderate inertia. The macroscopic observables (N, A, L, T res ) used in our work might help distinguish between particles with different properties in experimental cross-slot flows. To this end, a future investigation of the effect of channel geometry, in particular aspect ratio and Reynolds number, appears beneficial. Since most inertial microfluidic applications involve deformable cells, a study of deformable particles in cross-slot flows would provide further insight into the particle dynamics of real-world applications.

Lateral equilibrium positions in straight channels
Test case: effect of initial position and confinement on lateral migration in square duct. We determine the lateral equilibrium position of a rigid neutrally buoyant spherical particle in a Poiseuille flow through a straight duct with a square cross-section at Re = 100. The duct has a width and height W = 50Δx. The particle radius is varied such that confinement values χ = 0.1, 0.2 and 0.2857 are obtained. Periodic boundary conditions are used at the inlet and outlet of the duct. With L = 6W, the domain length is sufficiently long to ensure that periodic images of the particle do not interact. The viscosity is set to ν = 1/30Δx 2 /Δt.
We compare the equilibrium positions found using our IB-LB solver with those reported by Lashgari et al. 55 who originally proposed this case. Figure 9 shows excellent agreement of the equilibrium positions. As expected, equilibrium positions are independent of the initial position. Larger particles assume equilibrium positions closer to the duct centre at the mid-point of a duct edge.
Equilibrium positions in a straight channel equivalent to the inlet of the cross-slot. We perform simulations of a  Equilibrium positions versus confinement Fig. 10 Simulation of a rigid neutrally-buoyant spherical particle in a Poiseuille flows through a straight rectangular channel with the same aspect ratio as the inlet of the cross-slot geometry. The particle confinement χ is varied in the range [0.3, 0.75]. a Schematic of the set-up. The particle is initialised on the grey plane, which is located midway between the side walls. b Equilibrium positionsẑ eq ¼ z eq =ðH=2Þ as function of confinement χ. The particle remains on the mid-plane (y = 0) at all times rigid neutrally buoyant spherical particle flowing through a straight channel (along the x-axis) with the same crosssectional shape as the inlet of the cross-slot geometry to identify the lateral equilibrium positions for a given particle size. The simulation parameters are the same as in Section "Cross-slot geometry and simulation set-up"; in particular, the Reynolds number is Re = 120, and the channel aspect ratio is α = 0.6. Periodic boundary conditions are used at the inlet and outlet of the duct. With L = 4W, the domain length is sufficiently long to ensure that periodic images of the particle do not interact 56 . The schematic in Fig. 10a illustrates the initial configuration of the particle in the channel. Figure 10b shows the identified equilibrium positions z eq that we used as the initial positions for our main study in Section "Results and discussion". Our simulations confirm the known observation that, for the studied parameter values, the equilibrium position of the particle is located midway along the long cross-sectional axis (y = 0) 57 ; the particle remains on this mid-plane during its lateral migration.