# Viscous rebound of a quasi-2D cylinder on a solid wall

Alicia Aguilar-Corona1 ,Micheline Abbas2,∗,Matthieu J. Mercier3, Laurent
Lacaze3.

1Facultad de Ingenieria Mecanica, Universidad Michoacana de San Nicolas de
Hidalgo,
CP 58030, Morelia, Michoacan, Mexico
2Laboratoire de Génie Chimique, Université de Toulouse, CNRS, INPT, UPS,
Toulouse, France
3Institut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse,
CNRS, Toulouse, France

<EMAIL_ADDRESS>
##  1 Abstract

The purpose of the present study is to extend the simple concept of apparent
coefficient of restitution, widely approached in the literature for the case
of a single-contact-point between a sphere and a wall, to the case of bouncing
whose complexity is increased due to the shape of the contacting object. For
this purpose, experiments are carried out with a finite-length cylinder,
freely falling in a liquid at rest. A rigid tail is attached to the cylinder,
allowing to maintain a vertical trajectory and to keep the axis of the
cylinder parallel to the bottom wall down to small gap between them. Yet, more
complex 3D motions of the cylinder with respect to the wall occur during
bouncing, including multi-contact-points between the cylinder and the bottom
as well as cavitation. When the Stokes number $St$ (the ratio between
characteristic inertial forces experienced by the particle compared to viscous
forces in the fluid) is increased, the experimental results suggest that the
ratio of the apparent coefficient of restitution to the solid one increases
from $0$ (at low $St$) to 1 (at large $St$) with a critical Stokes number
$St_{c}$ below which no bouncing is observed, as usually obtained in the
literature for a single-contact-point bouncing sphere. In a conceptual
approach of the influence of a cut-off length scale prior contact and a
contact time scale on the observed experimental results, we investigated
numerical modelling of an idealized situation (2D infinite cylinder falling
parallel to the wall). To this end, we carried numerical 2D simulations where
the fluid equations of motion were coupled to the particle equation of motion
through an Immersed Boundary Method. The particle equation of motion was
coupled to an elastic force to model bouncing. This numerical model requires
parameterization of the cut-off length (interpreted as a roughness) and the
contact time (associated with the contact elasticity) used here to capture the
experimental observations. The simulations confirmed that i) the departure of
the coefficient of restitution from 0 is strictly dependent on the apparent
roughness and ii) the coefficient of restitution depends on the contact time.
Finally, in an effort to rationalize the experiments and the simulations on
such conceptual approach, we model the coefficient of restitution as the
product of two contributions to the mechanical loss of energy: the collision-
to-terminal velocity ratio ($V_{c}/V_{t}$) of the approach-phase and the
rebound-to-collision velocity ratio ($-V_{r}/V_{c}$) of the contact-phase. We
then interpreted the experimental measurements in light of this model, showing
evidence that the assumption of a global relationship between the contact time
and an apparent roughness (all being linked to the bouncing complexity
including multi-contact-points and cavitation in experiments) leads to a
reasonably good prediction of the coefficient of restitution in the
intermediate regime in $St$. This suggests the relevance of lumping the
complex details of physical phenomena involved during contact into a simple
concept based on the contact apparent roughness and elasticity.

##  2 Introduction

The interaction of solid particles evolving in a viscous fluid has been
considered in many studies for its obvious relevance in many geophysical and
industrial applications involving multiphase flows. The related key questions
are numerous ranging from the mesoscale dynamics, that is the collective
motion of the suspended particles, to the microscale dynamics at the scale of
the grain. In particular, if one simplifies such complex systems by
considering only two approaching particles, one handles local processes such
as dissipation induced by their interaction. This includes the viscous
dissipation due to fluid motion induced by their relative motion close to
contact as well as the mechanical dissipation during solid contact.

A main issue, when dealing with such apparently simple systems, is the
antagonism between geometrical and mechanical properties of modelled particles
and the ones of real particles. Then, two main aspects are usually modelled,
as opposed to solved, which are (i) the local deformation due to elasticity of
the particles and (ii) the surface roughness. Otherwise, infinitely rigid and
smooth objects would lead to a complete dissipation of the initial kinetic
energy, and no bouncing could occur. Such dissipation would be due to the
lubrication force induced by the viscous flow in the gap between particles
diverging when the gap goes to zero for perfectly smooth and rigid particles
as evidenced by [1] in the case of a sphere approaching a wall. However, a
full dissipation preventing bouncing has been shown through several
experimental studies to be unlikely, at least above some given threshold of
inertia close to contact [2, 3, 4, 5, 6, 7]. Therefore, the singularity of the
lubrication force at zero gap, that would only lead to solid contact on a
infinite time scale with zero relative velocity, should somehow be
regularised. This is why modelling the bouncing of a sphere is often sought as
a regularization of the physical system.

The elastohydrodynamic collision model proposed by [8] and [9] is a
regularization approach incorporating (i), which thus allows bouncing while
preventing solid contact. This led to several experimental studies dealing
with the bouncing of spherical particle on a rigid wall, for which the surface
roughness of particle could not be disregarded [2]. Then, on the opposite
side, the possibility of rigid grain to bounce only due to surface roughness
has been addressed in the literature [5]. As suggested by [10], (ii) proving
solid contact required sharp shape roughness, otherwise the elastohydrodynamic
collision is only scaled down to the roughness scale. Trying to separate (i)
from (ii), [6] has experimentally highlighted the possibility of solid contact
between rough rigid object on a finite time scale. An extension of the
elastohydrodynamic collision including (ii) could therefore be of relevance
[11]. Note that additional aspects can be also found in the literature, as for
instance the liquid-glass transition in the presence of highly viscous liquid
films between colliding surfaces [12, 13].

From a numerical modelling point of view, the scale of regularization, either
due to (i) or (ii), is usually too small to be reasonably resolved, in
particular for complex systems involving several particles. Actually,
significant advances have been recently made on the numerical modelling of
complex shape and non-rigid particles [see for instance 14]. However, these
methods still suffer from a limitation due to the computational cost and the
spatial scale of shape complexity, which usually lead to locally smooth
surface. Accounting for the shape complexity of a body does not necessarily
allow to capture the required change of scale of surface roughness that would
be at the origin of the regularization at contact mentioned previously. Then,
a specific attention is still required to incorporate the small-scale physical
processes as closure model into classical numerical approach such as Discrete
Element Method (DEM), for which (i) and (ii) are not fully resolved. This
approach has been extensively proposed in the literature using a coupled DEM-
DNS solver for the solid-fluid resolutions [15, 11, 16, 17, 18, 19, 20, 21,
22]. All these approaches are somehow similar, as the meshgrid of the fluid
solver covers the entire domain including the solid object to prevent from
complex and costly numerical algorithm when dealing with body-fitted grid
methods, which would moreover be incompatible with an objective multi-body
systems. For these fixed-grid methods, an extra term is required to enforce
the solution within solid body. Methods can slightly differ but mostly lead to
the same conceptual methodology. In this case, the above mentioned
singularity, or lack of resolution, occurs at the meshgrid of the fluid
solver, when the two solid surfaces close to contact reach a distance of the
order of the mesh-size. Due to numerical limitation, this length scale can
actually be quite important compared to the typical diameter of the particles
$D$, i.e. only one order of magnitude smaller $D/10$. As this scale is often
thought as large compared with the roughness of real grains, a subgrid
lubrication model is often added when the solids get closer. Moreover, a
separation between solid contact time scale and fluid time scale is often
considered, assuming strict rigidity of the solids. These two aspects
(lubrication model and rigidity) still require some attention. In particular,
adding a lubrication model allows to delay the singularity to a smaller scale,
and therefore to obtain a better quantitative agreement with experiments [as
the fluid dissipation is increased to a level closer to real situation, see
20, for instance]. However, this remains quite an empirical alternative. Using
Immersed Boundary Method (IBM) on a fixed grid, [21] attempted to refine the
meshgrid in the gap between solid object (sphere/wall) down to $\sim D/1000$
using local refinement and taking advantage of geometrical symmetry. Doing so,
a subscale lubrication model is not added, allowing to capture most of the
dissipation as the mesh reaches a length scale comparable to real roughness.
Then, roughness is seen as an apparent roughness based on the meshgrid [see
supplementary material in 23]. Even if the link between numerical dissipation
and real dissipation remains unclear, the scale issue requiring an additional
lubrication model is at least removed with such approach.

Quite surprisingly, whatever the ‘philosophy’ of regularization or physical
process at the microscopic scale, a robust observation obtained from
experiments and modelling either numerical or theoretical, is the evolution of
an effective coefficient of restitution including both viscous dissipation and
mechanical/structural elasticity, with the dimensionless Stokes number $St$,
measuring the ratio of particle inertia to fluid dissipation as explicitly
defined in eq. (2). In particular, a critical Stokes number $St_{c}$ is
obtained delineating a viscous region for which no bouncing is observed
($St<St_{c}$) and an inertial region where bouncing occurs ($St>St_{c}$). This
robustness allowed numerical simulations, as discussed later on, to be
relevant to predict this the transition from viscous to inertial regimes, even
if the local mechanisms allowing this bouncing (i vs ii) remain unclear. Then,
it would suggest that any complex mechanism at play during bouncing, due
either to the local complex shape of the bouncing particle or to deformation
and/or to phase change (in the liquid) induced by strong local pressure just
prior contact and/or low local pressure upon motion reversal, respectively,
can be captured by simplified modelling based on effective bouncing
characteristics. The terms elasticity and roughness, classically used in the
literature of IBM/DEM modelling, thus refer to the adjustable bouncing
parameters which are actually the time scale of contact and the regularization
length prior to contact, to be parameterized in an IBM/DEM bouncing
simulation. Then even if discriminating and modelling elasticity versus
roughness remain uncertain, their conceptual consequences on the effective
bouncing are promising and still deserve attention. One proposes to adopt this
approach to extend the bouncing models developed for the case of a sphere [20]
towards a conceptual approach of apparent elasticity and apparent roughness as
a generalization of more complex bouncing configurations.

For this purpose, experiments (section 3) and simulations (section 4) are used
to investigate the coefficient of restitution $e$ of a cylinder bouncing on a
‘smooth’ wall, as a function of its inertia characterized by the Stokes
number. Experiments are designed to provide the bouncing of a finite-length
cylinder onto a horizontal surface while simulations are designed to provide a
twin 2D situation of a 2D cylinder bouncing along a horizontal line, in which
roughness and elasticity are therefore subscale models which could mimic the
complexity of finite-length 3D bouncing.

In the case of an object falling under gravity in a viscous fluid, the Stokes
number is intimately related to the settling Reynolds number, defined here as

$Re=\frac{\rho_{f}DV_{t}}{\mu},$ (1)

The Stokes number corresponding to a falling object is built from the ratio
between the object relaxation time and the fluid motion at the scale of the
object. Instead of considering different expressions associated with the
cylinder settling in the experiments and numerical simulations (since the
shape is not exactly the same), we use for simplicity the commonly used
expression for a falling sphere:

$St=\frac{1}{9}\frac{\rho^{*}DV_{t}}{\mu}=\frac{1}{9}\left(1+\frac{\rho_{f}}{\rho_{p}}C_{M}\right)\frac{\rho_{p}}{\rho_{f}}Re.$
(2)

Here $\rho_{f}$ and $\mu$ correspond to the fluid density and viscosity.
$D=2R$ indicates the cylinder diameter, $R$ its radius.
$\rho^{*}=\rho_{p}\left(1+\frac{\rho_{f}}{\rho_{p}}C_{M}\right)$ is the
effective cylinder density taking into account the added mass effect and
$C_{M}$ represents the added mass coefficient. $V_{t}$ corresponds to the
cylinder terminal settling velocity, reached when the drag force balances the
cylinder apparent weight. Beyond the $St$, the bouncing model will then be
linked to a dimensionless equivalent elasticity parameter $\alpha$ and a
dimensionless apparent roughness parameter $\eta$, which will be defined
sections 4 and 4.2.

The paper is organized as follows. First, the results of experiments carried
with a cylinder of a circular base are reported in section 3. Despite the
setup apparent simplicity, it is shown that the values of the coefficient of
restitution as a function of the Stokes number are not easily rationalized
onto a master curve, but a shift remains between the reported curves. In order
to interpret the data, we postulate that i) the scatter is mainly associated
with the cylinder pitching at the onset of the collision process and the
generation of cavitation bublle at contact, and that ii) the physical
phenomena at play during the contact can be interpreted in terms of a cut-off
length at contact and a contact time, i.e. a contact roughness $\eta$ and a
contact elasticity parameter $\alpha$, as explained above. Numerical
simulations of a bouncing 2D cylinder are then provided in section 4 and
discussed according to a simple theoretical model. As the configuration
remains 2D all along the particle trajectory (collision along a contact line),
$\eta$ appears to influence mainly the value of the critical Stokes number,
while $\alpha$ is shown to lead to different levels of the $e-St$ curves. The
paper ends with a discussion section 5, where we show that using the
numerically and experimentally measured contact time, we capture the different
levels of $e-St$ curves in a reasonable way. This confirms the importance of
measuring contact time along the process of characterization of the complex
collision process.

##  3 3D bouncing: experiments

###  3.1 Experimental setup

To model an idealized rebound in a laboratory setup, we consider the
controlled settling of a finite width cylinder of diameter $D$ and length $L$
(larger than D), being constrained to fall between two vertical walls
separated by a distance $L(1+\epsilon)$ with $\epsilon$ very small compared to
$1$. This setup is aimed to be designed to be close enough to a 2D cylinder
bouncing as characterised in the previous section, anticipating an effective
2D description of the experimental rebound. For this purpose, the motion of
the cylinder is imposed to take place in a plane perpendicular to its long
axis and along the vertical direction. In practice, three-dimensional motion
are unavoidable, due to both the finite length of the cylinder and wake
instabilities in the vertical plane. These three-dimensional motions require
to be controlled enough. For this purpose, the gap between the cylinder and
the lateral walls (L$\epsilon$) is always sufficient to allow the cylinder
falling but small enough to reduce three-dimensional rotational to a range of
small angles (typically $L\epsilon\simeq 0.1cm$). Moreover, to minimize three-
dimensional motion, the cylinder is provided with a long thin tail of
dimension $\ell\gg D$ and with a similar length $L$.

(a) (b) (c) (d)

Figure 1: (a) Design and picture of the four cylinders used. (b) Schematic of
the experiment, the delimited area in red is the field of view of the camera.
(c,d) Definitions of the extracted quantities, from real images.

In order to vary the Stokes number $St$ in the experiments, several cylinders
and fluids have been used. The different cylinders are shown in Fig. 1(a), all
of them having an outer core of PVC but their inner part can be replaced by
another material to modify the object density but not its surface/contact
properties. The tail is a plate made of PVC as well. The cylinder properties
are summarized in table 1. The rebound of the cylinders have been studied in
four different fluids at room temperature ($\sim 20\deg$): air, salt water,
and a mixture of salt water with Ucon oil (DOW©) 75H900000 at $10\%$ and
$15\%$ in mass. Their densities and kinematic viscosities are summarized in
table 2.

| A | B | C | D
---|---|---|---|---
inner material | PMMA | PVC | Steel | Steel
inner diameter (cm) | 1.9 | N.A. | 1.0 | 1.8
$D$ (cm) | 2.074 | 2.071 | 2.072 | 2.072
$L$ (cm) | 3.95 | 3.96 | 3.91 | 3.96
$\ell$ (cm) | 9.0 | 9.0 | 9.0 | 9.0
$\rho_{p}$ (g/cm3) | 1.36 | 1.43 | 2.17 | 3.19
Table 1: Properties of the cylinders shown in Figure 1(a). fluid |  $\rho$ (kg/m3) |  $\nu$ ($10^{-6}$ m2/s)
---|---|---
air | 1.17 | 15.7
salt water | 1050 | 1.1
ucon/salt water ($10/90$) | 1050 | 8.4
ucon/salt water ($15/85$) | 1050 | 17.7
Table 2: Properties of the fluids used. Mixture ratio are percentage in mass.

As shown in Fig. 1(b), the glass tank horizontal dimensions are $18.0\times
4.0$ cm2 and its vertical extent of $74.0$ cm allows for a falling distance of
approximately $70$ cm. The cylinder, almost fully immersed, is released at the
top of the tank by opening a gate (fast mechanical aperture using a compressed
air piston) and falls due to gravity in a fluid at rest until reaching the
bottom of the tank where it encounters a fixed flat surface made of PVC (same
material as the outer part of the cylinder). We follow its dynamics with a
high-speed camera connected to a telecentric lens to avoid parallax effects,
the sampling frequency chosen is 2kHz and the field of view is $1000\times
2000$ px2 with a resolution of $\sim 0.1$mm/px observing the $(x,y)$ plane as
shown in Fig. 1(b,c). The gray levels due to the back-lightning of the tank
allows to precisely extract the contour of the object, and even its
orientation with respect to the $(x,y)$ plane. A series of experiments have
also been investigated using a Phantom©high-speed camera, with the sampling
frequency at 130kHz (in salt water only).

For each grey-level image, we extract the contour and the center ($P_{1}$) of
the cylinder, as well as the contour and the barycenter ($P_{2}$) of the whole
object. From the temporal evolution of those two points, we can follow the
orientation in the $(x,y)$-plane of the tail $\phi(t)$ which corresponds to
the angle of the line $(P_{1},P_{2})$ with the vertical; the velocity for the
trajectories of $P_{1}$ and $P_{2}$ and their angle with the horizontal are
found by computing a linear fit of the trajectory over $\sim 20$ time-steps.
We define more specifically $\alpha_{in}$ and $\alpha_{out}$ the angles of the
trajectory with the horizontal, just before and after rebound, and
$\theta=|\alpha_{out}-\alpha_{in}|$ the difference between these two angles.
The coefficient of restitution in the vertical direction, denoted $e_{y}$ in
the following, will be computed by making the ratio of the velocity of $P_{1}$
just before and after contact. As illustrated in Fig.1(d), by using the gray-
level to discriminate the front (black) and rear (lighter gray) faces of the
cylinder, we can extract the difference in height $\delta z$ of their centers
which relates to the orientation out-of-plane of the cylinder
($\sin^{-1}(\delta z/L)$).

The release of the cylinder is always done in a fluid at rest for several
minutes, with the initial position of its tail as vertical as possible.
Nevertheless, residual perturbations in the fluid or at release can generate a
small rotation of the tail angle with the vertical ($\phi$), which can stay
almost constant in viscous fluids, but might generate lift and some rotational
motions in general. A precise monitoring of all these angles are done for each
experiment.

###  3.2 Results

####  3.2.1 Contact time and inclination at rebound

The rebound on the bottom wall observed with this experimental setup was shown
to be nearly 2D. However, a deeper investigation indicates a slightly more
complex bouncing feature, highlighting 3D contribution in the out-of-plane
direction. The complexity of the bouncing is associated with the out-of-plane
orientation of the cylinder with respect to the horizontal, as discussed in
the previous section and quantified by the length $\delta z$. This leads to a
contact point of the cylinder (either front or rear section touching the
bottom plate) instead of an idealized contact line along the entire cylinder,
as would be obtained for a purely 2D bouncing. The full bouncing process of
such apparent 2D bouncing can actually be characterized by a succession of
bouncing of the front and rear sections of the cylinder on a very short time
scale, prior to a significant take off of the object indicating the end of the
apparent 2D bouncing. Such complex bouncing can be classified by the number of
visible contacts between the cylinder and the plate. It can have one (front or
rear part depending on its out-of-plane inclination), two (front and rear) or
even three (front/rear/front or vice-versa) contacts with the PVC plate before
going away from it. Accordingly, we can also define an apparent contact time
$t_{C}$ that corresponds to the delay between the first and last images
showing contact, with a lower bound being the frame rate of the camera
($0.5$ms) for the $1$-contact cases. Two examples are illustrated in Fig.
2(a,b), corresponding movies are given in supplementary material 111see
WATER_CYLA_1CONTACT_2000FPS.gif, WATER_CYLD_3CONTACTS_2000FPS.gif, and
WATER_CYLD_1CONTACT_CAVITATION_130000FPS.gif.

\begin{overpic}[trim=0.0pt 0.0pt 0.0pt
0.0pt,clip={true},width=411.93767pt]{all_contacts_v7.png} \put(-1.0,24.0){(a)}
\put(-1.0,10.0){(b)} \put(65.0,24.0){(c)} \end{overpic} Figure 2: Series of
images to illustrate rebounds with (a) 1-contact (front contact) for cylinder
A in salt water, and (b) 3-contact (rear/front/rear front contact) case for
cylinder B in salt water. The frames with contacts are indicated in red, a
dashed-dot white circle shows the front face of the cylinder. (c) Sorting of
rebound cases based on $\delta z/R$, the non-dimensional contact time and
number of contacts for cases in uconn/salt water (blue) and salt water (red).
Cylinders are sorted by symbols.

We have also registered the out-of-plane inclination at the rebound, or
equivalently $\delta z/R$, to relate it to the nature of the rebound. On
average, image analysis led to estimates of $\delta z/R$ smaller than $0.1$,
which corresponds to inclinations smaller than $5\deg$. Even if the non-
dimensional apparent contact time $t_{C}V_{t}/D$ roughly increases with
$\delta z/R$, dispersion of the results does not allow to assure that they are
strongly correlated, as illustrated in Fig. 2(c). One can notice however that
it is more likely to have two or three contacts at rebound when the out-of-
plane inclination is important, i.e. for increasing $\delta z/R$.

####  3.2.2 Influence of cavitation at rebound

For some of our experiments, the estimated values for $t_{C}$ is at the limit
of the sampling frequency of the camera (2kHz). In order to better resolve
those values, a series of experiments have been repeated with a camera working
at very high-speed (130kHz). As can be seen in Fig. 3(a) for a cylinder (type
D) in water, the overall vertical dynamics is well represented by two branches
of trajectories at constant speed. although some out-of-plane inclination is
visible after the rebound. If we focus on the instants around the rebound
defined as the origin of time, illustrated in the inset in (a), we notice that
the cylinder is actually staying almost steady very near the bottom for almost
$600\mu s$ after the contact, before starting to move away from it.

\begin{overpic}[trim=0.0pt 0.0pt 0.0pt
0.0pt,clip={true},width=411.93767pt]{cavitation_v1.png} \put(2.0,36.0){(a)}
\put(47.0,36.0){(b)} \end{overpic} Figure 3: Cylinder D in water. (a) Position
of the front and rear faces of the cylinder (as well as the average of the
two) with time, with the inset zooming on the $1ms$ around the rebound. (b)
Series of images from the camera (every $200\mu s$) with symbols indicating
the corresponding times in the inset in (a), the red arrow indicates the
cavitation bubble when visible.

Some images every $200\mu s$ of the rebound are shown in Fig. 3(b), with
symbols indicating the corresponding times in (a), a short movie of this
dynamics is also provided in supplementary material as well (see supplementary
material). One can observe the formation of a cavitation bubble at the
location of contact, in between the cylinder and the bottom, as indicated by
red arrows in 3 of the images. Its maximum horizontal extent is almost half
the diameter of the cylinder. The nucleation of this bubble originates from
very small bubbles at the surface of the cylinder and of the flat surface
(visible in the image with the disc, for $t=0\mu s$). This is in good
agreement with the crevice model for heterogeneous nucleation of bubbles in
water [24], with nuclei trapped at the plastic surfaces being destabilized at
rebound by pressure variations.

Similar observations have been obtained for all the rebounds with a cylinder
of type C or D, in salt water and in uconn/water mixtures as well. The
cavitation process seems to be related to the pressure drop generated at the
rebound for sufficiently dense cylinders. The velocity at rebound is not the
key element since no cavitation has been observed with cylinder of type B in
water (falling faster than C in ucon/salt water). Some tests on the influence
of degassing the water of the tank before the experiments have been done
without noticing a difference in the process. Further investigations are
needed, although they are out-of-scope for this manuscript.

To conclude, the cavitation bubble holds for about $200\mu s$ to $600\mu s$
(depending on the cases studied) and it is reproducible for similar
conditions. It is found here that cavitation is a process that is also
controlling the apparent time of contact $t_{C}$ for some cylinders (C and D),
preventing it to be shorter. Altogether, it is thus found that the concept of
contact time can be associated with several complex mechanisms.

####  3.2.3 Non-normal rebound

(a) (b)

Figure 4: (a) Distribution of the observed angles in air and in salt water at
rebound of the cylinder itself $\phi_{in}$, the initial orientation
$\alpha_{in}$, and of the change in the orientation of its trajectory with
respect to a perfect rebound with no friction $\alpha_{in}-\theta/2$. (b)
coefficient of restitution $e_{y}$ in air, as a function of the rotation at
rebound, sorted by the nature of the cylinder (symbols); the filled red
symbols with error bars are the estimates of $e_{y_{0}}$ since $St>10^{5}$ for
all cases.

The other parameters that can influence the rebound are related to the
orientation of the cylinder and its trajectory. For all experiments in air and
salt water, we represent the distribution of the orientation at contact
$\phi_{in}$, the orientation ($\alpha_{in}$) of its trajectory, and the change
in orientation with respect to a perfect rebound with no friction
$\alpha_{in}-\theta/2$ in Fig. 4(a). One can notice that the orientation at
contact $\phi_{in}$ is peaked near $[0,1]$ degree indicating that most cases
are nearly vertical, although some non-vertical orientation can occur. The
evolution of this orientation after rebound is almost unnoticeable and there
is no correlation of $\phi_{in}$ with the coefficient of restitution (not
shown). Related to this nearly vertical orientation, the distribution of angle
the trajectory of the cylinder before rebound is also very peaked near $[0,1]$
degree although perturbations in the fluid tank or at release can induce some
fluctuations that are not correlated with a lift force associated to
$\phi_{in}$ being non-zero. Finally, the most important quantity to measure
the non-normal aspect of the rebound is the distribution of the change in the
orientation of the trajectory of the cylinder, $|\alpha_{in}-\theta/2|$. The
pdf shown is nearly flat for values between $0$ and $10$ degrees, with few
cases larger than $10$ degrees. This quantity measures the importance of
friction in the $x$-direction that can induce not trivial trajectories at
rebound as discussed in [25].

####  3.2.4 Rebound in air

In Fig. 4(b), we focus on the experiments done in air. They correspond to very
large values of the Stokes number and are used to define a reference
coefficient of restitution $e_{y_{0}}$ solely controlled by elasticity (red
symbols). The estimates for $e_{y_{0}}$ are indicated with errorbars to
account for the influence of the rotation of the trajectory at rebound. The
values are slightly different for each object (decreasing from A to D),
indicating that the inner core might play a part on the elastic response of
each object. More specifically, the value of $e_{y}^{0}$ decreases with
increasing mass of the object, which can also be explained by the dissipation
increasing in the fixed plate at the bottom.

####  3.2.5 Synthesis

All the specificity of the apparent rebound of a cylinder discussed before can
influence the coefficient of restitution in the vertical direction, and
generate an important dispersion of the experimental observations.
Nevertheless, we have not clearly identified a correlation between $e_{y}$ and
the contact time, out-of-plane orientation or number of contacts at rebound,
except for 3-contacts that always correspond to weaker values of $e_{y}$. In
the following, to compare with idealized numerical simulations, we will
consider only observations for which $\alpha_{in}-\theta/2$ strictly smaller
than $6$ degrees, and with 1 or 2 contacts eventually. We could label these
cases ’nearly-normal’ rebounds of a cylinder.

The coefficient of restitution $e_{y}$ obtained in these selected experiments,
according to the above discussion, is shown as a function of $St$ in figure 5,
the shape and color of symbols indicating the cylinder and fluid properties as
described in tables 1 and 2. Error-bars in $e_{y}^{0}$ are discussed in the
previous section, but we can also notice some dispersion of the results in
salt water and in ucon/salt water. for which each symbol is a given
experiment, a signature of the many (uncontrollable) parameters at play
(cylinder inclination, pitching, cavitation,…). The main result here is the
estimation of a threshold $St_{c}$ above which rebound actually occurs. Above
$St_{c}=75\pm 25$, $e_{y}$ increases with $St$ until reaching $e_{y}^{0}$ for
high enough values. Furthermore, it can be noticed that the ”S-shape” that has
been already evidenced in the past studies on spheres bouncing on a wall is
still observed in the experiments on the falling cylinder. This suggests a
possible applicability of the concept of previous established models, provided
that those models are extended to take into account the specificity of the
contact in the actual problem, compared to the contact occurring between a
sphere and a plane. We will show in the following section, that it is
sufficient to use the apparent (measured) contact time in order to predict the
coefficient of restitution as a function of the Stokes number.

Figure 5: The coefficient of restitution as a function of the Stokes number
for experiments corresponding to nearly normal rebound. Colors of symbols
indicate the fluid considered, while the shape refers to the properties of the
cylinders. The critical value for rebound, $St_{c}=75\pm 25$, is indicated by
the black dashed line with gray shaded area.

##  4 2D effective bouncing

###  4.1 Numerical simulations of the rebound of a 2D cylinder

$0$$L_{y}$$L_{x}$(a)(b)(c) Figure 6: (a) 2D numerical simulation setup,
including the grid distribution and boundary conditions. Evolution in time of
(b) the separation gap between the cylinder surface and the wall and (c) the
cylinder velocity. The inset in (b) shows a zoom at small gaps. The colors
correspond to different Stokes numbers indicated in (c). The black dots
indicate the collision velocity $V_{c}$. The elasticity parameter used for
those simulations is $\alpha=4$.

The numerical setup is sketched in figure 6(a). It consists of a two-
dimensional closed vessel of dimension $(Lx,Ly)=(20R,100R)$ along the vertical
$y$ direction and the horizontal $x$ direction, respectively, filled with a
viscous fluid. No-slip velocity is imposed on three boundaries considered as
solid walls. The fourth boundary $x=0$ (along the $y$ direction) is a line of
symmetry (see figure 6(a) for details). A 2D cylinder of radius $R$ is freely
settling under an imposed force $g$ along the negative $y$ direction. The
cylinder center $(x_{p},y_{p})$initially placed at the line of symmetry is
constrained to move along $y$ and thus remains at $x=0$. The position of the
centre of the 2D cylinder is therefore denoted $(0,y_{p})$, with $y_{p}$
depending on $t$. Note that using advantage of the properties of symmetry of
this configuration, only half of the cylinder is represented and simulated
here.

The coupling between fluid and 2D cylinder equations of motion is based on the
Immersed Boundary Method, that we previously used to study the rebound on a
wall of i) a sphere settling under a constant force [20] and ii) a sphere
carried by a wall-normal flow toward the stagnation point [21, 26]. However in
the present work, unlike our work on the settling sphere [20], we do not
include any lubrication correction when the cylinder is close to the wall, as
the motion of the fluid squeezed between the cylinder surface and the wall is
resolved down to a scale relatively small compared with the cylinder radius
$R$ as explained in [21, 26]. In particular, near the contact region,
$(x,y)\approx(0,0)$, a fine grid resolution is used: almost 20 grid points are
employed in the region $y<0.01R$. This limits the minimal apparent roughness
scale of the particle to a small fraction of the particle radius, as will be
explained later on. Elsewhere, non-uniform Cartesian grid is used for
computational efficiency. The size of the grid elements increases smoothly
along positive $x$ and $y$ directions. At large distance from the wall, the
spatial discretization ensures 30 grid points per cylinder diameter in the $y$
direction, and slightly more in the $x$ direction.

Simulations are carried with constant cylinder diameter, and constant fluid
and cylinder densities. Cylinder inertia is varied by changing the fluid
viscosity and thus the Reynolds number while the density ratio is kept
constant. The Reynolds and Stokes numbers defined from 1 and 2 are varied in
the range $[1-500]$. Note that using a symmetrical domain leads to a symmetric
wake behind the cylinder. Even if wake instability shall be expected at the
highest $Re$, this assumption remains acceptable as the settling time scale of
interest for the bouncing is much smaller than the time scale required for the
wake instability to take place, and therefore to affect the results presented
here. During the settling stage, and after rebound, the time step is set to a
small fraction of the settling time $dt\approx 0.001d/|V_{t}|$. During the
collision stage, a smaller time step, equal to $1/40$ of the collision time
estimated from eq. (4) as explained later.

The 2D cylinder settles from rest, its center being initially located at
$y=90R$. It accelerates until the velocity reaches a terminal velocity
$V_{t}$, balancing drag and apparent weight, prior to sudden deceleration due
to a large hydrodynamic resistance close to the bottom wall (see the temporal
evolution of the cylinder velocity in figure 6(b)). From the simulations, the
cylinder terminal velocity $V_{t}$ is obtained from the temporal signal of the
cylinder velocity (as shown in figure 6(b)). The plateau corresponding to the
balance between drag and apparent weight is clearly observed for $St<100$,
between the acceleration and deceleration stages. At larger $St$, the cylinder
falling time is not sufficient to reach the plateau. In this case, the
terminal velocity is then defined as the maximum (negative) velocity before
the cylinder starts to decelerate.

At the end of the deceleration stage the cylinder comes to rest at low
inertia, i.e. small $St$. However at high inertia, larger $St$, the cylinder
surface becomes critically close to the wall while the cylinder settling
velocity remains finite. When the gap between the cylinder surface and the
wall becomes smaller than a threshold $\eta R$, referred to as the apparent
roughness length ($\eta$ is the contact roughness), solid contact is assumed
to occur. The finite velocity of the particle at the contact onset called
$V_{c}$ hereafter, will be taken as the particle velocity when
$\delta_{n}=y_{p}-(1+\eta)R$ becomes negative. While the value of $\eta$ is
set to few percents in the numerical simulations, we carefully verified that
the fluid motion in the separation gap is fully resolved with more than 20
grid points during solid contact, in a way to capture correctly the
lubrication effect without using sub-grid models. Solid collision is then
accounted for in the numerical simulation by adding a contact force to the
cylinder equation of motion (in the wall-normal direction), as soon as
$\delta_{n}$ becomes negative. The contact force is modeled as a linear
elastic force $F_{c,n}=-k_{n}\delta_{n}$ with spring stiffness $k_{n}$. The
solid dissipation during the contact is neglected in order to maintain viscous
lubrication as a unique source of energy dissipation. We use $\delta_{n}$ to
represent the solid apparent overlapping during the contact which mimics
elastic deformation. The spring stiffness per unit length is thus modeled as

$k_{n}=\frac{m\ \pi^{2}}{L\left(\alpha T_{H}\right)^{2}},$ (3)

where $m/L$ denotes the cylinder mass per unit length. $\alpha T_{H}$ sets the
characteristic time scale of collision (during which the spring force is
turned on). Inspired by the Hertzian collision time scale used to interpret
the experimental results for spherical particles colliding with a wall [27]
(based on reference [28]), $T_{H}$ adapted to cylindrical particle is written
as:

$T_{H}=\left[\left(\frac{m}{L}\
\frac{1}{E_{el}}\right)^{2}\frac{R}{2V_{c}}\right]^{0.2},$ (4)

where $E_{el}$ denotes the Young Modulus of the cylinder, assumed equal to
that of the wall. If the collision between the cylinder and the wall was
taking place in dry conditions (in the absence of viscous liquid), the contact
duration would have been equal to $\alpha T_{H}$. However in viscous liquid,
the contact duration is larger than this value , the difference being a
decaying function with the Stokes number. The parameter $\alpha$ allows to
tune the contact elasticity in a practical way in the simulations. Thus the
effective Young modulus is equal to $\frac{E_{el}}{\alpha^{2.5}}$. Softer
(resp. stiffer) collisions occur for larger (resp. smaller) $\alpha$. A
dimensionless number, $S=\frac{\mu V_{c}/(\eta R)}{E_{el}/\alpha^{2.5}}$, is
built in a way to compare the shear stress of the flow in the separation gap
during the collision process (estimated by $\mu V_{c}/(\eta R)$ and varies
with the Stokes number) and the solid elasticity. The meaning of this
dimensionless number is similar to that defined in the elastohydrodynamic
theory of [8]. Figure 7 shows the contour plot of $S$ corresponding to the
conditions where numerical simulations have been carried. The range of $S$
obtained here is close by construction to the Capillary number corresponding
to slightly deformable drops. Figure 7 suggests that $S$ depends weakly on the
Stokes number at small $\alpha$, which means that the contact softness is
mainly tuned by the elasticity parameter $\alpha$.

Figure 7: Contour plot (in natural log-scale) of the softness dimensionless
number $S=\frac{\mu V_{c}/(\eta R)}{E_{el}/\alpha^{2.5}}$ for all numerical
simulations used for this study. The value of this parameter ranges between
$10^{-3}$ (the darkest color) and $O(10)$ (the lightest color). The dots
indicate the parameters used in the different simulations carried with
$\eta=0.01$.

From the simulations we extract the terminal velocity $V_{t}$, the velocity
$V_{c}$ at the onset of collision process (at the instant where $\delta_{n}$
becomes negative) and the velocity $V_{r}$ at the end of the collision process
(at the instant when $\delta_{n}$ becomes positive again). In the simulations,
the ”dry” coefficient of restitution is equal to one (the contact force
$F_{c,n}$ does not contain any damping term). The effective coefficient of
restitution is defined as the ratio $e=-V_{r}/V_{t}$ which can be written as
the product of two contributions, prior to solid contact $V_{c}/V_{t}$ and
during contact $-V_{r}/V_{c}$. These two contributions are plotted as a
function of $St$ in figure 8. Panel a) in this figure shows that $V_{c}/V_{t}$
increases with $St$, which suggests that, at the onset of rebound, the
cylinder kinetic energy increases with $St$. Panel b) shows that
$-V_{r}/V_{c}$ is close to 1 and almost independent of $St$ when the
elasticity parameter is small. The magnitude of $-V_{r}/V_{c}$ is
significantly decreased when the elasticity parameter is increased (softer
contact).

To summarize, the contact in the numerical system is controlled by two
parameters linked to two different physical processes: the elasticity
parameter $\alpha$ which mainly controls the effective contact time of the
collision and therefore the contact softness/rigidity, and the cut-off length
$\eta$ that sets the collision onset, akin an apparent roughness. It will be
shown in the following that both play a role on the coefficient of
restitution.

###  4.2 Model for the restitution coefficient

This section explains the basis of a simple model constructed to allow the
prediction of both the cylinder velocity at the onset of collision $V_{c}$,
and the rebound velocity $V_{r}$ when the solid contact is completed. This
model follows closely the one outlined in the work of Izard et al. [20]
(considering the rebound of a sphere on a wall).

(a)(b) Figure 8:  (a) Collision-to-terminal velocity ratio $V_{c}/V_{t}$ and
(b) rebound-to-collision velocity ratio $-V_{r}/V_{c}$ as a function of the
Stokes number. $V_{c}$ is the particle velocity at the collision onset, for
$\eta=0.01$. The symbols correspond to numerical simulations. The solid line
in panel a) predicts $\frac{V_{c}}{V_{t}}$ from eq. 7. In panel (b) different
symbols correspond to different $\alpha$: stars, plus, up triangles, down
triangles and squares correspond to $\alpha=0.5$, 1, 4, 7 and 10, respectively
(from rigid to soft cylinder). The solid lines in the inset correspond to
$-\frac{V_{r}}{V_{c}}$ obtained from eq. 9 using a collision time
$\tau=\pi\sqrt{\frac{m^{*}}{k_{n}}}$ and $\alpha=$ 1, 4, 7 and 10 (from top to
bottom).

First, let us attempt to estimate the ratio between the collision velocity
$V_{c}$ and terminal velocity $V_{t}$, by writing the equation of motion of
the cylinder during the deceleration stage. While unsteady forces that the
cylinder might experience during the deceleration stage are not readily
available in the regime of interest (Stokes number ranging between 10 and
100), we write a simplified equation of motion of the cylinder, in a way to
predict the motion at the leading order. For this, we assume that the
hydrodynamic force experienced by the particle is the superposition of the
drag experienced by a cylinder moving in an unbounded fluid, the viscous
lubrication associated with the pressure divergence as the gap between the
cylinder surface and the wall $\zeta=(y_{p}-R)\rightarrow 0$ and an added mass
contribution

$m^{*}\frac{dV}{dt}=F_{lub}.$ (5)

In eq. (5), $V$ denotes the cylinder velocity and $m^{*}=\rho^{*}\pi R^{2}$
corresponds to the cylinder mass per unit length accounting for the added mass
effect. $\rho^{*}=\rho_{p}\left(1+C_{M}\rho_{f}/\rho_{p}\right)$ corresponds
to the cylinder apparent density accounting for the added mass effect ($C_{M}$
refers to the added mass coefficient). Here, the gravity force and steady drag
corresponding to a cylindrical body moving in an unbounded fluid are omitted,
assuming that their balance is independent of the cylinder position. In other
words, eq. (5) describes the departure from the equilibrium (where the
settling velocity equals the terminal velocity $V_{t}$) due to the presence of
the wall. The lubrication force (per unit length) experienced by a 2D cylinder
moving toward a wall with a given velocity and assuming the flow motion is
quasi-steady follows (see appendix 7.1):

$F_{lub}=-\frac{6}{\sqrt{2}}\pi\mu\left(\frac{R}{\zeta}\right)^{3/2}\frac{d\zeta}{dt},$
(6)

where $d\zeta/dt$ corresponds to the particle instantaneous velocity in a
Eulerian frame of reference, i.e. $d\zeta/dt=V$. Although this expression is
strictly valid at $\zeta<<R$, we assume that it applies continuously from a
distance $\zeta=O(10R)$ (where $F_{lub}$ tends to zero) until very small
$\zeta$. Of course in the region $\zeta=O(R)$ where the particle decelerates,
this approximation is not valid. An accurate force balance on a cylindrical
particle approaching a wall, where unsteady contributions are accounted for,
would be necessary, however not available currently.

Integration of (5)-(6) from a distance where the cylinder starts to decelerate
($V\approx V_{t}$ at $\zeta\approx R$) to the wall where collision occurs
($V=V_{c}$ at $\zeta=\eta R$), assuming $\eta\ll 1$ leads to (see appendix
7.2):

$\frac{V_{c}}{V_{t}}=\left[1-\frac{St_{c}}{St}\right],\quad
St_{c}=\frac{8}{3\sqrt{2}}\eta^{-1/2},$ (7)

$St_{c}$ being a critical Stokes number for bouncing, i.e. over which
$V_{c}\neq 0$. Then, when $St<St_{c}$, the cylinder energy is entirely damped
away before it reaches the wall, i.e. $V_{c}/V_{t}\rightarrow 0$ when
$\zeta\rightarrow\eta R$. However above the rebound onset, i.e. $St>St_{c}$,
the ratio $V_{c}/V_{t}$ increases with particle inertia. Its value extracted
from the numerical simulations is compared with the solution of eq. (7) for
varying $St>St_{c}$ in figure 8(a) (symbols correspond to numerical simulation
and solution of (7) is plotted as solid line). Despite its obvious simplicity
to obtain a theoretical solution, the estimation of $\frac{V_{c}}{V_{t}}$ by
the model (7) is acceptable.

(a)(b) Figure 9: Coefficient of restitution as a function of the Stokes
number; the symbols are from numerical simulations and the dashed lines
correspond to $e_{RB}$ (eq. 10) assuming that the cylinder is infinitely rigid
and that the near-wall hydrodynamic force during settling is dominated by
lubrication. Panel (a) shows evidence of the role of apparent surface
roughness: the blue plus and empty red circles correspond to $\eta=0.01$ and
$0.1$, respectively, obtained with the elasticity parameter $\alpha=1$. Panel
(b) evidences the role of elasticity: the stars, plus, up triangles, down
triangles and squares correspond to simulations carried with $\alpha=0.5$, 1,
4, 7 and 10, respectively and apparent roughness $\eta=0.01$. The dashed line
is reported from panel (a). The solid line corresponds to $e_{Iz}$ (eq. 12)
assuming a balance between elastic deformation and lubrication forces during
contact. The dotted lines are calculated from
$e=-\frac{V_{c}}{V_{t}}\times\frac{V_{r}}{V_{c}}$, where $\frac{V_{r}}{V_{c}}$
is obtained from eq. 9 with $\tau=\pi\sqrt{\frac{m^{*}}{k_{n}}}$ like in
figure 8b. These dotted lines correspond to $\alpha=$ 1, 4, 7 and 10, from top
to bottom (from rigid to soft contact).

Above the rebound onset, the cylinder velocity at the end of the collision
$V_{r}$ is obtained from the equation of motion of the cylinder, accounting
for viscous lubrication and contact force. For the latter, we consider the
elastic force $F_{c,n}$ identical to the one used in the numerical
simulations. Again, dissipation in the solid is neglected. This leads to the
damped oscillator equation, written in terms of the overlapping distance with
respect to the contact cut-off length scale $\delta_{n}=\zeta-\eta R$:

$m^{*}\frac{d^{2}\delta_{n}}{dt^{2}}+\lambda\frac{d\delta_{n}}{dt}+k_{n}\delta_{n}=0\quad\mbox{with}\quad\lambda=\frac{6}{\sqrt{2}}\pi\mu\left(\frac{1}{\eta}\right)^{3/2}.$
(8)

The rebound velocity $V_{r}$ can then be estimated as $V_{r}=d\delta_{n}/dt$
at $t=\tau$, where $\tau$ denotes the time at which the overlapping distance
is again equal to the cut-off length $\eta R$. Based on this definition,
$\tau$ represents the solid contact time. The obtained rebound-to-collision
velocity ratio reads

$\frac{V_{r}}{V_{c}}=-exp\left[-\frac{\lambda\tau}{2m^{*}}\right].$ (9)

According to this model, the ratio $V_{r}/V_{c}$ depends on the particle
stiffness via the solid contact time $\tau$. The longer is the contact time,
the smaller is the magnitude of the rebound velocity. Except very close to the
rebound threshold, where viscous effects are important, the contact time is
controlled by contact elasticity and it is reasonable to approximate it by
$\tau=\pi\sqrt{m^{*}/k_{n}}$. Figure 8(b) shows the ratio $V_{r}/V_{c}$
obtained from both numerical simulations and eq. (9), for different elasticity
parameter $\alpha$, ranging from $0.5$ to $10$. Overall, the model captures
the trend found in the simulations, i.e. the ratio $V_{r}/V_{c}$ increases
with $St$ and decreases with the elasticity parameter. The overestimation of
$V_{r}/V_{c}$ by the model is likely due to inertial effects that are not
accounted for in eq. (8), as for instance the contribution of the wake behind
the cylinder to the dynamics during motion reversal.

Finally, the effective coefficient of restitution is defined using the product
of both velocity ratios, i.e.
$e=-\frac{V_{c}}{V_{t}}\times\frac{V_{r}}{V_{c}}$. Figure 9 shows the separate
effect of the cut-off length $\eta R$ (panel a) and elasticity parameter
$\alpha$ (panel b) on $e$, while increasing the Stokes number. Figure 9(a)
obtained with relatively rigid contact ($\alpha=1$) highlights that the cut-
off $\eta$ mostly affects the critical Stokes $St_{c}$, i.e. the onset of
solid contact. The larger the apparent roughness, the earlier the no-rebound
to rebound transition takes place. In order to support this observation,
numerical results are compared with the coefficient of restitution of an
infinitely rigid body bouncing on a wall, in which case
$V_{r}/V_{c}\rightarrow-1$ (no loss during collision), corresponding to an
infinitely small solid contact time, i.e. $\tau\approx 0$. In this limit, the
coefficient of restitution called becomes equal to:

$e_{RB}=\left(1-\frac{St_{c}}{St}\right).$ (10)

$e_{RB}$ is displayed with dashed lines in figure 9(a), with $\eta=0.01$
(blue) and $\eta=0.1$ (red). One clearly observes that $e(St)$ obtained from
numerical simulations follows closely the rigid body model, as long as the
appropriate $St_{c}$, set by the apparent contact roughness modeled by $\eta$,
is used.

Figure 9(b) shows that the cylinder elasticity (varied by the parameter
$\alpha$) impacts significantly the evolution of the coefficient of
restitution at a given Stokes number. For the most rigid cylinders,
$\alpha=0.5$ and $1$, the curves of the coefficient of restitution almost
match together and approach the rigid body limit ($e_{RB}$). An increase of
$\alpha$, i.e. softer contact, leads to longer contact duration and thus to a
more significant contribution of the viscous dissipation to the overall energy
loss during the collision process ($\frac{V_{r}}{V_{c}}$ decreases).
Consequently the slope of the $e-St$ curve is lowered, and the transition of
the coefficient of restitution from $0$ to $1$ is expected to take place on a
much wider range of Stokes numbers. Simulation results are compared to the
model $e=-V_{r}/V_{t}$ (the product of (7) and (9)) displayed in dotted lines.
The relatively good agreement, suggests that the proposed simple model
contains the required ingredients to capture the energy loss during the 2D
cylinder bouncing.

We finish this section by reminding the reader that the elasticity parameter
used here in both the model and numerical simulations remains a model of
elastic effects during bouncing. The analogy with experiments is not trivial
(as will be discussed in the next section). In order to anticipate a
parameterization of an apparent elastic process for complex bouncing, which
could be induced by more complex hydro-mechanical process, but still affecting
the time scale of bouncing, we test here a simple scaling of the solid contact
time $\tau$ following [20]. For this purpose, if one assumes a balance between
the force of elastic deformation that scales like $k_{n}\eta$ and the
lubrication force that scales like $\lambda|V_{c}|$, the ratio $V_{r}/V_{c}$
can be written in a way that does not depend on the cylinder rigidity, but
only on particle inertia and roughness which sets the critical Stokes number
$St_{c}$

$\frac{V_{r}}{V_{c}}=-exp\left[-\frac{\pi}{2\sqrt{2}}\sqrt{\frac{St_{c}}{St-
St_{c}}}\right].$ (11)

This model can be coupled with (7) to obtain $e_{Iz}$, leading to another
interesting limit for the coefficient of restitution (and for the solid
contact time, shown in appendix 7.3) which is compared against numerical
simulations and experiments (figures 9 and 10),

$e_{Iz}=\left(1-\frac{St_{c}}{St}\right)exp\left[-\frac{\pi}{2\sqrt{2}}\sqrt{\frac{St_{c}}{St-
St_{c}}}\right]\,.$ (12)

It is interesting to note at this stage that both limits $e_{RB}$ and $e_{Iz}$
allow to delineate the region of transition of the ensemble of numerical
simulations, i.e. for a large range of contact stiffness.

##  5 Discussion and conclusion

The study carried in the previous section with a 2D cylinder suggests that (i)
the increase of dissipation associated with the contact elasticity is well
captured by a model based on solid contact time scale $\tau$, leading to
relation (9) for the velocity ratio $V_{r}/V_{c}$, and that (ii) the cut-off
length $\eta$ impacts mainly the velocity ratio $V_{c}/V_{t}$ following
relation (7). We will examine here to which extent these concepts apply to our
experimental measurements. The panels in figure 10 display the variation of
the coefficient of restitution (a) and the contact time (b) as a function of
the Stokes number, from experiments and numerical simulations. The Stokes
number is scaled by the critical Stokes number corresponding to the onset of
bouncing in different cases. Regarding the results from experiments, we remind
the reader that only those obtained from observations when the rebound is
nearly-normal, and when the contact dynamics is with 1 or 2 contacts, are kept
in this figure 10. This induces some dispersion of the results, due to the
various processes at contact between the falling object and the wall, as
classified in section 3.2. This is different from measurement uncertainties,
with vertical error-bars associated with uncertainties in $e_{y}^{0}$ and
horizontal ones with uncertainties in $St_{c}$, when changing its value from
$50$ to $100$.

(a) (b)

Figure 10: (a) The coefficient of restitution normalized by its value in air
as a function of the Stokes number and (b) the non-dimensional contact time as
a function of the Stokes number scaled by $St_{c}$. Black symbols are similar
to figure 9(b), obtained with different elasticity parameters $\alpha$. Red
symbols correspond to experiments as in figure 5. Vertical errorbars for
experiments in (a) are associated to uncertainties in $e_{y}^{0}$ and
horizontal ones to uncertainties in $St_{c}$ (not reproduced in panel b).
Empty (resp. filled) symbols are for 1-contact (resp. 2-contact) rebounds. The
black dashed line in (a) is for $e_{RB}$. In (a,b), the blue solid line is for
$e_{Iz}$ and $\tau_{Iz}$ with $\eta\in[3.6-15]\,10^{-4}$, while the green and
cyan lines correspond to models $e_{m1}$ and $e_{m2}$ with $St_{c}=75$
($\gamma=2$) and apparent non-dimensional contact time being constant at
$(6\pm 3)\,10^{-2}$ or fitted by $(8\pm 4)\,10^{-3}\sqrt{St/St_{c}-1}$
respectively. Shaded areas around the models describe the uncertainties on the
fitting parameters or the value of $St_{c}$.

Let’s examine in figure 10(a) the coefficients of restitution from experiments
and simulations. Despite the dispersion, we can derive a general trend based
on previous modelling. The departure of $e_{y}/e_{y}^{0}$ from zero near
$St_{c}$ is assumed to be controlled by the apparent contact roughness,
associated in the experimental observations with the surface properties of the
material considered (PVC in our case). If contact elasticity is first ignored,
the coefficient $e_{y}/e_{y}^{0}$ could be estimated from the rigid body
approximation, i.e. following eq. (7), shown by the black dashed line. This
model is only characterized by the threshold in Stokes number; the
experimental estimate of $St_{c}=75$ corresponds to a cut-off length of
$\eta=6.3\times 10^{-4}$ according to eq. (7). It is worth noting here that
this expression for $St_{c}$ (obtained with the definition of the Stokes
number in eq. (2) inspired from the falling sphere) should be adjusted to
account for the shape of the experimental object, i.e. a cylinder with tail.
However, the form of the relationship between $\eta$ and $St_{c}$ still holds
as

$St_{c}=\gamma\eta^{-1/2},$ (13)

where $\gamma$ shall be a shape dependent constant of order unity
($\gamma\simeq 1.88$ for the cylinder without tail). Whatever the value of
this constant, the value of $\eta$ does not seem to reflect the range of
$\delta z/R$ reported in figure 2(c). Varying $St_{c}$ from $50$ to $100$
(with $\gamma=2$) would lead to $\eta\in[4.0\times 10^{-4}-1.6\times
10^{-3}]$, that is of the order of the roughness measured for instance for
Nylon particle surface [4], but still below the range of $\delta z/R$ observed
in the experiments. Only values of $\gamma$ larger than $10$ could allow for a
better match. Hence the apparent roughness is more likely to relate to the
material roughness, which is identical for all the cylinders (outer shell in
PVC bouncing on a PVC plate).

Nevertheless, the rigid body approximation is clearly an upper bound for the
evolution of $e_{y}/e_{y}^{0}$ as a function of $St/St_{c}$. This observation
suggests that the elasticity parameter $\alpha$ in the model, and its apparent
counterpart in experiments, has to be considered. We recall here that by
apparent contact elasticity, we refer to the fact that the contact time
occurring during this overlap scale is finite.

Figure 10(b) displays jointly the contact time measured from numerical
simulations (the time where the separation gap is smaller than $\eta R$) and
the apparent contact time $t_{C}$ from the experiments, as a function of
$St/St_{c}$. The contact time is scaled by the settling time scale $D/V_{t}$.
This figure suggests that i) the dimensionless contact time obtained from the
simulations (black symbols) varies significantly with the elasticity parameter
$\alpha$ but more weakly with the Stokes number, and that ii) it falls in the
same range of contact time measured experimentally. The latter can be grouped
in two families. One group can be described by $t_{C}V_{t}/D=(6\pm 3)\times
10^{-2}$, corresponding mainly to the collisions with 2-contacts which lead to
a large apparent contact at high Stokes numbers (filled red symbols in figure
10b). As for the other group (empty red symbols in 10b), the contact time
seems to be adequately described by a law of the form $t_{C}V_{t}/D=(8\pm
4)\times 10^{-3}\sqrt{St/St_{c}-1}$. Note that the empty symbols correspond to
1-contact rebounds and cavitation has been observed for the cylinders of type
C (squares) and D (triangles), even in the low $St/St_{c}$ range. In the first
group where the dimensionless contact time can be considered as constant (and
relatively large), the contact time is close to, but exceeds the largest
contact time measured in the simulations at large $\alpha$ (soft contact).
This observation agrees with the experimental coefficient of restitution being
globally smaller than the numerical value at a given $St/St_{c}$. However in
the second group, the contact time increases by more than one order of
magnitude in the considered range of $St/St_{c}$, in a way consistent with the
influence of the elasticity parameter $\alpha$ on the contact time obtained
from numerical simulations (black symbols).

We shall discuss in the following the possibility of modelling the coefficient
of restitution given a law for the contact time. From eqs. (7) and (9), if one
uses the relation between cut-off length $\eta$, the critical Stokes number
following eq. (13) and the definition of $St$ in (2), one obtains

$e=\left(1-\frac{St_{c}}{St}\right)\exp\left[-\frac{\pi}{3\sqrt{2}}\left(\frac{St_{c}}{\gamma}\right)^{3}\frac{1}{\tilde{S}\,St}\frac{t_{c}V_{t}}{D}\right]$
(14)

where $\tilde{S}=S_{obj}/D^{2}$ with $S_{obj}$ the 2D surface of the falling
object. For a disc, $S_{obj}=\pi D^{2}/4$ while for the experimental cylinder
with a tail, one has $S_{obj}=\pi D^{2}/4+S_{tail}\approx 2D^{2}$ for the
geometry of tail used here. Let’s first consider the results that can be
described by a constant contact time $t_{C}V_{t}/D=(6\pm 3)\times 10^{-2}$.
Following eq. (14), this translates into

$e_{m1}\simeq\left(1-\frac{St_{c}}{St}\right)\exp\left[-\frac{\pi}{6\sqrt{2}}\left(\frac{St_{c}}{\gamma}\right)^{3}\frac{(6\pm
3)\,10^{-2}}{St}\right]$ (15)

The green solid lines in (a,b) describing this model, with $St/St_{c}=75$ and
$\gamma=2$, offer a good description of the experimental results for values of
$St/St_{c}$ large enough, with a strong mismatch for values between $1$ and
$20$. Accounting for a variations of this apparent contact time with the
Stokes number allows to better capture the experimental data. Indeed, the cyan
lines are associated with the time scale of bouncing described by
$t_{C}V_{t}/D=(8\pm 4)\,10^{-3}\sqrt{St/St_{c}-1}$, which implies for the
coefficient of restitution

$e_{m2}\simeq\left(1-\frac{St_{c}}{St}\right)\exp\left[-\frac{\pi}{6\sqrt{2}}\left(\frac{St_{c}}{\gamma}\right)^{3}\frac{(8\pm
4)\,10^{-3}\sqrt{St/St_{c}-1}}{St}\right]\,.$ (16)

Like for the green model, we considered $St_{c}=75$ and $\gamma=2$ to plot the
cyan model as well. The shaded areas indicate the influence of the
uncertainties in estimating the fitting parameters for a value of $St_{c}=75$.
In appendix 7.4, figures 12(a) and (b) show the same model with different
$St_{c}$, keeping $\gamma=2$. The values 50 and 100 were used respectively in
panels (a) and (b). Comparison with figure 10(a) suggests that using
$St_{c}=75$ leads to a better matching between the model and our data.

In summary, the cyan and green lines allow to describe the features of the
rebound in two different limits: the low $St/St_{c}$ limit where 1-contact
with relatively short dimensionless time and sometimes cavitation are
observed, and the high $St/St_{c}$ where rather collisions with 2-contact
occur due to cylinder pitching. Before ending this section, we discuss the
blue lines in figure 10 that correspond to the model of Izard et al. [20] (eq.
(12) and eq. (22)). According to this model, the elasticity and thus the
contact time is variable, assuming a balance between the pressure and elastic
stress during contact. The contact time resulting from this assumption (blue
line in fig. 10(b)) increases with the Stokes number, but it under predicts
the range of measured of contact time. In a consistent way, the coefficient of
restitution obtained from this model (blue line in fig. 10a) follows the same
S-shape but over predicts the experimental results. Thus, while the simple
model of [20] captures qualitatively the features of the rebound and is closer
to the experimental results than the rigid-body model, it fails to predict the
wide range of measured $e=f(St/St_{c})$ that seems to be associated with more
complex apparent elasticity during contact. The latter ingredient is therefore
essential to allow a general description of the rebound of non-spherical
objects.

To conclude this work, we have shown through figure 10 that using a critical
Stokes number and contact time scale (parameterized with $\eta$ and $\alpha$)
allows to capture most of the characteristics of the bouncing: coefficient of
restitution and contact time scale. Even if a perfect matching is not obtained
for both characteristics, its evolution with Stokes and its range of
dispersion are quite consistent in view of the difficulty to measure such
quantities in experiments. Note again that mechanisms at play during
experimental contact are numerous: pitching leading to multiple contact,
cavitation, among other, for which the implementation of a contact elasticity
in models remain an empirical parameterization. A finer description requires
larger database, which is left to future investigation.

##  6 Acknowledgments

Computational resources were provided by the computing meso-centre CALMIP
under project no. P1002. The research federation FERMaT is acknowledged for
the optical measurement facilities. AAC thanks CONACYT AND CIC-UMSNH for their
financial support.

##  7 Appendix

###  7.1 Lubrication force between a cylinder and a wall

Figure 11: Configuration used for the calculation of the viscous lubrication
force in the gap between an infinite circular cylinder (of radius R) and a
rigid plane.

Eq. 6 was obtained in the frame of the thin lubrication approximation [29].
The cylinder is assumed to approach the wall with a given normal velocity
$\textbf{U}_{c}$. The coordinates $x$ and $y$ are equal to 0 at the center
line at the wall (see the geometry in figure 11). We assume that the problem
is quasi-steady (the time variation of the flow velocity u and pressure $p$ is
small compared to their spatial variation). The separation gap $h$ varies in
time due to the cylinder motion according $\frac{\partial h}{\partial
t}=U_{c,y}$, where we recall that the separation gap is not uniform along x,
i.e. $h(x,t)=b(t)+\frac{R}{2}\frac{x^{2}}{R^{2}}$ at the leading order, with
$b$ denoting the minimum height. In the frame of the thin gap lubrication
assumptions, the flow in the gap is dominantly along the $x$ direction, while
the pressure does not vary significantly along the $y$ direction. After
finding the flow velocity profile and assuming a constant flow rate along $x$,
one obtains the following relation between the gap $h$ and the pressure $p$.

$\frac{\partial h}{\partial
t}=\frac{d}{dx}\left(\frac{h^{3}}{12\mu}\frac{dp}{dx}\right)$ (17)

This leads to a differential equation for the pressure, which after
integration along $x$ becomes:

$p=12\mu\frac{\partial h}{\partial
t}\int_{0}^{x}\frac{x}{\left(b+\frac{x^{2}}{2R}\right)^{3}}dx$ (18)

We assume that at a certain distance $R_{0}$ from the center line
($0<R_{0}<R$), the pressure tends toward the outer pressure $p\rightarrow
p_{0}$ and the lubrication assumptions do not hold anymore. The pressure
profile along $x$ direction becomes then:

$p(x)-p_{0}=12\mu\frac{\partial h}{\partial
t}\left[\frac{R}{2\left(b+\frac{R_{0}^{2}}{2R}\right)^{2}}-\frac{R}{2\left(b+\frac{x^{2}}{2R}\right)^{2}}\right]$
(19)

The wall-normal lubrication force can be obtained by integrating the pressure
along the $x$ direction, and assuming small gaps
$\epsilon=\frac{b}{R}\rightarrow 0$:

$F_{y}=-12\pi\mu\frac{\partial h}{\partial
t}\left(\frac{1}{2\epsilon}\right)^{3/2}$ (20)

###  7.2 Calculation of the collision-to-terminal velocity ratio

Integration of (5)-(6) from a distance where the cylinder starts to decelerate
($V\approx V_{t}$ at $\zeta\approx R$) to the wall where collision occurs
($V=V_{c}$ at $\zeta=\eta R$) leads to
$m^{*}\int_{V_{t}}^{V_{c}}dV=-\frac{6}{\sqrt{2}}\pi\mu R^{3/2}\int_{R}^{\eta
R}\frac{d\zeta}{\zeta^{3/2}}$. As the roughness $\eta R$ is very small
compared to the cylinder radius $R$ the equality then approximates to
$m^{*}(V_{c}-V_{t})\simeq\frac{12}{\sqrt{2}}\pi\mu
R\left(\frac{1}{\eta}\right)^{1/2}$, or similarly, $\frac{V_{c}}{V_{t}}\simeq
1-\frac{8}{3\sqrt{2}}\frac{1}{St}\left(\frac{1}{\eta}\right)^{1/2}$. At this
stage, we can define the critical Stokes number, that corresponds to the
minimum inertia required for the cylinder to reach the wall with a finite
velocity that allows the rebound. This critical Stokes number corresponds to
$V_{c}=0^{+}$, which leads to
$St_{c}=\frac{8}{3\sqrt{2}}\left(\frac{1}{\eta}\right)^{1/2}$. Therefore the
collision-to-terminal velocity ratio is given by:

$\frac{V_{c}}{V_{t}}=1-\frac{St_{c}}{St}.$ (21)

###  7.3 Estimation of the contact time from simple scaling analysis

Above the rebound onset, the contact time corresponding to the damped
oscillator harmonic can be approximated by
$\tau=\pi\sqrt{\frac{m^{*}}{k_{n}}}$, as explained below eq. 9. Yet, the
assumption of equilibrium between elastic and viscous lubrication forces
during the contact stage leads to $\lambda|V_{c}|\approx k_{n}\eta R$. Here,
$\lambda$ corresponds to the coefficient of the lubrication force as defined
in 8, the distance between the particle surface and the wall is approximated
by $\eta R$. After replacing the contact velocity $V_{c}$ from eq. 7, one can
estimate the ratio $m^{*}/k_{n}$, and thus obtain the contact time, that we
will call here $\tau_{Iz}$:

$\frac{\tau_{Iz}V_{t}}{D}=\frac{\pi}{2}\left(\frac{3}{\sqrt{2}}\frac{\rho_{p}}{\rho^{*}}\frac{St}{1-St_{c}/St}\right)^{1/2}\eta^{5/4}=\frac{\pi}{2\Gamma}\eta^{5/4}\,,$
(22)

with $1/\Gamma=\left(3\rho_{p}St/\sqrt{2}\rho^{*}(1-St_{c}/St)\right)^{1/2}$

###  7.4 Influence of $St_{c}$ on models $e_{m1}$ and $e_{m2}$

In Figure 12, we compare the influence of the value of $St_{c}$ on the models
for the coefficient of restitution, $e_{m1}$ and $e_{m2}$, described in
equations (15) and (16) respectively.

(a) (b)

Figure 12: Same as in Fig.10(a) but with different values of the threshold for
rebounds in models $e_{m1}$ and $e_{m2}$ with $\gamma=2$, (a) $St_{c}=50$ and
(b) $St_{c}=100$.

## References

  * [1] H. Brenner.  The slow motion of a sphere through a viscous fluid towards a plane surface.  Chemical engineering science, 16(3-4):242–251, 1961.
  * [2] G. Barnocky and R. H. Davis.  Elastohydrodynamic collision and rebound of spheres: experimental verification.  Phys. Fluids, 31(6):1324–1329, 1988.
  * [3] G. G. Joseph, R. Zenit, M. L. Hunt, and A. M. Rosenwinkel.  Particle–wall collisions in a viscous fluid.  J. Fluid Mech., 433:329–346, 2001.
  * [4] P. Gondret, M. Lance, and L. Petit.  Bouncing motion of spherical particles in fluids.  Phys. Fluids, 14(2):643–652, 2002.
  * [5] T. Chastel, P. Gondret, and A. Mongruel.  Texture-driven elastohydrodynamic bouncing.  J. Fluid Mech., 805:577–590, 2016.
  * [6] S. K. Birwa, G. Rajalakshmi, R. Govindarajan, and N. Menon.  Solid-on-solid contact in a sphere-wall collision in a viscous fluid.  Phys. Rev. Fluids, 3(4):044302, 2018.
  * [7] A. Mongruel and P. Gondret.  Viscous dissipation in the collision between a sphere and a textured wall.  J. Fluid Mech., 896, 2020.
  * [8] R. H. Davis, J.-M. Serayssol, and E. J. Hinch.  The elastohydrodynamic collision of two spheres.  J. Fluid Mech., 163:479–497, 1986.
  * [9] G. Lian, M. J. Adams, and C. Thornton.  Elastohydrodynamic collisions of solid spheres.  J. Fluid Mech., 311:141–152, 1996.
  * [10] C. J. Cawthorn and N. J. Balmforth.  Contact in a viscous fluid. part 1. a falling wedge.  J. fluid Mech., 646:327–338, 2010.
  * [11] F.-L. Yang and M. L. Hunt.  A mixed contact model for an immersed collision between two solid surfaces.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 366(1873):2205–2218, 2008.
  * [12] C. M. Donahue, W. M. Brewer, R. H. Davis, and C. M. Hrenya.  Agglomeration and de-agglomeration of rotating wet doublets.  J. Fluid Mech., 708:128–148, 2012.
  * [13] R. H. Davis.  Simultaneous and sequential collisions of three wetted spheres.  J. Fluid Mech., 881:983–1009, 2019.
  * [14] G. Mollon.  A multibody meshfree strategy for the simulation of highly deformable granular materials.  International Journal for Numerical Methods in Engineering, 108(12):1477–1497, 2016.
  * [15] A. M. Ardekani and R. H. Rangel.  Numerical investigation of particle–particle and particle–wall collisions in a viscous fluid.  J. Fluid Mech., 596:437–466, 2008.
  * [16] Z.-G. Feng, E. E. Michaelides, and S. Mao.  A three-dimensional resolved discrete particle method for studying particle-wall collision in a viscous fluid.  J.Fluids Eng., 132(9):091302, 2010.
  * [17] X. Li, M. L. Hunt, and T. Colonius.  A contact model for normal immersed collisions between a particle and a wall.  J. Fluid Mech., 691:123–145, 2012.
  * [18] T. Kempe and J. Fröhlich.  Collision modelling for the interface-resolved simulation of spherical particles in viscous fluids.  J. Fluid Mech., 709:445–489, 2012.
  * [19] JC Brändle de Motta, W-P Breugem, Bertrand Gazanion, J-L Estivalezes, Stéphane Vincent, and Eric Climent.  Numerical modelling of finite-size particle collisions in a viscous fluid.  Phys. Fluids, 25(8):083302, 2013.
  * [20] E. Izard, T. Bonometti, and L. Lacaze.  Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid.  J. Fluid Mech., 747:pp–422, 2014.
  * [21] Q. Li, M. Abbas, and J. F. Morris.  Particle approach to a stagnation point at a wall: Viscous damping and collision dynamics.  Phys. Rev. Fluids, 5(10):104301, 2020.
  * [22] A. Wachs, M. Uhlmann, J. Derksen, and D. P. Huet.  Modeling of short-range interactions between both spherical and non-spherical rigid particles.  In Modeling Approaches and Computational Methods for Particle-Laden Turbulent Flows, pages 217–264. Elsevier, 2023.
  * [23] L. Lacaze, J. Bouteloup, B. Fry, and E. Izard.  Immersed granular collapse: from viscous to free-fall unsteady granular flows.  J. Fluid Mech., 912:A15, 2021.
  * [24] Anthony A. Atchley and Andrea Prosperetti.  The crevice model of bubble nucleation.  The Journal of the Acoustical Society of America, 86(3):1065–1084, 09 1989.
  * [25] G. G. Joseph and M. L. Hunt.  Oblique particle–wall collisions in a liquid.  J. Fluid Mech., 510:71–93, 2004.
  * [26] Q. Li, M. Abbas, J. F. Morris, E. Climent, and J. Magnaudet.  Near-wall dynamics of a neutrally buoyant spherical particle in an axisymmetric stagnation point flow.  J. Fluid Mech., 892:A32:1–29, 2020.
  * [27] R. Zenit, M. L. Hunt, and C. E. Brennen.  Collisional particle pressure measurements in solid-liquid flows.  J. Fluid Mech., 353:261–283, 1997.
  * [28] W. Goldsmit.  The theory and physical behavior of colliding solids.  Edward Arnold Publishers, London, 1960.
  * [29] L. G. Leal.  Advanced transport phenomena: fluid mechanics and convective transport processes, volume 7.  Cambridge University Press, 2007.


[a,b,c,d]Rainer J. Fries

# Excited Hadron Channels in Hadronization

Jacob Purcell  Michael Kordell II  Che-Ming Ko

###### Abstract

The proper treatment of hadronic resonances plays an important role in many
aspects of heavy ion collisions. This is expected to be the case also for
hadronization, due to the large degeneracies of excited states, and the
abundant production of hadrons from their decays. We first show how a
comprehensive treatment of excited meson states can be incorporated into quark
recombination, and in extension, into Hybrid Hadronization. We then discuss
the quantum mechanics of forming excited states, utilizing the Wigner
distribution functions of angular momentum eigenstates of isotropic 3-D
harmonic oscillators. We further describe how resonance decays can be handled,
based on a set of minimal assumptions, by creating an extension of hadron
decays in PYTHIA 8. Finally, we present first results by simulating
$e^{+}+e^{-}$ collisions using PYTHIA and Hybrid Hadronization with excited
mesons up to orbital angular momentum $L=4$ and radial quantum number 2. We
find that states up to $L=2$ are produced profusely by quark recombination.

Hadronization of partons is a longstanding problem for Monte Carlo (MC) event
generators. The process of forming bound states of quarks and gluons can not
be described from first principles. Instead, several models have been
developed over the years to successfully describe certain aspects of
hadronization. Among them are the string fragmentation model, that is for
example deployed in the PYTHIA event generator [1, 2], and the quark
recombination model [3, 4, 5, 6]. String fragmentation is based on the idea of
QCD strings forming between color charges at large distances, and has been
successfully applied to all kind of "small" collision systems, i.e. those not
involving nuclei. On the other hand, quark recombination has had success
describing certain aspects of hadronization in nuclear collisions, in
particular large baryon/meson ratios and the constituent-quark-number scaling
of elliptic flow. The idea behind Hybrid Hadronization is to combine these two
models so a comprehensive and consistent description of hadronization in all
collision systems at a broad range of collision energies can be achieved [7,
8].

The overarching idea is that systems of partons that are ready to hadronize
are first allowed to recombine by sampling the recombination probabilities of
all quark-antiquark pairs and all quark and anitquark triplets in the system.
Gluons are assumed to have decayed into quark-antiquark octet pairs for this
step. The recombination probabilities are computed in a phase-space formalism
briefly sketched below. Some partons, preferentially those close in phase
space, thus recombine directly into hadrons. All remnant partons are then
assumed to be connected by strings and these string are subjected to a string
fragmentation model, in our case the one implemented in PYTHIA 8. If the
original system of partons has color tags assigned, e.g. in output from PYTHIA
8 for $e^{+}+e^{-}$ or $p+p$ collisions, these color tags are used in the
computation of recombination probabilites, and they are again used to form the
remnant strings, preserving color flow in the parton system.

The focus of this work is the implementation of proper physical hadronic
resonances in the recombination process. In the orignal work [7] excited
hadrons were not mapped onto the proper physical states in the Particle Data
Book [10] since the recombination probabilities into eigenstates of orbital
angular momentum in the Wigner phase-space formalism were not known. Our study
of angular momentum eigenstates of the 3-D isotropic harmonic oscillator [9]
remedies this shortcoming. We briefly discuss the work done in [9] and then
present first results from an application to $e^{+}+e^{-}$ collisions.

We assume the potential between color singlet quark-antiquark pairs to be be
modelled by a 3-D isotropic harmonic oscillator. Our work restricts itself to
mesons for now for simplicity. The widths of the potentials can be fixed to
data by computing the squared charge radii for stable mesons and comparing
those results to measured values of $\langle r^{2}\rangle$, as laid out in
[7]. The inverse length scale of the harmonic oscillator is denoted by $\nu$
in the following. We compute the Wigner distributions
$W_{kl}(\mathbf{r},\mathbf{q})$ of eigenstates of the potential with radial
quantum number $k$ and orbital angular momentum quantum number $l$. The
magnetic quantum number $m$ is averaged over since we do not wish to consider
the polarization of hadrons. In that case the distributions only depend on the
magnitudes of the position and momentum vectors $\mathbf{r}$ and $\mathbf{q}$,
and the angle $\theta$ between those vectors. In our work we reduce the
complicated problem of 3-D phase-space distributions to the known Wigner
distributions of the 1-D harmonic oscillator [11]. These distributions have
been computed before in [12] using a very different approach. The lowest
energy distributions are

$\displaystyle W_{00}$
$\displaystyle=\frac{1}{\pi^{3}\hbar^{3}}e^{-\frac{q^{2}}{\hbar^{2}\nu^{2}}-\nu^{2}r^{2}}\,,$
(1) $\displaystyle W_{01}$
$\displaystyle=W_{00}\left(-1+\frac{2}{3}\nu^{2}r^{2}+\frac{2}{3}\frac{q^{2}}{\hbar^{2}\nu^{2}}\right)\,,$
(2) $\displaystyle W_{10}$
$\displaystyle=W_{00}\left(1+\frac{2}{3}\nu^{4}r^{4}-\frac{4}{3}\nu^{2}r^{2}-\frac{4}{3}\frac{r^{2}q^{2}}{\hbar^{2}}+\frac{8}{3}\frac{(\mathbf{r}\cdot\mathbf{q})^{2}}{\hbar^{2}}-\frac{4}{3}\frac{q^{2}}{\hbar^{2}\nu^{2}}+\frac{2}{3}\frac{q^{4}}{\hbar^{4}\nu^{4}}\right)\,.$
(3)

The distributions for $k=1$, $l=0$ are also visualized in Fig. 1.

Figure 1:  Phase-space distributions $W_{10}$ for quantum numbers $k=1$, $l=0$
as functions of $r=|\mathbf{r}|$ and $q=|\mathbf{q}|$. The white lines
indicate nodes where $W_{10}=0$. Distributions are shown for several values of
the angle $\theta$ given by $\cos\theta=\mathbf{r}\cdot\mathbf{q}/rq$.

We procede by computing the probabilites for two Gaussian wave packets with
given widths $\delta$, representing the quark and antiquark, to coalesce into
a bound state given by the Wigner distributions above. The probabilities
$P_{kl}(\mathbf{x},\mathbf{p})$ depend on the relative distance vectors
$\mathbf{x}$ and $\mathbf{p}$ of the wave packets in position and momentum
space, respectively. Here the magnetic quantum number $m$ is summed over to
capture all polarization states. The probabilities for the lowest energy
states are

$\displaystyle\mathcal{P}_{00}$ $\displaystyle=e^{-u}\,,$ (4)
$\displaystyle\mathcal{P}_{01}$ $\displaystyle=e^{-u}u\,,$ (5)
$\displaystyle\mathcal{P}_{10}$
$\displaystyle=\frac{1}{2}e^{-u}\left(\frac{1}{3}u^{2}-\frac{1}{3}t\right)\,$
(6)

for the simplest case $2\delta\nu=1$. For a discussion of other ratios between
the length scales see [9]. The distrubutions have been written in terms of
variables $u=\nu^{2}x^{2}/2+p^{2}/(2\hbar^{2}\nu^{2})$ and
$t=(\mathbf{x}\times\mathbf{p})^{2}/\hbar^{2}$, which have straightforward
physical interpretations as the dimensionless squared distance of the initial
partons in phase space, and their dimensionless squared angular momentum,
respectively. It is then interesting to analyze the mapping of initial angular
momentum $L^{2}$ in the parton system onto the choice of angular momentum
quantum numbers $l$, in particular if states of several different $l$ with
degenerate energy quantum number $n=2k+l$ are available [9]. Maps of the
coalescence probabilities $P_{10}$ are shown in Fig. 2

Figure 2:  Coalescence probabilities $P_{10}$ for quantum numbers $k=1$, $l=0$
for two Gaussian wave packets interacting through an isotropic 3-D harmonic
oscillator potential. Probabilities are shown as functions of relative
coordinates $x=|\mathbf{x}|$ and $p=|\mathbf{p}|$ for several values of the
angle $\cos\theta=\mathbf{x}\cdot\mathbf{p}/xp$ between the vectors.

We now turn to the implementation of this formalism for the recombination of
mesons within the Hybrid Hadronization framework. The phase-space
probabilities $P_{kl}$ are supplemented with probablities $P_{s}$ to form spin
single or triplet states and the probability $P_{c}$ to form color singlets.
We have implemented light excited mesons up to $n=4$. This means we include
mesons up to $G$-wave for the radial ground state $k=0$, and up to $D$-wave
for a single radial excitation $k=1$. This set of states exceeds the states
listed as confirmed by the Particle Data Group [10]. To preserve unitarity we
opt to keep the experimentally unconfirmed states. However, we are in need of
the masses and strong decay branching ratios of these states. Missing masses
are estimated by using the phenomenologically established scaling of squared
masses with radial and spin quantum numbers, thus using linear interpolation
and extrapolation. We consider decays into the simplest possible sets of
stable hadrons plus additional pions, up to a maximum of five hadrons in
total, which are allowed by the quantum numbers. Branching ratios are
determined using phase-space weights and isospin algebra. As a random example,
we pick the $1^{3}F_{2}$ isospin triplet state which would be known as
$a_{2}(1918)$. We estimate its mass to be $1.918$ GeV and allow its decay into
two 3-pion states (63.3%) and three 5-pion states (36.7%). The necessary data
on the full set of meson states considered by us is collected in a XML
particle data file following PYTHIA 8 standards, which enables Pythia to read
in and decay the excited states [13].

With the formalism set we look at a first example. We use PYTHIA 8 to generate
$e^{+}e^{-}$ di-jet events at $\sqrt{s}=91.2$ GeV. The partons after final-
state showers are extracted and fed into a standalone version of the Hybrid
Hadronization code that supports the generation of the new physical meson
resonances. Baryons are still treated as described in [7]. Hybrid
Hadronization calls another instance of PYTHIA 8 to hand over remnant strings
to fragment and includes hadrons from recombination. Both recombination and
fragmentation hadrons can then decay according to settings in PYTHIA 8.

We focus on one important preliminary result from this study. If one analyzes
the quantum numbers of hadrons from recombination before any secondary decays.
One finds the result shown in Fig. 3. The three panels depict the relative
abundances of mesons as functions of radial and orbital angular momentum
quantum numbers $k$ and $l$, as well as the total angular momentum quantum
number $j$ which comes from the addition of orbital and spin angular momentum
operators, $\mathbf{J}=\mathbf{L}+\mathbf{S}$. The distribution of radial
quantum numbers is strongly peaked at the ground state $k=0$ with the next
excited stated adding about one third of that strength. On the other hand,
$P$\- and $D$-wave states ($l=1$ and $2$, respectively) are as numerous as
$S$-wave states, showing the importance of orbital angular momentum
excitations in this jet systems. Even $F$-wave states are created with roughly
half the strength of $S$-waves. As a result, vector and tensor mesons are, by
far, the dominant channels in recombination, with $j=3$ mesons still being
more than doubly as numerous as spinless mesons.

Figure 3:  Relative yields of mesons of different radial quantum number $k$
(left panel), orbital angular momentum quantum number $l$ (center panel) and
total angular momentum quantum number $j$ (right panel) from recombination in
$e^{+}e^{-}$ collisions at $\sqrt{s}=91.2$ GeV using PYTHIA 8 and Hybrid
Hadronization, before decays of excited states.

In summary, we have developed a formalism that allows for the inclusion of
physical excited mesons into the quark recombination model and, in extension,
the Hybrid Hadronization model. To this end we have computed the probabilties
for the coalescence of Gaussian wave packets into angular momentum eigenstates
of 3-D isotropic oscillator potentials, approximating the dynamics of quark-
antiquark systems at distances which are not too large. We have utilized a
phase-space formalism which allows us to compute quark-antiquark coalescence
probabilities from the mean distance of the partons in coordinate and momentum
space. We have compiled masses and branching ratios for strong decays for
mesons up to $G$-waves, using estimates for those not found in the Particle
Data Book. We find that in $e^{+}e^{-}$ collisions at $\sqrt{s}=91.2$ GeV,
using PYTHIA 8 and Hybrid Hadronization, mesons states with excited orbital
angular momenta dominate recombination. In the future, we plan to extend our
formalism to heavy flavor mesons and baryons.

This work was supported by the U.S. National Science Foundation under awards
1812431 and 2111568, and under award 2004571 through a subcontract with Wayne
State University. In addition, this work was supported by the U.S. Department
of Energy under Award No. DE-SC0015266. RJF would like to thank the ExtreMe
Matter Institute EMMI at the GSI Helmholtzzentrum für Schwerionenforschung,
Darmstadt, for support and the Institute of Theoretical Physics at the
University of Frankfurt and the Frankfurt Insitute for Advanced Studies for
their hopitality.

## References

  * [1] B. Andersson, G. Gustafson, G. Ingelman and T. Sjostrand, Phys. Rept. 97, 31-145 (1983) doi:10.1016/0370-1573(83)90080-7
  * [2] C. Bierlich, S. Chakraborty, N. Desai, L. Gellersen, I. Helenius, P. Ilten, L. Lönnblad, S. Mrenna, S. Prestel and C. T. Preuss, et al., doi:10.21468/SciPostPhysCodeb.8, [arXiv:2203.11601 [hep-ph]].
  * [3] R. J. Fries, B. Müller, C. Nonaka and S. A. Bass, Phys. Rev. Lett. 90, 202303 (2003), doi:10.1103/PhysRevLett.90.202303, [arXiv:nucl-th/0301087 [nucl-th]].
  * [4] R. J. Fries, B. Müller, C. Nonaka and S. A. Bass, Phys. Rev. C 68, 044902 (2003), doi:10.1103/PhysRevC.68.044902, [arXiv:nucl-th/0306027 [nucl-th]].
  * [5] V. Greco, C. M. Ko and P. Levai, Phys. Rev. Lett. 90, 202302 (2003) doi:10.1103/PhysRevLett.90.202302 [arXiv:nucl-th/0301093 [nucl-th]].
  * [6] V. Greco, C. M. Ko and P. Levai, Phys. Rev. C 68, 034904 (2003) doi:10.1103/PhysRevC.68.034904 [arXiv:nucl-th/0305024 [nucl-th]].
  * [7] K. C. Han, R. J. Fries and C. M. Ko, Phys. Rev. C 93, no.4, 045207 (2016), doi:10.1103/PhysRevC.93.045207, [arXiv:1601.00708 [nucl-th]].
  * [8] R. J. Fries and M. Kordell, PoS HardProbes2018, 046 (2019) doi:10.22323/1.345.0046 [arXiv:1901.08157 [nucl-th]].
  * [9] M. Kordell, II, R. J. Fries and C. M. Ko, Annals Phys. 443, 168960 (2022), doi:10.1016/j.aop.2022.168960, [arXiv:2112.12269 [quant-ph]].
  * [10] P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020) doi:10.1093/ptep/ptaa104
  * [11] T. Curtright, T. Uematsu, C. K. Zachos, J. Math. Phys. 42, 2396 (2001). doi:10.1063/1.1366327, arXiv:hep-th/0011137.
  * [12] S. Shlomo and M. Prakash, Nucl. Phys. A 357, 157-170 (1981), doi:10.1016/0375-9474(81)90631-X.
  * [13] R. J. Fries, J. Purcell, M. Kordell, II, and C. M. Ko, in preparation.