Olga Stamati

Benjy Marks

Edward Andò

Stéphane Roux

Nathanaël Machicoane
email: nathanael.machicoane@univ-grenoble-alpes.fr

X-ray radiography 4D particle tracking of heavy spheres suspended in a turbulent jet

Keywords: particle suspension, particle-laden turbulence, X-ray radiography, 3D tracking, 4D kinematics

The suspension of a heavy sphere by an upward jet is a classical fluid mechanics experiment to demonstrate the fluid forces acting on an object.

In the range of the parameter space where the sphere can be suspended, the dynamics can either be regular, i.e., with oscillations around an equilibrium position, or chaotic, with extreme events leading to large deviations from that equilibrium region. The existence and characteristics of suspension regimes of several heavy spheres in such flow configurations remain open questions. Spheres compete for the equilibrium position and come very close to each other, resulting in large local particle concentrations that prevent direct imaging. Relatively high speed X-ray radiography along with the radioSphere analysis technique is leveraged here to study the time-resolved 3D trajectory of each individual sphere in a vertical jet. radioSphere is

Introduction

The suspension of a sphere by a jet of a fluid lighter than the sphere is a problem whose study dates back to O. Reynolds [START_REF] Reynolds | On the suspension of a ball by a jet of water[END_REF]. Despite its widespread use in fluid mechanics education and outreach [2; 3; 4], this situation received fairly little attention. The sphere's suspension is explained in terms of the Coandȃ effect [START_REF] Wille | Report on the first european mechanics colloquium, on the coanda effect[END_REF], through an attractive force exerted on the sphere resulting from the deflection of the jet by the sphere [START_REF] Vil | Gel'Mi, Jet flow around a sphere[END_REF]. Since the early studies, experimental and theoretical approaches have been led for upward vertical and tilted jets with heavy spheres and a vertical downward jet with light spheres [7; 8; 9; 10], up to recent studies [START_REF] Mejia-Alvarez | A levitating/rotating ball-basic considerations[END_REF]. This work investigates heavy spheres suspended in a vertical upward jet, but with several spheres interacting in the fluid, a situation which has not been explored in the literature to the best of the authors' knowledge.

Determining the existence and location of the equilibrium position of a single sphere in a jet poses an apparently simple and canonical problem, which is not easily resolved. Already with a vertical jet, i.e., when the buoyancy force and the drag force are in the same plane with opposite directions, the flow inhomogeneity around the sphere prevents the use of readily available drag coefficients (e.g., Schiller-Naumann model [START_REF] Clift | Bubbles, drops, and particles[END_REF]).

Using the empirical knowledge of the incoming flow below the sphere, [START_REF] Davoust | The suspension of a sphere in a turbulent jet[END_REF] derived a formulation to predict the equilibrium position of the sphere in the vertical jet, but failed at capturing the dynamics when extreme events are present. While it is generally accepted that the equilibrium position cannot exist very far from the jet exit, they classified the existing equilibrium positions in different regimes. In the first two regimes, the sphere oscillates around an equilibrium position, with a larger frequency in the horizontal plane than in the vertical plane, as predicted by [START_REF] Mcdonald | Levitating beachballs[END_REF]. Both regimes are distinguished by the size of the sphere compared to the jet's width at the sphere equilibrium location. This ratio decreases as the equilibrium position moves away from the jet exit, e.g., by increasing the mean jet exit velocity or by changing the sphere-to-fluid density ratio toward lower values. In the first regime, termed large sphere behavior, the sphere is larger than the local jet width, while it is smaller in the second regime, referred to as the regular regime. As the sphere equilibrium position moves further away and the attractive force further decreases, extreme events are observed, resulting in large deviations with respect to the equilibrium point; this regime is termed "chaotic". The decay of the attractive force along the jet axis can be seen as a decrease in the equilibrium region's characteristic size. At large jet velocities, the sphere is hence suspended far from the nozzle and turbulent fluctuations can easily push the sphere away from the stable region and create chaotic dynamics. The first motivation of this work is to study the competition of several spheres for this equilibrium region, as the number of spheres increases and as the size of this region decreases, spanning the regular and chaotic regimes of a single sphere.

The current study poses the experimental issue of measuring the timeresolved 3D dynamics of multiple objects that evolve in close vicinity. While tracking many thousands of particles in 3D is possible with modern highspeed cameras and algorithms (e.g., [START_REF] Schanz | Shake-the-box: Lagrangian particle tracking at high particle image densities[END_REF]), as long as the inter-particle distance remains relatively small with respect to their apparent diameter and their frame-to-frame displacement, this is not the case here. When several spheres are suspended in a turbulent jet, contacts between spheres occur often and typically involve more than just two spheres. The approach chosen here is hence to leverage a recently proposed X-ray radiography algorithm, radioSphere [START_REF] Andò | Single-projection reconstruction technique for positioning monodisperse spheres in 3d with a divergent x-ray beam[END_REF], that retrieves the full 3D information from a single 2D projection under certain assumptions. In X-ray radiography, an X-ray beam is partially attenuated as it passes through a sample and this attenuation is recorded onto a 2D plane. While X-ray radiography is often resorted to for the study of multiphase flows [16; 17], it is mostly used and devel-oped in solid mechanics and material science, with applications ranging from manufacturing to medicine. In particular, tomography reconstructs the 3D map of the internal structure of complex systems through many 2D projections obtained by X-ray radiography at various angles around the sample (e.g., [18; 19]). However, the number of projections needed in this method greatly limits the acquisition rate; reconstruction the slice radiographies into a single 3D volume is under the assumption of no movement. To overcome this limitation, radioSphere uses prior knowledge of the sample (spheres) and imaging geometry to retrieve a 3D particle position from a single 2D radiograph. In particular, it is suited to the study of mono-dispersed spheres of known diameter in a typical laboratory X-ray set-up, in any complex arrangement, and uses the geometrical magnification resulting in the combination of a divergent X-ray beam and the relative distance between each sphere and the detector. The second motivation of this work is to validate and use this approach in a dynamic particle-laden flow setting, as radioSphere was so far only validated using synthetic data and a static granular assembly [START_REF] Andò | Single-projection reconstruction technique for positioning monodisperse spheres in 3d with a divergent x-ray beam[END_REF]. Section 2 presents the experimental set-up and measurements methods.

The application of radioSphere for the analysis of the radiographs and the characterization of the suspension regimes are detailed in Section 3. Section 4 reports the dynamics of the suspended spheres and changes in regime from freely advected spheres to collision-driven motions. The results are discussed in Section 5, drawing a general qualitative picture and introducing a framework for the establishment of suspension regimes of a collection of spheres suspended in a jet. This is followed by a conclusion (Section 6).

Experimental methods

Materials

The experiment aims to study the 4D kinematics of the suspension of spheres in a turbulent jet inside a cylindrical tank. The choice of the fluid and the sphere properties is a compromise between their relative density and the jet velocity. From an imaging point of view, a material that attenuates more X-rays compared to its surrounding fluid is needed, so as to be clearly visible in the radiographs. From a fluid mechanics point of view, the drag force applied from the fluid to the spheres should be large enough to exceed Archimedes' force and in turn suspend the spheres. The best compromise between the above lead to the selection of 10 mm soda lime glass spheres suspended in water, although compared to Refractive Index Matched Scanning (RIMS) this is more flexible considering that contrast agents can be used in particles and there is no need to match refractive indices.

Experimental set-up

The in-situ acquisition is performed inside the X-ray scanner of Laboratoire 3SR (Grenoble, France), as shown in the schematic set-up of Fig. 1. A cylindrical PMMA (i.e., low X-ray absorption) tank of 50 cm in height and 6 cm inner diameter is bolted on top of the rotation table and sealed with an O-ring. The top of the tank is left open, while a hole is drilled at 4/5 of its height to allow the circulation of water in the system. A 3D-printed piece is placed at the bottom, which has a 4.5 mm cylindrical hole in the middle. This base controls the input velocity, while it retains the beads inside the container. It is also slightly tapered, with an angle of 6.3 • to prevent stag- Note that the system is not rotated and radiographs are acquired only in the orientation illustrated. Spheres that are at altitudes z < d s are labeled as grounded, while spheres at higher altitudes are deemed suspended.

nation of spheres at the bottom of the tank against the cylinder's wall. This ensures a constant number of spheres during the entire duration of a given experiment.

A magnetically driven centrifugal pump is used, connected in parallel, to produce the water flow. The input velocity is measured through a flowmeter, placed before the tank entrance. Two flow rates are examined: 1.7 L/min and 2.2 L/min, resulting in jet mean exit velocities U j = 1.8 and 2.3 m/s and jet Reynolds numbers Re j = U j d j /ν l = 8950 and 11600, with ν l the kinematic viscosity of water, exploring two suspension regimes in the case of a single sphere (regular and chaotic, see Section 3.2). Note that the maximum number of spheres that can be suspended at 1.7 L/min is 12, and this will be the configuration with the maximum number of spheres tested for both injection rates. Other experimental parameters remain fixed, so the results will be presented in terms of jet Reynolds numbers, but Table 1 provides alternative dimensionless numbers for comparison. In addition, the sphereto-fluid density ratio is K = ρ s /ρ l = 2.51 and the sphere-to-jet diameter ratio is D = d s /d j = 2.22. The sphere Reynolds number Re s = U s d s /ν l = 7010 is based on the sphere diameter d s and terminal velocity based on free-fall in quiescent water U s = 4/3d s /C d g(ρ sρ l )/ρ l ≃ 0.63 m/s, where the drag coefficient C d is approximated to be 0.5 and g is the gravitational acceleration.

Scanning geometry and acquisition parameters

The current and voltage of the X-ray source are set to 500 µA and 150 kV, respectively, while the source operates at a large spot mode. The beam is strongly hardened with a 2.4 mm Cu filter. Radiographic acquisition is performed at the detector's highest speed setting, which is 60 Hz, imposing a 4 × 4 binning and an effective pixel size of 0.598 mm on the detector panel.

The time resolution limitation unavoidably causes motion blur artifacts in the acquired radiographic projections, especially for the highest injection rate.

The treatment of the motion blur is specifically taken into account during the 4D sphere tracking, as discussed in Section 3.1.

For both input velocities, the source-detector distance is 766.3 mm, with the source-object distance (i.e., center of cylindrical tank) set to 272.7 mm for the lowest and to 400. field of view since spheres are suspended at higher distances from the tank's bottom. This means that a smaller magnification level is achieved for this configuration, resulting in smaller variation in the projected sphere sizes.

In particular, for the lowest velocity the minimum and maximum projected diameters are 49 and 60 px, while for the highest one they are 31 and 35 px, respectively. These small size variations are leveraged by the radioSphere technique to position the spheres in 3D space, as discussed in Section 3.1.

Test procedure

In the beginning of each test, the tank is filled with water and a "flatfield" image of the entire system without spheres is acquired by averaging 100 projections. Spheres are then inserted into the tank, and with the pumps still off, another set of 100 images is recorded. This averaged projection with optimal imaging conditions (very long effective exposure for noise reduction, no motion artifacts) enables a good first 3D guess of the static spheres' position, as discussed in the following Section. Pumps are then turned on and 2000 radiographs are continuously recorded at 60 Hz.

The acquisition encompassed both the transient from static spheres at rest in the quiescent fluid and the steady state, where the spheres are suspended in the turbulent jet. The acquisition duration in the steady state is approximately 17 s. In this experimental configuration, this means that the very first part of the acquired data corresponds to a transient, where the jet velocity starts at zero until reaching a constant value. This transient is removed from the analysis, which focuses on the steady-state dynamics.

Analysis

Application of the radioSphere technique

A recently developed X-ray radiography-based technique is employed to reconstruct the 3D position of each sphere: radioSphere [START_REF] Andò | Single-projection reconstruction technique for positioning monodisperse spheres in 3d with a divergent x-ray beam[END_REF]. It is a twostep approach that enables the 3D positioning of assemblies of spheres from only a single radiograph. radioSphere is based on the strong a-priori knowledge of the particles' shape (spherical) and size (mono-dispersed), combined with the imaging geometry. In the tested granular assemblies, naturally, some spheres are positioned closer to the X-ray source with respect to the others. Thanks to the geometrical magnification resulting from the divergent X-ray beam, this means that despite having the same physical size (monodispersed) their projected sizes on the detector are different. radioSphere leverages these small size variations to position spheres along each beam ray.

A flowchart of the technique is shown in Fig. 2 for the assembly of 12 stationary spheres at the lower input velocity. First, the recorded projection is converted into a calibrated path length image (in mm) through a calibration procedure described in Sec. 6. The first step of radioSphere is based on a set of fast Fourier transform (FFT) deconvolutions between the measured radiograph and a series of structuring elements. For a given magnification level, a structuring element is a full-scale projection of a single disk (same physical size as the tested spheres) centered on the detector panel (see Fig. 2d). The structuring element can be considered as a template to be matched, or a shape function (which here only varies in its size) for the Fourier deconvolution. The position of a sphere along the X-ray beam can be detected by performing a number of FFT-based deconvolutions for different magnification levels, which requires a series of structuring elements of varying projected sizes accounting for positions from closer to the source to closer to the detector1 . Here, the limits of the examined zoom positions are directly given by the scanning geometry and the physical limits of the cylindrical tank. At the end of this first step, a 3D initial guess of the particle centers is obtained which, based on the scanning geometry, is projected to yield an estimated synthetic radiograph (see Fig. 2 top right).

The second step consists in refining this initial guess through an iterative minimization algorithm. The quantity to be minimized is the pixel-wise squared residual, defined as the difference between the measured radiograph and the gradually updated synthetic one. More precisely, the iterative algorithm aims to minimize the current residual as a linear combination of three synthetic residual fields, which are the perturbations of the current guess of each sphere in each orthogonal direction. The convergence criterion is set as the norm of the difference between the 3D position guess of two successive iterations, which is set here as 1/10 of the sphere diameter. Given a good initial guess, the algorithm converges after a few iterations.

As already mentioned, the minimization step requires an initial guess of the particles' positions. A source of initialization can be a 3D X-ray tomographic scan, which here would have been very impractical given the flow loop arrangement. Another source of initialization can be the deconvolution step, or in the case of a time series analysis, the previous time frame. For the presented experiments, an initial guess coming from the previous frame is sufficient to achieve a good convergence for the current one. This means that for each tested configuration, the deconvolution step is only run once, for the initial frame containing the stationary spheres. This choice is justified by the fact that there is no motion blur present, but also a higher signal-to-noise ratio (SNR), since 100 images are averaged to produce the standing frame.

For all the remaining frames the optimization step is directly run, using as an initial guess the position of the previously converged frame.

The identification of particle positions in subsequent time steps directly leads to 3D time-resolved particle traces. The measured particle displacements can be used as additional constraints in the minimization process to overcome the potential motion blur artifacts caused by the detector speed limit. For this, a second optimization round is run, where before calculating the pixel-wise residual, the modelled projection of each sphere is convoluted with a unique step (or hat) kernel, the size of which corresponds to the magnitude of the sphere's displacement, while its direction corresponds to the one of the displacement vector. This second run of the optimization accounting for the spheres' displacements leads to clearly reduced residuals reaching values of 1/10 of the sphere diameter, as shown in Fig. 3.

The described application and adaption of the radioSphere technique is its first implementation in a dynamical system. It should be noted here, that mainly for the highest input velocity, the tracking procedure required a certain degree of manual verification in its implementation. More specifically, for frames in which some spheres exhibited large displacement increments, a trial and error approach was necessary to select the right amount of the a)

x-ray attenuation applied perturbation around the sphere's current position to ensure the convergence of the minimization, and, in turn, a successful tracking. A potential development of the technique for similar dynamic cases will be to include the velocity vector directly inside the minimization functional of the pixel-wise residual, improving the generality and robustness of the method.

Applying radioSphere on the radiograph time-series results in a 3d trajectory for every individual sphere. The output data, N s time-series of x, y, and z for each operating condition, is available at https://doi.org/10.5281/zenodo.7438422 [START_REF] Stamati | X-ray radiography 4D particle tracking of heavy spheres suspended in a turbulent jet[END_REF].

Suspension regimes

A parametric study considered the suspension of spheres in a vertical upward jet with variations of the jet and spheres diameters, mean jet exit velocity, and exploring two different spheres densities [START_REF] Davoust | The suspension of a sphere in a turbulent jet[END_REF]. They reduced the number of parameters into a single non-dimensional group, a jet-based Froude number F j = U j / gd j (different from the Froude number F r defined above, see Table 1), varied in the range 100 < F j < 500, to establish the following suspension regimes: 120 < F j < 175: large sphere behavior, 175 < F j < 250: regular regime,

F j > 250: chaotic regime.
The trajectory of a single sphere in the considered configuration reveals a regular regime for Re j = 8950 and a chaotic one for Re j = 11600 (Fig. 4a)).

While the mean altitudes, respectively ⟨z/d j ⟩ = 9 and 15.5, corresponding to ⟨z/d s ⟩ = 4 and 7, are compatible with the ones reported in [START_REF] Davoust | The suspension of a sphere in a turbulent jet[END_REF] for these regimes (see Fig. 4 in [START_REF] Davoust | The suspension of a sphere in a turbulent jet[END_REF]), the regimes occur in a very different region of the parameter space. Here, the sphere-to-jet density ratio is K = ρ s /ρ l = 2.51 and diameter ratio is D = d s /d j = 2.22, while they considered K = 85 or 700 and 1.5 < D < 7, resulting in high jet Reynolds and Froude numbers.

In the current study, F j = 5.7 results in a regular regime, and a value of 7.4 leads to a chaotic regime, well below the transition values indicated above.

Note that a regular regime is found in a similar region of the parameter space for a downward jet and light spheres [START_REF] Barois | Equilibrium position of a rigid sphere in a turbulent jet: A problem of elastic reconfiguration[END_REF]. This suggests that for small values of K, the jet Froude number may not be the most pertinent nondimensional group to distinguish transitions among suspension regimes (or that the transition values have a non-trivial dependency on K). In what follows, the jet Reynolds number will be used for simplicity to distinguish between regular and chaotic regimes, as only two values are explored here, For the radial position r, a broad peak around f r = 6.75 Hz is highlighted for the lower Reynolds number, while a small peak for the vertical position z is shown around f z = 2.1

Hz.

Re j = 8950 and 11600 respectively, and the focus is on the effect of varying the number of spheres in both cases.

In the regular regime at a jet Reynolds number Re j = 8950, the single sphere remains in the close vicinity of the equilibrium position, and oscillations are observed. Figure 4 reports the power spectral densities (PSD) of the radial r and vertical z time series. The radial position is obtained from the two horizontal positions x and y resulting from the 3D tracking, through a change from a cartesian (z, y, x) to a cylindrical (z, r, θ) coordinate system. The PSD of the radial position shows a broad peak centered around f r = 6.75 Hz. The vertical position spectrum appears as monotonously decreasing, but a small peak can be visible at f z = 2.1 Hz which corresponds to the vertical oscillations. While a vertical oscillation frequency lower than the horizontal one is expected, the ratio f r /f z is approximately twice that in the prediction of [START_REF] Mcdonald | Levitating beachballs[END_REF] (3.2 versus 1.4), but smaller than the ratio of about 5

found in [START_REF] Davoust | The suspension of a sphere in a turbulent jet[END_REF].

For a higher jet Reynolds number value, Re j = 11600, the sphere explores a large region around the equilibrium positions, and at times exits the stable jet region completely until being entrapped again (see for instance the leftmost part of the trajectory in Fig 4a ) where the sphere exited the jet's core at a high altitude, fell almost vertically, rebounded, and joined the jet's axis again). These extreme events are the signature of the chaotic regime and strongly affect the position power spectral densities. While the vertical PSD may have retained a small peak around f z = 2.2 Hz, the broad largeamplitude peak on the radial PSD is completely lost. The spectrum shows a plateau at low frequency, characteristic of fully decorrelated dynamics at long times, followed by a steep decay at high frequencies. Note that both radial and vertical spectra contain more energy at high frequencies in the chaotic regime than in the regular one.

Collective suspension dynamics

Observations and global metrics

The use of radioSphere on the time-resolved X-ray radiography measurements yields the 3D trajectory for each sphere present in a given measurement run. Fig. 5a) presents such trajectories, arbitrarily color-coded to distinguish different spheres, in the case of 3 spheres at a jet Reynolds number Re j = 11600. The presence of more than one sphere results in a much broader exploration of the flow volume. While not directly visible, already with only 3 spheres, the spheres are often found in close vicinity and even in direct contact, which will be further discussed in Sec 4.2.

As visualizing many concomitant 3D trajectories over long durations is not straightforward (long meaning tens of oscillation periods 1/f z ≃ 0.5 s here), two complementary approaches are proposed. On the one hand, Figure 5b-c) displays a single 2D trajectory from each measurement run with 1, 3, 6, and 12 spheres being suspended. This number is later referred to as the number of spheres, noted N s . In the example of N s = 6, one trajectory out of 6 is selected at random to be plotted in Fig. 5b-c), for both jet Reynolds numbers Re j = 8950 (b) and 11600 (c). On the other hand, Figure 6 only plots the time series of the vertical position but for every sphere present in a given measurement run. Note the change of scales between sub-figures of and 11600 (c) (F r = 6.3 and 8.16). The vertical axis is adapted as the jet velocity is increased, but the horizontal axes are fixed (extended up to the cylinder walls minus the sphere radius, so a marker on the axis corresponds to a sphere contacting the wall).

apparent change is that extreme events are present as soon as N s ⩾ 3, indicating a transition to a chaotic regime, induced by the competition for the equilibrium position. For N s = 3, while vertical positions below the equilibrium position (⟨z/d j ⟩ = 9 for a single sphere) and at larger radial positions (| y |> d j ) are observed, high altitudes seem easily reached. However, the vertical exploration happens in a smaller range for N s = 6, with the onset of a radial exploration at moderate altitudes. These changes are even more apparent when N s = 12, where only half of the previous altitude range is explored for that trajectory, and a large part of it is spent on the cylinder's bottom wall. This change of behavior, with reduced vertical exploration and increased radial exploration, is also observable for Re j = 11600. However, the larger jet velocity associated with an already chaotic dynamic for a single sphere leads to extreme events where high altitudes (comparable or greater than ⟨z/d j ⟩ = 15.5 for a single sphere) can be reached at times even for conditions with a large number of spheres.

By focusing on the vertical coordinate, Figure 6 shows the simultaneous dynamics of every sphere present in a given run, and qualitative comparisons can be made as the number of spheres (rows) and the jet Reynolds number (column) increase. The maximum altitude reached by a single sphere (dotted line) is only surpassed by an extreme event at N s = 3 and Re j = 8950 resulting from the collision between two spheres. At this condition, a sphere can still oscillate for a certain duration around the equilibrium position, but even here it is most often seen to visit its vicinity only shortly before falling back toward lower altitudes region (if not to the cylinder's bottom wall).

For the other conditions, with the exception of several extreme events, the This behavior is not observed for the vertical position, where the decreasing trend with an increasing number of spheres is followed by every sphere. In addition, this trend is independent of the jet velocity, as seen by the renormalization in Fig. 7c). The mean altitude is known to be driven by the drag force exerted by the jet on the sphere, which scales along F D ∼ U 2 j . The mean altitude rescaled by (U s /U j ) 2 (where U s is merely introduced to keep the quantity non-dimensional, as it is a constant here) collapses for the two studied jet Reynolds numbers, showing, in addition, independence to the suspension regime. In what follows, a cylindrical coordinate system, (z, r, θ), is adopted and the motions in the horizontal plane are discussed in terms of the radial coordinate r. Figure 7d) reports the mean radial position, which, as described, increases with N s . The change in the coordinate system brings a similar trend for the maximum and minimum of the time-average values, that surround well the ensemble average trend, despite showing a spread that increases with the number of spheres, in agreement with previous comments. The mean value of the radial and vertical coordinates can serve as a proxy for the location of the equilibrium position for the collection of spheres.

Similarly, the associated standard deviations can inform on the spatial extent of these regions. For a single sphere at Re j = 8950, this region is stable so the probable position is in the vicinity of the cylinder's bottom wall, except in the case of a single sphere which shows a bell curve centered on the equilibrium position. The PDF for N s ≥ 3 presents roughly decreasing exponential tails, with more and more pronounced slopes with increasing N s . In opposition, at the higher value of the jet Reynolds number, the PDF for a single sphere is characterized by a broad peak around the equilibrium position and a small one at a much lower altitude. This peak increases in value for N s ≥ 3 and moves closer to the bottom wall, and a plateau is observed for larger altitudes.

This plateau, showing equiprobability, extends up to 15d j for N s = 3, where a steep decay is found. As the number of spheres increases, this plateau ends at lower and lower altitudes, and the decay becomes slightly less steep.

This results in a probability of finding the spheres around z = 15d j that is respectively 0.23, 0.07, and 0.02 times the one of a single sphere for N s = 3, 6, and 12 (similar comparisons can be made for Re j = 8950 around z = 9d j , showing even smaller ratios as the collection of spheres does not reach high altitudes in that case).

While the radial position PDF are very different between both regimes for a single sphere, they have a similar shape for several spheres independently of N s and Re j , with a slow monotonous exponential decrease. Note that the left-most point, centered at r = 0, is probably underestimated due to measurement and binning biases. For a single sphere, half a bell curve centered around r = 0 is observed in the regular regime, while a much broader distribution is found for the chaotic regime, almost flat up to r = d j , followed by a steep decay. The exponential decay for several spheres is the steepest for N s = 3 at Re j = 8950 and the rate decreases with increasing both N s and Re j , so that 12 spheres in the chaotic regime have a probability of exploring r = 5d j that is only 4 times less than exploring r = 0.5d j (to be compared to a ratio of 90 for N s = 3 in the regular regime).

For completeness, Figure 8c-f) shows the power spectral densities (PSD)

of the vertical and radial positions in both regimes. The vertical position PSD show little change with respect to the number of spheres, with the exception of an increase of energy in the high-frequency range for Re j = 8950, which is reminiscent of the behavior observed for a single sphere as Re j increases (Fig. 4b)). This is the signature of the onset of a chaotic regime when several spheres are present in the flow configuration that yields a regular regime for a single sphere. This is also seen in Fig. 8d), where the high amplitude peak is lost for N s ≥ 3, as the spectra show a plateau followed by a power law decay. This decay becomes slightly less steep as N s increases so that the plateau is found in a lower and lower frequency range (to the point of not being observed for N s = 12, indicating that the dynamics remain correlated for longer duration). At a higher jet Reynolds number, the decay is even more moderate, and the behavior along N s is also similar, but in a lesser extent, and no plateau is observed over the acquisition duration explored.

The change in the particles' dynamics and in their sampling of the flow volume can be understood in the light of the velocity statistics. The mean velocity does not carry much information in a closed flow, as it should be zero if the spheres sample a large enough portion of the flow volume (by conservation of mass and as the mean particle velocity is independent of the particle characteristics [21; 22]). Despite a moderate spread, this holds true here (see Fig. 9a-b)). The track-by-track disparities around the mean value are a growing function of the jet Reynolds number, and also of the number of spheres for the r coordinate. The standard deviation of the spheres' velocity show less disparities, with the exception of the vertical coordinate at high

Re j and N s (Fig. 9c-d)). The standard deviation of the vertical velocity shows a trend that is very similar to the one of the vertical position (Fig. 7e)), with higher values in the chaotic regime and a peak for N s = 3 followed by a fast decrease. The slope is less steep at Re j = 11600, resulting in values of v ′ z almost three times larger than in the regular regime for N s = 12. The trend displayed by the radial coordinate is overall similar, but the slope at Re j = 8950 is much more moderate, and the chaotic regime presents a plateau. The values of v ′ r are hence larger in the regular regime for N s = 3 and 6 but become smaller than in the chaotic regime when 12 spheres are present. The high values of the radial velocity standard deviation are at the origin of the extreme events that lead to the wide radial exploration of the vessel, which leaves a clear signature with the monotonously growing radial position standard deviations (Fig. 7e)) or alternatively with the radial position PDF presenting wider and wider exponential tails (Fig. 8b)).

The velocity of the radial and azimuthal coordinates, r and θ, are approximately Gaussian and their shapes do not show a strong dependence on the number of spheres (not shown here for conciseness). Because the upward drag exerted by the jet on the sphere is strongly inhomogeneous in the cylinder, the competition between suspension (drag) and settling (buoyancy) strongly vary spatially and hence the vertical velocity distributions are asymmetric, except for the single sphere in the regular regime, where both forces remain in balance at all times (Fig. 9e-f)). While this asymmetry persists when the jet Reynolds number is increased, the distribution for N s = 12 resembles the ones at lower numbers of spheres (while a slight ordering of the distributions is still visible with respect to N s , especially for positive velocity values). The suspension for Re j = 8950 and N s = 12 being only partial, the vertical velocity PDF presents much narrower tails in that case, as the spheres do not get advected to high altitudes, resulting in both lower upward velocity increases during rising events and lower velocity decreases during falling events. This explains the strong differences in vertical position PDF observed in Fig. 8a), which is not necessarily observable in the of the average vertical position with N s , a monotonously decreasing function (Fig. 7a)) that results from the ordering of the vertical velocity PDF.

Residence times, collisions, and contacts

Velocity auto-correlation functions are used to probe the particle time dynamics, to complement high-frequency range of power spectral densities presented in Fig. 8c-f). As indicated the spectra, the short-time dynamics different behaviors along of N s for the different coordinates jet Reynolds numbers. This is illustrated by the radial velocity auto-correlation functions (Fig. 10a-b)), defined as R vr (τ ) = ⟨v r (t)v r (t + τ )⟩, where τ is a time lag and ⟨•⟩ is an ensemble average. Symbols marking 75% decorrelation are added to help illustrate the short-time decorrelation dynamics. For Re j = 8950, a short-time decorrelation, as well as oscillations, are observed for a single sphere, due to the oscillatory motion in the regular regime. The behaviors are almost indistinguishable for N s ≥ 3, as only a slight ordering towards longer decorrelation is visible as N s increases (much smaller than the increase in correlation time compared to the single sphere case). On the contrary, at Re j = 11600, the short-time decorrelation becomes faster as N s increases, starting from the single sphere case that shows the longest correlation time. The decrease of the correlation time seems monotonous but appears to saturate for N s = 12.

However, the auto-correlation change shapes, much than for Re j = 8950, so that a correlation time extracted from 100% decorrelation (i.e., first instance of R vr = 0), through the value of τ for this given threshold or through an integral up to this point, would give a much different result.

This means extracting a timescale of the particle dynamics would yield trends that heavily depend on the decorrelation threshold.

To quantitatively investigate the partial suspension of the spheres for N s = 12 and Re j = 8950, associated with narrow-tailed vertical velocity PDF (Fig. 9e)), the long-time dynamics is explored through a residence time analysis. For that purpose, two "regions" are defined, namely the ground and the rest of the flow volume. A loose definition of a grounded sphere is taken, to account for the slightly slanted bottom wall of the cylinder and for the possibility of having two spheres on top of each other, with z < d s (while a sphere sitting at the very bottom of the vessel presents z = d s /2). This definition is chosen for its simplicity and as it only biases towards slightly longer residence times on the ground (by only a couple times the temporal resolution as a sphere lands on or is suspended away from the ground), instead of missing grounded events by exclusion. The remainder of the time series is set to belong to the suspended sphere category. solid lines. No data is available for the single sphere case in the regular regimes as the sphere remains suspended as long as the jet is on, and only three grounded events are observed in the chaotic regime. The time intervals of the suspended events in between these grounds events are too long to be visible on the axes displayed here, except for one event, and the count for each time is reported instead of a PDF for the same reason. This shows two ground events of the shortest duration possible (dt = 1/60 s, the inverse of the acquisition frequency), corresponding to direct rebounds, separated by a suspension of 0.1 s, and a later ground event that lasts 10 times longer than the other rebounds (dt = 1/6 s), before the sphere gets lifted up by the jet.

When several spheres are present, many grounded and suspended sphere events occur. The most probable residence time, in both cases, remains the shortest time measurable. However, while for Re j = 8950 and N s ≤ 6 secondary peaks are present at long times (around 1.5 and 1/f z for N s = 3 and 6), the PDF decays from this value onward for N s = 12, showing the lack of suspension of the spheres in that case. This is accompanied by a broad secondary peak at long times for grounded events (1.5 ≤ dtf z ≤ 4.5) (while no such peak is visible for lower numbers of spheres, as only a small peak at 1.5/f z appears for N s = 6). A peak for dt/f z = 5 and at a probability of 0.15 % also exists but is not displayed for visibility. For Re j = 11600, while the suspension probability is smaller at long times and higher at small times for N s = 12 compared to 3 and 6, no sphere remain grounded for more than dt/f z = 2.5. Figure 6e) seemed to imply that longer residence times on the ground could be found for N s = 3 than N s = 12. However, this was caused by many rebound events that translate to higher probabilities of short-time grounded events (that can happen subsequently) compared to the dynamics at N s = 12, while this latter case presents as expected much higher probabilities of long grounded events.

The overall changes in suspension capacity are measured in Fig. 10ef), showing respectively the average number of suspended spheres and the percentage of the time where N spheres are suspended. The average number of suspended spheres, normalized by N s , decreases almost linearly with N s , with a steeper slope at the lower jet Reynolds number value. In addition, the dashed lines report the maximum number of spheres suspended, showing that all the spheres can be suspended up to respectively N s = 3 and 6 for Re j = 8950 and 11600, while only approximately 60 and 80% of the spheres can be suspended at once when 12 spheres are present. To complement this, the dashed lines in Fig. 10f) report the percentage of the measurement duration that presents a suspension of at least half of the spheres, a value almost constant and always above 90% for the high jet Reynolds value for N s ≤ 6, while a linear decay is seen for the lower Re j value over the whole range of N s . The suspension percentage of all the spheres is on the contrary a very quickly decaying function, already lower than 5% at N s = 6 independently of the jet Reynolds number. An alternate probing of the suspension capacity of the jet is done by studying the statistic of the cumulative altitude, defined as Z c = Ns n=1 z n with z n the altitude of each sphere. The cumulative altitude presents statistics that are close to a Normal distribution, with a standard deviation that mostly depends on the jet velocity, while the average value increases both with N s and Re j (not shown here for conciseness). For both Reynolds numbers, the increase of Z c with N s is approximately linear and equal to more than twice the value of ⟨z⟩ in the case of a single sphere. This means that the total potential energy of the collection of spheres increases with N s , which is thought to be due to a better harnessing of the jet's kinetic energy by the spheres as they experience an upward drag force surpassing their buoyancy a region of the flow that is larger than the vicinity of the equilibrium position.

In addition to grounded sphere events, the collective dynamics change the inter-particle distances as N s increases. At each frame, the distance in 3D between pairs among the N s spheres, dr, is computed, and the probability density functions are reported in Fig. 11a-b). The probability of dr = d s , i.e., two spheres touching, increases drastically with N s , only after N s = 6 for Re j = 8950 and monotonously for Re j = 11600. For long distances, the probability is larger at N s = 3 in both cases, but is followed by N s = 12 at low Reynolds number instead of the expected N s = 6 found for the higher value of Re j . This can be explained by the fact that 12 spheres occupy each other vicinity in a more compact fashion, so several spheres in contact are separated by distances greater than dr = d s . This is well illustrated by the peak of probability at dr = 2d s , indicating three spheres aligned and touching (e.g., A is in contact with B which is in contact with C). As more suspension is found for that case at Re j = 11600, this secondary peak is very broad and almost merges with the peak at dr = d s . In addition, for this jet Reynolds number, the PDF is almost flat at N s = 3, indicating an equiprobability up to distances of 6d s that corresponds to a dynamics almost void of particle-to-particle interactions.

Such interactions are explored by analyzing the probability of spheres touching each other, divided into two categories: contacts and collisions.

Contacts are defined as events happening for grounded spheres, where low velocity values are expected, while collisions occur between two suspended spheres, often associated with large momentum transfers (e.g., collision between the black and brown colored spheres in Fig 6a) resulting in a high altitude reached by the black sphere while the brown one becomes grounded).

The numbers of contacts and collisions occurring at a given time, N c , are computed and their PDF are reported in Fig. 11c-d). Contacts are found to be always less likely (or as likely for N s = 12 at Re j = 11600) than collisions, with the exception of N s = 12 at Re j = 8950. In this case, there is a high probability of having several spheres on the ground, in contact with each other (showing almost an equiprobability up to 3 contacts before a decay)

and hence unlikely to be easily resuspended. This explains the occurrence of very long residence times with non-negligible probabilities. The PDF of the number of simultaneous contacts or collisions N c are approximately exponentially decreasing, and the slopes are decreasing functions of N s .

On average, the inter-particle distance is constant along N s in the regular regime, at a low value (⟨dr/d s ⟩ ≃ 2), and is found higher and only slightly decreasing in the chaotic regime (between 4.5 and 3 d s , Fig. 11e)). However, the numbers of contacts and collisions are linearly increasing with N s . The rate of increase is approximately the same for collisions at both jet Reynolds numbers (with collisions being always about twice more frequent than at the lower Re j value). On the contrary, the slope along N s is much bigger for contacts at a low Reynolds number than at the higher value, so the frequencies are about equal for N s = 3, but three times larger when 12 spheres are present. In this case, contacts are 1.45 times more likely than collisions, while they are as likely for Re j = 11600.

Discussion

The increase of probability in larger (negative or positive) vertical velocity values compared to the single sphere case deserves to be discussed further (Fig. 9e-f)). For the regular regime, a very narrow PDF is observed for v z as the sphere only oscillates at small amplitudes. The typical velocity can be estimated as z ′ , the oscillation amplitude, times f z , its frequency, resulting

in v O z /U j = z ′ f z /U j ≃ ±0.
02. This value is in very good agreement with the largest (absolute) values displayed in Fig. 9e), up to 10 times smaller than the ones found when more than one sphere is present. For N s ≥ 3, a large decrease of the zero vertical velocity probability is found, and both large negative and positive values become much more probable. Collisions lead to fast and drastic upward velocity increase, leading to very high altitudes (Fig. 6a)), that are followed by falling events with slow growth of the downward velocity. This difference in timescale between both events explains partially the wider negative tail on the velocity PDF, compared to the positive one.

In addition, this is also due to the fact that the jet, where upward fluid drag is imparted on the spheres, occupies only a small portion of the flow volume 2 . At this low jet Reynolds number value, spheres can spend longer periods of time on the bottom wall of the vessel for N s ≥ 6 (Fig. 10c)), so 2 A cylinder of height 40 mm and radius 30 mm is 1.1 10 -4 m 3 , while a cone of the same height and of full angle of 22 • , based on twice the spreading rate of the velocity half-width radius, see [START_REF] Pope | Turbulent flows[END_REF], only contains 2.7 10 -6 m 3 , i.e., a bit less than 2% of the cylinder's volume.

that the zero vertical velocity probability increases drastically. The dynamics become largely modified by the many contacts at N s = 12 (Fig. 11), 50% more likely than collisions in the same flow regime, and 3 times more likely than contacts at Re j = 11600. This explains the much narrower tails on the vertical velocity PDF in this case (Fig. 9e)). The behavior along N s in the chaotic regime observed in Fig. 9f) shares some similarities, with broad tails for N s ≥ 3 and a lower zero vertical velocity probability compared to the single sphere case, but the conditions at N s = 12 do not stand out as all the spheres remain moderately well suspended (Fig. 10e-f)). The PDF are all ordered along N s , with a narrowing of the tails and an increase of the zero velocity probability. Note that due to the chaotic dynamics at Re j = 11600, the single-sphere vertical velocity PDF already displays moderate tails, wider for negative values, similar to the case of several spheres (while in a less pronounced fashion).

The studied system, unlike the case of a single sphere, is multistable, with each sphere evolving around a metastable position, until it switches for another one, through large fluctuations in the fluid's forces or by a collision with another sphere. In addition, when more and more spheres are present, contacts increase the stability of positions on the ground, in the sense that it requires events of even larger magnitudes to move the spheres.

The global metrics report on the sum of the metastable positions, so this framework retains the system axisymmetry (where the word global indicates obtained through an ensemble average here). However, some deviations are observed, as on the mean velocity where a zero value is expected for a closed flow (Fig. 9a-b). This is not due to a lack of statistical convergence, but rather to the fact that only a few tens of metastable positions are summed, and even a slight dissymmetry in the number of such events can result in deviations from axisymmetry. The difficulty in retaining symmetries in the experimental measurements of such multistable systems is a well documented challenge, that is very common in the general contexts of physics, chemistry, or biology [24; 25], and has been many times reported in fluid mechanics with fluctuation-driven bi-stable systems (e.g., [26; 27; 28; 29; 30]) as well as other particle-laden flows [31; 32; 33; 34; 22].

One objective of this study is to explore the competition between several spheres for a given equilibrium position in the jet and to characterize its influence on the suspension regimes. While both regular and chaotic regimes are explored for a single sphere in this configuration, the dynamics are always chaotic as soon as 3 spheres are present. However, drastic changes in behavior were highlighted by change in the jet Reynolds number and the number of spheres, requiring the definition of other suspension regimes, that adequately describe the collective suspension dynamics. Since the main mechanism modifying the dynamics is the interaction between spheres, this definition should be based on collisions and contacts. The change in the exploration of the flow volume, the velocity statistics, and the residence times all appear to be linked to the occurrence of these events.

When almost no contact or collision exists, the spheres evolve independently of each other, by definition in an uncoupled way, as the probability of spheres being far from each other is large (Fig. 11b). This is only found for N s = 3 and Re j = 11600 and the dynamics remain chaotic, as even a single sphere shows extreme events here (due to hydrodynamics fluctuations). This defines an "uncoupled chaotic" regime. When the frequency of contacts, divided by the number of spheres, becomes close to or larger than the oscillation frequency, this onsets a "collision-driven" regime. It occurs for the lower jet Reynolds number with 3 and 6 spheres and for the higher Reynolds number at 6 and 12 spheres. In this regime, collisions occur as much or more than contacts, and the suspension of half the spheres occurs more than 50% of the time (Fig. 10f). In the case of respectively N s = 3 and 6 for Re j = 8950 and 11600, this rate is in the vicinity of 90% and all the spheres can be suspended at times, defining a high-suspension sub-regime, while the other cases are only moderately suspended. Finally, when 12 spheres are present at the lower jet Reynolds number value, contacts are more likely than collisions, and long residence times of grounded spheres are not rare, as groups of spheres occupy the bottom wall of the cylinder and only a quarter of the spheres are suspended in average, and even suspending only half the spheres is very rare (Fig. 10e-f). This defines a "contact-driven" regime, and the various regimes are summarized in Table 2. The uncoupled chaotic and contact-driven regimes are illustrated in the Supplementary Materials (videos).

For the single-sphere case, the occurrence of a chaotic regime at higher jet velocities explored by [START_REF] Davoust | The suspension of a sphere in a turbulent jet[END_REF] is generalized by the current study. The exploration of the parameter space differing widely, the onset of this regime cannot be presented along a unified non-dimensional-group (e.g., the Froude number as proposed by [START_REF] Davoust | The suspension of a sphere in a turbulent jet[END_REF]), and a parametric study would be required to answer this question. The orders of magnitude of the mean altitudes in both the regular and chaotic regimes are however in good agreement. When a collection of spheres is suspended in the jet, the drastic changes in dynamics lead to more time spent on the ground, both with more frequent and longer events, in opposition to only quick rebounds in the single-sphere case. This means that the geometry of the studied system can play an important role (comparatively to the role played by the jet), in particular for the dynamics of grounded events. For instance, changes in the current geometry (e.g., value of the slope of the cylinder base, change in the cylinder's diameter) could result in differences in the values reported here such as the distributions of the ground residence times or the contact frequencies. In addition, when many spheres are present, such as in the contact-driven regime, the size of the container becomes only slightly larger than the pile formed by the spheres, and confinement effects play a role as well. However, beyond the exact values of the metrics and their scaling laws, which are prone to such biases and are unlikely to be directly transposable to other systems, the qualitative picture is presumably very general. In the case of a single sphere, the changes of regimes are explained by a competition between the apparent size of the jet at a given altitude and the sphere diameter, and such confinement balance is solely caused by the base flow geometry and therefore independent of the surrounding system geometry. Similarly, the transition to the collision-driven regime is not likely to be influenced by the system geometry, as it is governed by how many spheres can be sustained in the vicinity of the stable region, a function of the sphere diameter and the jet apparent size (at a given altitude, which adds a dependence on the fluid's and spheres' densities). The threshold of the onset of the contact-driven regime is probably heavily influenced by the geometry, but its existence is most likely universal. In addition, the experimental methodology and the analysis framework presented here are general and can be used to study such systems in various variations of the parameters space and geometry, as well as broader aspects of particle-laden flows.

Conclusion

X-ray radiography was used with a divergent X-ray source and a planar detector to record time series of attenuation maps resulting from a collection of mono-dispersed glass spheres immersed in water. With the knowledge of the system geometry and the spheres' diameter, radioSphere was leveraged to retrieve the 3D trajectory of each sphere, yielding the 4D kinematics of the collection of spheres. This stands as the first implementation of radioSphere in a dynamical system, as its use had been limited to synthetic data and a static granular assembly [START_REF] Andò | Single-projection reconstruction technique for positioning monodisperse spheres in 3d with a divergent x-ray beam[END_REF]. The combination of the X-ray measurement technique and the developed analysis method allowed for the investigation of the suspension of a collection of heavy spheres in an upward vertical jet. Such configuration is well described for a single sphere, and the use of two different jet velocities led to the exploration of two suspension regimes, namely the regular and chaotic regimes [START_REF] Davoust | The suspension of a sphere in a turbulent jet[END_REF]. In both cases, the sphere has an equilibrium position at a certain distance from the jet and along the jet's axis, but can present either regular oscillations or large deviations associated with extreme events. The number of spheres was varied in the range 3 ≤ N s ≤ 12 at each jet Reynolds number value, exploring 6 conditions for a collection of spheres, to investigate the competition of many spheres for an equilibrium position.

While the trajectories could have been retrieved using multiple cameras in direct visible-light imaging in the case of a few spheres, the conditions explored led to many collisions and prolonged contacts between several spheres which necessitated the use of X-rays. In addition, a cylindrical vessel, with a slightly downward slanted base, contained the working fluid and maintained a constant number of spheres throughout the measurements. While the specifics of the geometry (strong confinement and slanted base in particular) may prevent direct transposition of the results to other systems, the qualitative picture is presumably very general, and that the methodology to characterize the different regimes is also of broad applicability.

The exploration of the flow volume by the spheres was investigated using the statistics of the spheres' positions, namely with the first two moments and the probability density functions, while the dynamics were presented in terms of power spectral densities and auto-correlation functions. Wide changes occur as N s increases, which can be explained by the spheres' veloc-ity statistics. In particular, the vertical velocity probability density functions present broad tails that narrow with N s , resulting in a decrease in the altitudes reached and an increase in the radial exploration. The tails originate from collisions between spheres, yielding large positive and negative vertical velocity values. As more spheres are present in the vessel, some spheres exit completely the suspension region created by the jet for some time, and the probability of finding a sphere at the bottom wall of the vessel increases.

Spheres can be found in contact with each other there, decreasing the likelihood of a resuspension. The suspension rate and the occurrence of collisions and contacts led to the introduction of three suspension regimes of a collection of spheres, listed in decreasing degree of agitation: "uncoupled chaotic" regime, "collision-driven" regime (split into a high or moderate level of suspension), and "contact-driven" regime. As only two parameters were varied in the current study, extending this framework and the existence of these suspension regimes in a wider parameter space would pose a natural future direction. In particular, variation of the sphere-to-jet diameter and density ratios are possible with the current approach, but they would require an adaptation of the imaging configuration, namely of the distances between the source, object, and detector, as higher altitudes would be reached and the apparent magnification range would change.

Appendix: Calibration of raw attenuation into path length images

An essential initial step of radioSphere is the transformation of the raw attenuation radiographic projections into calibrated path length images (expressed in mm). For the presented experiments, this was done along a twostep approach. First, a set of 100 radiographies of a single sphere centered on the detector panel were averaged, and the relationship between the path length inside the sphere and the measured attenuation on the detector panel was calculated, as shown in Fig. .12a). This calibration curve for a single sphere follows the Beer-Lambert law. However, for an assembly of spheres the linearity of the simple Beer-Lambert attenuation law is violated, mainly due to beam-hardening and scattering artifacts [35; 36].

To account for a non-linear attenuation law beyond one particle diameter, a second calibration step is performed. For this, the stationary projection In summary, for each experimental configuration, a non-linear calibration curve is estimated and used to transform the set of acquired radiographs to path length images. This law corresponds to the studied spheres material, specific acquisition settings, and imaging geometry, accounting at the same time for any non-linearity caused by beam-hardening and scattering artifacts.

It should be noted here that alternatively a more automatic parametric function fit between the scatter plot of the measured detector attenuation vs.

path length for each pixel can be performed through a least-squares method in the style of [START_REF] Fragnaud | Cad-based x-ray ct calibration and error compensation[END_REF].

Figure 1 :

 1 Figure 1: (color online) Schematic representation of the set-up inside the X-ray scanner.

  2) and their associated nondimensional parameters. Re j = U j d j /ν l is the jet Reynolds number based on the jet diameter d j and the mean exit velocity U j . Re s = Re j d j / ⟨z⟩ is a lower bound of the sphere Reynolds number, where the slip velocity is estimated as the fluid velocity at the equilibrium position of a single sphere ⟨z⟩ through a linear decrease starting at z = 0 (ignoring the potential core). A higher bound is approximated by DRe j ≃ 2Re j which applies to the at-rest sphere meeting the jet. V = U s /U j is the sphere-to-jet velocity ratio, with U s = 4/3d s /C d g(ρ sρ l )/ρ l terminal velocity based on free-fall in quiescent water, where the drag coefficient C d is approximated to be 0.5, and g is the gravitational acceleration. The Froude number, comparing the free-fall timescale of the sphere to the jet's timescale is F r = d s U j /d j U s = D/V . The fluid is distilled water, at an ambient temperature of 25 • C, with a kinematic viscosity of ν l = 8.96 10 -7 m 2 s -1 and density of ρ g = 1.18 kg m -3 . The soda lime glass spheres have a density of ρ g = 2.5 kg m -3 .

Figure 2 :

 2 Figure 2: (color online) Schematic representation of the radioSphere technique. The input projection is transformed into a calibrated path length image, which is then deconvoluted with a series of structuring elements. The estimated projection is then refined through a residual minimization iterative algorithm. The x and y axes on insets c) and e) correspond to the frequencies per image length, centered on the middle of the image, which means the axes span from [-L X /2, L X /2 -1] and [-L Y /2, L Y /2 -1], where L X and L Y are the size of the image, in pixels, in the x and y directions. Since the recorded images are not on a square detector the rings created in Fourier space are distorted by this aspect ratio.

Figure 3 :

 3 Figure 3: (color online) a) Measured radiography (normalised by the flat-field), b) calibrated path length image, (c) residual after the first run of optimization, d) residual after the second run of optimization accounting for motion blur.

2 Figure 4 :

 24 Figure 4: (color online) a) 2D projections of the 3D trajectory of a single sphere at a jet Reynolds number of Re j = 8950 (•) and 11600 (□) respectively (alternatively F r = 6.3 and 8.16). The dark to light colors indicate the trajectory evolution with time, from start to end. b) Corresponding power spectral densities (PSD) for a cylindrical coordinate system, normalized by the variance (square of the standard deviation i ′ , where i = r, z).

Fig. 5 and 6 when

 6 Fig. 5 and 6 when Re j is increased. In addition, Figure 5 uses a cartesian frame of reference to represent the trajectories.

FigureFigure 5 :

 5 Figure 5b-c), while only showing a subset of the 4D kinematics, readily allows for observations of the change in dynamics as the number of spheres N s increases at a fixed jet Reynolds number. For Re j = 8950, the most

Figure 6 :

 6 Figure 6: (color online) Time series of the vertical position normalized by the jet diameter, z/d j , along the non-dimensional time tf z , where f z = 2.1 Hz is the vertical oscillation frequency in the regular regime (Re j = 8950). Each sphere trajectory is represented as the number of spheres N s is increased along the row and for both jet Reynolds numbers considered (columns; alternatively F r = 6.3 and 8.16). A sphere laying on the base of the cylinder is at z/d j ≃ 1.1 and the dashed and dash-dotted lines respectively represent the mean and maximum altitude reached in the case of a single sphere suspended in the jet.

Figure 7 :

 7 Figure 7: (color online) Average vertical z (a), transverse y (b), and radial r (d) positions normalized by the jet diameter d j as a function of the number of spheres N s . c) is a normalization of a) by the square of the ratio of the sphere terminal velocity U s and the jet mean exit velocity U j . Similarly, e) and f) display the standard deviations of the vertical and radial positions, z ′ and r ′ . Once a mean or standard deviation is obtained from each individual sphere trajectory, it is averaged to produce the solid lines, while the dashed lines are obtained from the minimum and maximum of these sets of values. The symbols indicate the jet Reynolds number values: Re j = 8950 (•) and 11600 (■) respectively (alternatively F r = 6.3 and 8.16).

Figure 8 :

 8 Figure 8: (color online) Probability density functions (PDF) of the vertical z (a) and radial r (b) positions normalized by the jet diameter d j at different number of spheres N s and jet Reynolds number Re j . c-f) Power spectral densities (PSD) of the vertical z (left) and radial r (right) positions, normalized by the corresponding variances, at different number of spheres N s , and for jet Reynolds numbers Re j = 8950 (middle) and 11600 (bottom), alternatively Froude numbers F r = 6.3 and 8.16.

Figure 9 :

 9 Figure 9: (color online) Mean vertical (a) and radial (b) velocities normalized by the jet mean exit velocity U j as a function of the number of spheres N s , and similarly with the standard deviation of the vertical (c) and radial (d) velocities. The mean and standard deviations values have respectively been multiplied by 10 2 and 10 3 for visibility. Once a mean or standard deviation is obtained from each individual sphere trajectory, it is averaged to produce the solid lines, while the dashed lines are obtained from the minimum and maximum of these sets of values. Probability density functions (PDF) of the normalized vertical velocity for jet Reynolds numbers Re j = 8950 (e) and 11600 (f), alternatively Froude numbers F r = 6.3 and 8.16.

FigureFigure 10 :

 10 Figure 10c-d) reports the probability density functions of consecutive periods of time spent on the ground or suspended, respectively in dashed or

Figure 11 :

 11 Figure 11: (color online) Probability density function (PDF) of the distance r between pairs of spheres, normalized by the sphere diameter d s for Re j = 8950 (a) and 11600 (b) (alternatively F r = 6.3 and 8.16). PDF of the number of collisions (solid lines) and contacts (dashed line) happening simultaneously for Re j = 8950 (c) and 11600 (d). e) Mean pair separation as a function of the number of spheres N s . f) Collision (solid lines) and contact (dashed lines) frequency, normalized by the vertical oscillation frequency in the regular regime f z = 2.1 Hz and N s , as a function of N s . Frequencies are obtained as the number of events divided by the measurement duration.

Figure . 12 :

 12 Figure .12: Fitted calibration of normalized attenuation versus path length for a single 10 mm sphere (a) and 12 stationary spheres at lowest input velocity (b).

  of each experimental configuration is transformed into a path length image (expressed in mm) based on the single sphere linear attenuation fit obtained previously. radioSphere is then run to reconstruct the 3D position of the stationary spheres, which are, in turn, used to produce a synthetic projection of the spheres assembly, expressed in mm. The synthetic and the acquired projections of the stationary spheres are then used to fit a new attenuation law, as shown in Fig..12b).

  [START_REF] Feng | The motion of a solid sphere suspended by a newtonian or viscoelastic jet[END_REF] mm for the highest. The latter requires a larger Regime U j (m/s) Re j

			Re s	V	F r
	Regular	1.8	8950 2560 0.35 6.30
	Chaotic	2.3	11600 1780 0.27 8.16

Table 1 :

 1 Suspension regimes explored (see Section 3.

Table 2 :

 2 Collective suspension regimes as functions of the jet Reynolds number Re j = U j d j /ν l and the number of spheres N s : the regular and chaotic regime are defined only for a single sphere, while several spheres can be in an "uncoupled chaotic" regime, a "collisiondriven" regime (associated with a high suspension percentage, marked here by the

*

, or with moderate suspension), and a "contact-driven" regime. The regimes are termed and distinguished based on the quantitative measurements of suspension rates, contacts and collisions from Fig.

10 and 11

The structuring element series can be either synthetic (as used here) or experimental by acquiring a set of radiographs of a single centered sphere while varying its position along the X-ray beam.

Acknowledgments

This work has been supported by a LabEx Tec21 grant (Investissements d'Avenir -Grant Agreement # ANR-11-LABX-0030).