Analogue tuning of particle focusing in elasto-inertial flow

We report a unique tuneable analogue trend in particle focusing in the laminar and weak viscoelastic regime of elasto-inertial flows. We observe experimentally that particles in circular cross-section microchannels can be tuned to any focusing bandwidths that lie between the “Segre-Silberberg annulus” and the centre of a circular microcapillary. We use direct numerical simulations to investigate this phenomenon and to understand how minute amounts of elasticity affect the focussing of particles at increasing flow rates. An Immersed Boundary Method is used to account for the presence of the particles and a FENE-P model is used to simulate the presence of polymers in a Non-Newtonian fluid. The numerical simulations study the dynamics and stability of finite size particles and are further used to analyse the particle behaviour at Reynolds numbers higher than what is allowed by the experimental setup. In particular, we are able to report the entire migration trajectories of the particles as they reach their final focussing positions and extend our predictions to other geometries such as the square cross section. We believe complex effects originate due to a combination of inertia and elasticity in the weakly viscoelastic regime, where neither inertia nor elasticity are able to mask each other’s effect completely, leading to a number of intermediate focusing positions. The present study provides a fundamental new understanding of particle focusing in weakly elastic and strongly inertial flows, whose findings can be exploited for potentially multiple microfluidics-based biological sorting applications.

weakly viscoelastic regime, where neither inertia nor elasticity are able to mask each other's effect completely, leading to a number of intermediate focusing positions. The present study provides a fundamental new understanding of particle focusing in weakly elastic and strongly inertial flows, whose findings can be exploited for potentially multiple microfluidics-based biological sorting applications.
Keywords Elasto-inertial Á Particle focussing Á Analog tuning Á Reynolds number Á Weissenberg number 1 Introduction We study particle focusing in strongly inertial and weakly elastic flows, whose importance is defined by two critical dimensionless numbers: the Reynolds and Weissenberg numbers. The Reynolds number, quantifying the importance of inertia over viscous effects, is defined as Re ¼ qUcH=l, where q is the fluid density, Uc is the centreline flow velocity, H is the characteristic length across the microfluidic channel cross-section, and l is the shear viscosity of the fluid. The Weissenberg number is defined as Wi ¼ kUc=H, where k is the relaxation time of the polymer additives. The ratio between these two parameters gives the elasticity number El ¼ Wi=Re. In Newtonian fluids, it has been well recognized that in the laminar flows typical of microfluidic channels (Re [ 1 and Re\2300) [1], inertia does play a role and particles tend to migrate to equilibrium positions at distances of order 0:6 of the channel hydraulic diameter, R, from the center, closer to channel centerlines, as initially predicted by Segré and Silberberg [2]. Numerous studies exploited this effect in microfluidic applications, frequently referred to as ''inertial focusing'' of particles [3][4][5][6][7]. The equilibrium position in inertial microfluidic flows result from the balance of two opposing forces: (1) the wall-induced lift force arising out of the interaction between the particle and the adjacent wall, which directs the particle away from the wall and (2) the shear-gradient force, induced by the velocity profile curvature, pushing the particle away from the flow channel centreline and towards the wall [1,8,9].
Particle migration and focusing have also been extensively investigated for viscoelastic fluids in microfluidic channels [10][11][12][13]. It has been observed that in purely viscoelastic flows, the particles tend to flow through a narrower focused beam close to the centreline of the channel [14][15][16][17]. The phenomenon of single line focusing of particles has been observed in a wide range of microfluidic channel geometries that include circular, [18][19][20] and square channels [15,16,21] typically at low Reynolds numbers (Re\1) elastic flows. The added effect of inertia coupled with the viscoelastic effect leads to particle migration to the channel centreline and has been utilized in numerous applications for selective separation of particles based on physical parameters such as size and shape. As an example, under constant viscoelasticity but comparatively higher inertia for a bigger sized particle, the larger particle moves first towards the channel centreline and can thus be separated from smaller sized particles in sample. Various applications have been developed on the phenomenon of viscoelastic focusing such as microflow cytometry, [22,23] cell and particle focusing, [24] as well as droplet [25] and bubble [5] manipulation. A recent review by Stoecklein and Di Carlo [26] covers the underlying physics as well as an exhaustive set of applications related to elasto-inertial microfluidics.
A clear observation from all the studies involving the combined effect of inertia and elasticity, clearly points towards a bistability scenario, i.e. the particles either focus at the Segre-Silberberg annulus when inertia dominates or at the channel centreline when elastic effects gain an upper hand [27]. In this study, we explored a weak elasto-inertial regime experimentally, by diluting PBS solution with PEO concentrations ranging from 1ppm till 50ppm (El = 0.001-0.5). We found that particle focusing in elasto-inertial microfluidics is tuneable at this regime, and the particle focusing bandwidth is not restricted to the bistability scenario, but rather can be engineered to any analogue value that lies between the centreline focusing and the Segre-Silberberg annulus positions. The phenomenon has been reported in the past in spiral microchannels, where the presence of Dean forces is exploited for particle manipulation [28]. However, the same particle focussing behaviour has not been observed in straight channels (e.g. rectangular and circular cross-section straight channels) in the weak elasto-inertial fluid regime (El * 0.1 or less). This regime of viscoelasticity is clearly not very well understood as previous reports on viscoelastic fluids have focussed on elasticity numbers 1 or higher [29,30], which means the elasticity component was not negligible. Hence, particle tuneability in straight channels in the weak elasticity regime can lead to new insights on the focussing mechanisms for elastoinertial microfluidics in a straight circular geometry.
We support our hypothesis with both experimental observations as well as numerical simulations. For the numerical simulations, an Immersed Boundary method is used to predict particle velocities and trajectories [31], whereas the elastic component of the stress, due to the presence of the polymers in the fluid, is modelled by the Finite Element Non-Linear Extensibility Peterlin model (FENE-P) [32]. One of the reasons for choosing the FENE-P model over the popular Oldroyd-B and Giesekus models [33][34][35], which also have been used for studying viscoelastic flows, is the robustness of the FENE-P model at finite inertia even for Weissenberg numbers greater than 1. The three-dimensional fully resolved direct numerical simulations provide access to viscous and elastic stresses in the flow, which enable us to explain the entire dynamics of particle migration, particle stability and equilibrium for low elasticity flows in a circular pipe, in the intermediate Reynolds number regime. We believe our current study addresses a crucial aspect in the understanding of weakly viscoelastic flows in elasto-inertial microfluidics, and more importantly demonstrates the possibility to tune the particle position in microchannels by to addition of minute amounts of elasticity.

Experimental setup
The experiments were carried out in a single-inlet and single outlet silicon capillary tube with a diameter of 56 lm. A solution containing 10 lm fluorescent particles is pumped through the silicon capillary tube by a syringe pump (Nemesys, Centoi Gmbh). The applied flow rates reach up to 240 ll=min. Fluorescent polystyrene particles, 10 lm diameter (Fluoro-Max, Thermo Scientific) at a volume fraction of 0.1% were suspended in Newtonian and viscoelastic fluids. The Newtonian fluid comprised of Phosphate Buffer Saline (PBS 1x) aqueous solution containing 0.1% Tween20. The viscoelastic fluid is prepared by mixing a range of different concentrations (from 0:0001 to 0:05 wt%) of Polyethylene Oxide (PEO) (with an average molecular weight equal to 1 Â 10 6 Da) in a Phosphate Buffered solution (PBS). A homogenization process for up to 2 h follows the preparation of the solution, after which the experiments were carried out immediately. All the experiments are done at room temperature, i.e., 25 degree centigrade.
In the present work, a systematic parametric study is performed by varying the Reynolds and Weissenberg numbers. To this end, we need to introduce several parameters, e.g., the zero shear viscosity (g 0 ), the infinite shear viscosity (g 1 ), the shear-rate ( _ c) and the relaxation time (k) of the resulting PEO solution, as well as the channel dimensions and volumetric flow rates. The Bird-Carreau model [36,37] is used to calculate the effective viscosity g eff of the PEO solution at different shear-rates _ c, which is given by In a previous study, Ebagninin et al. [38] reported that, in dilute regimes of PEO concentrations, the zero shear viscosity g 0 varies as g 0 $ c 2:17 for concentrations c less than 0:1 wt%. Also, based on other results from the same study, it can be deduced that the zero to infinite viscosity ratio, g 0 g 1 , is approximately equal to 2 for PEO concentrations up to 0:1 wt% and g 0 g 1 % 3:33 for PEO concentrations between 0:2 and 1 wt%. Finally, the flow index n in the Bird-Carreau model is found by comparing viscosities at two different shearrates. The relaxation time k was computed based on a recent study [39] on weakly viscoelastic flows, which gave an empirical law: 2.2 Imaging technique For imaging the particle flow in the circular capillary, we used an inverted microscope (Nikon Eclipse TI, USA) with a CMOS camera (Andoe Zyla,USA) and a LED lighting system (Lumenor Spectra X LED,USA).
To capture the images and control the microscope movements, the Micro-Manager Open Source Microscopy Software was used, while the image processing was performed with the MATLAB software (R2018a, Mathworks, USA). In particular, we measured from the images the particle focusing bandwidth F (see Fig. 1), which is the difference between the width of the streaks of particles captured in the images and their diameter, normalized by the entire cross section of the channel. For example, say a spherical particle of dimension 0.18 times the diameter of the circular cross section focusses at the equilibrium position in a channel inertial flow, i.e. at a focusing position that is 0.3 times the diameter. In this case the center of this particle will be at 0.2D, and if viewed from the top the particle streak will appear across 0.11D to 0.29D. As the equilibrium position of particle inertial focusing is symmetrical about the centre in a circular cross section, so there should be a particle streak with its center at 0.8D as well, that stretches from 0.71D till 0.89D. Thus, the value of F (center to center distance) in this case will be 0.8D-0.18D = 0.62D. F will be 0 when the particle is focused at the center and have a maximum value of 0:82 if the particle sticks by the wall.

Numerical setup
We study the motion of rigid neutrally buoyant spherical particles in a straight pipe with a circular cross section. The fluid is incompressible and Non-Newtonian, and its motion is governed by the Navier-Stokes equations where u i is the velocity vector field, p denotes the hydrodynamic pressure, s ij the polymer stress tensor, f i the immersed boundary force to account for the presence of the suspended particles, P is the ratio of the solvent viscosity to the total viscosity (function of the concentration of polymers), and Re is the Reynolds number as introduced above.
The additional stress tensor s ij describes the Non-Newtonian behavior of the carrier fluid, modeled with the Finite Element Non Linear Extensibility Peterlin model (FENE-P) as commonly done in numerical simulations of polymeric fluids [40]. . In the FENE-P model, s ij is written as a function of the configuration tensor C ij defined as where L is the dumbbell extensibility, o ij the Kroeneker delta and Wi the Weissenberg number. The configuration tensor is a symmetric second-order tensor, solution of the following dynamic equation for its 6 independent components The previous equation represents a balance between the advection of the configuration tensor on the left-hand side, and the stretching and relaxation of the polymer, represented by the first two terms and the last one on the right-hand side, respectively.
The motion of a neutrally buoyant rigid spherical particle is dictated by the Newton-Euler equations [31] given as: where m p and I p are the mass and moment of inertia of the particle, and U p c and X p c its center velocity and rotation vectors. The integrals are performed over the Fig. 1 Fluid flow setup for the experimental and numerical simulations in a circular microcapillary (drawn to scale). The ratio of the particle to pipe diameters is approximately 1 : 5 for both the experimental and numerical cases. The cross-sectional view shows the existence of possible multiple particle equilibrium states in very dilute pressure driven flows in a circular microchannel. The top view defines F, the focusing bandwidth surface of the particles oVp with unit normal vector n; r is the vector connecting the particle center to its surface. The no-slip and no-penetration boundary conditions on the suspended particle are enforced through the body force f i found by a direct forcing Immersed Boundary Method [41,42].
The boundary conditions on the curved pipe surface are enforced by a volume penalization method [43], while periodic boundary conditions are imposed in the streamwise direction. The numerical domain size is 15D Â 5D Â 5D, Dbeing the particle diameter, and is discretized with 480 Â 164 Â 164 grid points in the streamwise and crossflow directions, respectively. The equations of motion are solved with a fractional-step method with the low-storage third-order Runge-Kutta scheme for time integration and the second-order finite difference scheme for the spatial derivatives. Note that, the fifth-order Weighted Essentially Non-Oscillatory scheme [44] is used for the advection terms in the evolution equation for the conformation tensor. This numerical method has been used in the past for several applications with appropriate modifications, which include simulating spherical and oblate particle trajectory in a Newtonian fluid in microchannels [31], as well as elastic particles [45,46] and cell membranes in shear flows [47]. The interested reader is referred to the work by Izbassarov et al. [48] and Rosti and Brandt [49] for a detailed description of the numerical methods used to simulate multiphase flows in various Non-Newtonian fluids. Finally, we perform box-size and resolution studies to ensure that our results are not affected by any numerical artefact in a previous study [31].

Experimental analysis
The particle behaviour in circular capillaries is studied by examining the amount of focusing obtained in a Newtonian and Non-Newtonian flow, covering a wide range of flow-rates and ppm concentrations. Past studies have reported in detail that due to symmetric cross sectional stresses in a circular geometry, the particles settle on a symmetric annular radius, both in Newtonian [2] and viscoelastic flows [23]. As it can be seen in Fig. 2, we measured the area occupied by the particles in the channel and found that it corresponds to 78% of the entire channel. The maximum particle streak centre to centre distance is denoted by F and is called the focussing bandwidth, as mentioned earlier.
F is computed as the distance between the particle streaks, minus the particle diameter, normalised by the channel size. Since the ratio of the particle to channel dimension (10 lm and 56 lm, respectively) is 18%, this implies that the maximum particles centre to centre distance is F ¼ 60% (78 À 18%) for the Newtonian case. F ¼ 0:6D corresponds well to the theoretical predictions and to the previous experimental observations reported in the. literature for the equilibrium position of a particle in a Newtonian inertial flow [2].F ¼ 0:6D translates to an equilibrium position at 0:6R with R ¼ D 2 being the radius of the channel. Further, we observe that the focussing bandwidth of a particle shrinks (shrinking focusing bandwidth) towards the center due to the (5 ppm) reported in Fig. 2b, we observe that the particle focussing bandwidth is 0:5F instead of 0:6F for the same flow rate of 180 ll=min. Thus, the elastic effects of Non-Newtonian fluids are visible even with a slight increase in the elasticity of the fluid at high Reynolds number. This is because with an increasing Reynolds number, the Weissenberg number also becomes significantly prominent due to an increasing shear rate. A similar trend is observed for even larger elasticity concentrations, where at 50 ppm a more compressed particle focusing bandwidth (Fig. 2c) is found at 0:4D for the same flow rate as the Newtonian case (180 ll=min; 50 ppm). Finally, Fig. 2d shows a centreline focussing scenario where the Reynolds number and Weissenberg number are comparable, but both are very low (\1), which results in the dominance of elastic effects. Similar results on centreline particle focussing have also been reported in previous studies related to elasto-inertial fluids [10,50].
We tested the particle behaviour for various concentrations of PEO in PBS and for different flow rates in order to study the dual effects of the Reynolds and Weissenberg numbers on the particle focusing. It has previously been reported that in Newtonian fluids for circular cross sections, the particles achieve a focused state at a Reynolds number of 25 onwards [51]. We varied the concentration of PEO from 1 to 50 ppm and the flow rate from 5 to 240 ll=min (Re = 1-100). In Fig. 3 we report the focusing bandwidth F as a function of the Weissenberg number, at four different Reynolds numbers (Re ¼ 25, 50, 75 and 93).
We observe that the focusing bandwidth reduces as the Weissenberg number increases, i.e. particles tend to approach the microcapillary center, for all the Reynolds numbers. On the other hand, at a fixed Weissenberg number, the particles move away from the center with increasing Reynolds number for all the cases considered. While imaging the cross sectional distribution of the particle is out of the scope of this work, a previous report has demonstrated the distribution of particles to be symmetric about the channel cross-section for all Reynolds number [51] in viscoelastic fluids.

Numerical simulations
In order to gain more physical insight on our experimental observations, we performed three-dimensional direct numerical simulations of a rigid spherical particle in a polymeric flow modelled by the FENE-P model. This model allows to reach higher Wi in dilute polymer flows than other viscoelastic models such as the Oldroyd-B [52] and better prediction of shear viscosity at high shear rates than Giesekus fluids [53] and Phan-Tien Thanner models [54].
The simulations consider a single particle, in both the Newtonian and Non-Newtonian case, and are started from a fully developed single-phase flow obtained at the selected Re and Wi. In our simulations, a dilute particle suspension is therefore studied, with a single rigid spherical particle with diameter to channel size ratio equal to 1:5, similar to the experimental conditions. Also, in the viscoelastic cases, the dumbbell extensibility L is fixed equal to 60 and P to 0:9.
The Reynolds number is varied between 35 and 600, and the Weissenberg number between 0 and 3. In the Newtonian case, we found that the particle equilibrium position stabilizes at 0:6R for a Reynolds number of about 100, as shown in Fig. 4. This is consistent with the experimental measurements and the Segre-Silberberg theory which predicts that in an inertial laminar flow and in the absence of elasticity, the particles focus at a distance 0:6R from the center of the channel. From our results we can observe slight variations of the particle focusing position as Re increases, with the equilibrium position approaching more towards the wall at moderate Fig. 3 Effect of the Reynolds and Weissenberg numbers on the particle focusing. Particle focusing as a function of Wi at four differentRe. In all cases, the particle moves closer to the center as Wi increases Re Re ¼ 10 À 150 ð Þ and then reducing again for Re [ 150. Note that, the particle focuses at approximately 0:53R at Re ¼ 35, which is in very good agreement with the theoretical prediction by Asmolov et al. [55], who reported the same equilibrium position for a finite size neutrally buoyant particle of the same size at the same Reynolds number, and also agrees with the experimental results shown in Fig. 3. In a Non-Newtonian solution, an inverse Segre-Silberberg scenario has been proposed in the past to explain particle migration towards the walls, as well as channel centreline, at low Re and Wi [14]. This is due to the distribution of first normal stress difference which is minimum at the channel centre in case of circular cross section channels and minimum at the centre and walls in case of square and rectangular cross sections [50]. In an elasto-inertial regime with appreciable elasticity and inertia, the corner positions disappear and only centerline focusing is observed. However, we observe that under appreciable inertia and weak elasticity, the elastic effects no longer exert complete dominance on inertia.. Instead, multiple equilibrium positions exist that depends on competitive domination of inertia and elasticity of the carrier fluid.
The simulation results (Fig. 4) are in good agreement with our experimental measurements (Fig. 3). The particles move away from the center until a critical Reynolds number of about 200 for the Newtonian cases as well as for Wi ¼ 1, with a resulting particle focusing position close to 0:6R.
With increasing Weissenberg numbers (Wi = 0, 1, 2, 3), the particles move closer to the center at all Reynolds number which matches with the experimental observations in Fig. 3. While at Wi ¼ 0 and 1 inertial effects dominate over elasticity and the particle settles close to 0:6R, above Wi ¼ 2 elasticity effects start dominating displacing the particle consistently closer to the center as Re Increases (see the Supplementary movies S1 and S2). Note that, for each Wi there is a critical Reynolds number Re M where the particle has a maximum distance from the center, i.e., the. particle is closer to the walls, and that Re M is smaller for higher Weissenberg numbers. This is not so visible till Re ¼ 100 which is our experimental observation range, but the trend is clearly seen in the simulation results. Beyond this critical Reynolds number Re M , the particles start moving closer to the center with increasing effect of inertia, and the probable cause of this discrepancy beyond the critical Reynolds number ðRe M Þ can be due to dominance of particle size effect as well as the curvature of the pipe [56]. Figure 5 shows the first normal stress difference as a function of the Weissenberg number, and also the interaction of the particle and the polymeric component of the fluid to generate a distorted and nonsymmetric field. A stretching similar to that of a rubber band is observed due to the first normal stress difference around the particle, which causes the particle to move closer to the center in an oscillatory fashion with increasing elastic effects. The first normal stress difference increases with the Weissenberg number, i.e., with the fluid elasticity, thus inducing larger particle displacements. Next, we study the entire particles trajectory, from their initial position (0:3R) to their final equilibrium positions, for a fixed Reynolds number Re300 and different Wi (Fig. 6a).  For all the Weissenberg numbers, the particles right upon introduction have a short transient motion towards the center of the channel. The short transient motion can be attributed to the stress build-up around the particle right upon introduction in a fully developed baseflow, as well as to inertial effects. Interestingly, while the particle moves smoothly in the Newtonian fluid, in the polymeric flow, the final equilibrium position is reached through oscillations which slightly increase in amplitude with Wi and Re, and are in the order of R=10. Note that, the particle in the Non-Newtonian fluid oscillates also when its final equilibrium position is reached (Fig. 6a, Wi = 1,2,3). This is due to the distorted first normal stress difference that builds up around the particle leading to an oscillatory stretching around the particle as shown in Fig. 5c, see also [57]. The final equilibrium position of the particle is also shown to be independent of its initial condition for both the Newtonian and the Non-Newtonian fluid, as shown in Fig. 6b, where the results obtained from two different initial conditions for the particle are compared (0:3R and 0:7R).
Finally, we also performed preliminary numerical simulations at fixed finite Reynolds and Weissenberg numbers, Re300 and Wi ¼ 2, in a channel with square cross section, in order to compare the results with those in a circular pipe. In the square channel, different equilibrium positions are reached when the particle is initialized at two different azimuthal positions (vertical axis and diagonal), even if they are at the same distance from the center. This is due to the loss of periodicity in the azimuthal direction induced by the presence of the four corners.
In Newtonian fluids, particles tend to settle in one of the four positions that are at right angle from the centre, i.e., along the horizontal and vertical lines passing through the centre, even at non-zero and finite Reynolds numbers.
On the other hand, it has been previously shown that in viscoelastic fluids with no inertia, particles tend to focus in the centre of the channel and also in one of the four corners [50]. The differences between the two geometries can be explained by considering the first normal stress difference, reported in Fig. 7a and Fig. 7b. While in the circular cross section (Fig. 7a) case, the first normal stress difference is uniform in the azimuthal direction, this is not the case in the square cross section (Fig. 7b) case, where it shows a strongly non uniform distribution. In particular, it is maximum at the wall in the intersection with the horizontal and vertical lines passing through the center, while it is minimum in the channel center and in the four corners. Figure 7c shows a comparison of the final equilibrium positions in the cases of square and circular cross sections in presence of both finite inertia and elasticity. In both the geometries and for all the initial positions considered, the particle first moves towards the center of the channel and then towards the walls, eventually reaching the equilibrium positions.
The initial motion towards the center can be explained due to stress build-up around the particle upon introduction as well as the effect of inertia, as in the circular capillary.
Note, in the square duct (Fig. 7c) two different equilibrium positions are found: one on the horizontal and vertical lines at right angle from the center, similarly to Newtonian fluids, and the other one on the main diagonals, similar to viscoelastic fluids. The latter case is really interesting as in a similar case for a Newtonian fluid it was earlier observed that the particle first reaches an equilibrium manifold and then undergoes slow lateral migration to the vertical or horizontal axis [9,31]. However, due to the finite inertia considered in our simulations, the particle does not reach the corner, but settles at a certain distance from it, and the symmetrical first normal stress difference about the diagonal axis prevents the slow migration of the particle along the equilibrium manifold [31]. This result is similar to what observed in the case of deformable capsules in a Newtonian fluid [52]. The interested reader is referred to the work by Trofa et al. [58] and Lashgari et al. [31] for further analysis on the different migration dynamics of particles in square channels.

Conclusions
We have reported results of a combined experimental and numerical effort to characterize the dynamics of very dilute suspensions of particles in a circular microcapillary for both Newtonian and Non-Newtonian fluids. In particular, we study the phenomenon of particle elasto-inertial focusing in presence of moderate inertia and weak fluid elasticity and found a general tendency of the particle to assume multiple equilibrium positions based on competitive effects of inertia and elasticity. The results demonstrate for the first time, the possibility of analog tuning of the particle focusing positions in elasto-inertial flows.
By introducing in the carrier fluid a small amount of viscoelasticity of the order of a few ppm of PEO solution, it is possible to tune the position of the particle at any desired position in the microcapillary that lies between F ¼ 0:6R and F ¼ 0.
The experimental evidence is supported by wellresolved three-dimensional numerical simulations, where the FENE-P model is used to account for the viscoelastic behavior of the flow. From the numerical simulations, we show the entire particle migration process and extend the predictions to moderate viscoelastic regimes where the inertial effect is predominant (Re [ 100Þ. We believe the possibility of particle tuning in weakly elasto-inertial regimes will greatly aid in a better fundamental understanding of elasto-inertial focusing in microfluidic channels in a new perspective.
Author contributions IB, MER, LB and AR are responsible for the conceptualization of the idea. IB, MER and TK are responsible for the investigation, specifically performing the experiments and simulations and to develop the methodology for the data analysis. AR and LB are responsible for funding and resources acquisition. AR is responsible for the project administration and for the work supervision. IB did the data curation and formal analysis, and is responsible for the preparation, creation and presentation of the published work, as well as for the data presentation. IB and MER wrote the original draft and all authors are responsible for its review and editing. All authors have read and approved the manuscript.