[1]Ruben Tous

[1]Department of Computer Architecture, Universitat Politècnica de Catalunya
(UPC), Jordi Girona, 1-3, Barcelona, 08034, Spain

# Lester: rotoscope animation through video object segmentation and tracking

<EMAIL_ADDRESS>*

###### Abstract

This article introduces Lester, a novel method to automatically synthetise
retro-style 2D animations from videos. The method approaches the challenge
mainly as an object segmentation and tracking problem. Video frames are
processed with the Segment Anything Model (SAM) and the resulting masks are
tracked through subsequent frames with DeAOT, a method of hierarchical
propagation for semi-supervised video object segmentation. The geometry of the
masks’ contours is simplified with the Douglas-Peucker algorithm. Finally,
facial traits, pixelation and a basic shadow effect can be optionally added.
The results show that the method exhibits an excellent temporal consistency
and can correctly process videos with different poses and appearances, dynamic
shots, partial shots and diverse backgrounds. The proposed method provides a
more simple and deterministic approach than diffusion models based video-to-
video translation pipelines, which suffer from temporal consistency problems
and do not cope well with pixelated and schematic outputs. The method is also
much most practical than techniques based on 3D human pose estimation, which
require custom handcrafted 3D models and are very limited with respect to the
type of scenes they can process.

###### keywords:

Animation, Rotoscoping, Segmentation, Computer Graphics, video, Artificial
Intelligence

##  1 Introduction

In industries such as cinema or video games, rotoscoping is employed to
leverage the information from live-action footage for the creation of 2D
animations. This technique contributed to the emergence of a unique visual
aesthetic in certain iconic video games from the 80s and 90s, such as “Prince
of Persia” (1989), “Another World” (1991) or “Flashback” (1992), and it
continues to be employed in modern titles like “Lunark” (2023). However, as in
the case of Lunark, this technique still requires manual work at the frame
level, making it very expensive. This article describes a novel method, called
Lester, to automate this process. A user enters a video and a color palette
and the system automatically generates a 2D retro-style animation in a format
directly usable in a video game or in an animated movie.

(a)

(b)

Figure 1: Example synthetised animation frame (right) and its related original
frame (left, face has been anonymized) from an example input video.

The method approaches the challenge mainly as an object segmentation problem.
Given a frame depicting a person, a mask for each different visual trait
(hair, skin, clothes, etc.) is generated with the Segment Anything Model (SAM)
[1]. To univocally identify the masks in all the frames, a video object
tracking algorithm is employed using the approach of ([2]). The masks are
purged to solve segmentation errors and their geometry is simplified with the
Douglas-Peucker algorithm [3] to achieve the desired visual style. Finally,
facial traits, pixelation and a basic shadow effect can be optionally added.

This being essentially a specific case of unsupervised video-to-video
translation, a different strategy could involve employing conditional
generative models, like diffusion models. Yet, as we will explore
subsequently, the proposed method proves to be a superior alternative due to
its enhanced efficacy in handling the specific visual style and addressing
temporal consistency issues.

Figure 2: Overview of the system’s workflow (inputs on the left, output on the
right)

##  2 Related Work

Max Fleischer’s invention of the rotoscope in 1915 transformed the ability to
create animated cartoons by using still or moving images. Variants of this
technique have been used to this day in animated films (e.g. “A Scanner
Darkly”, (2006)), television animated series (e.g. “Undone” (2019)) and video
games (e.g. “Prince of Persia” (1989), “Another World” (1991), “Flashback”
(1992), “Lunark” (2023)). Contemporary productions integrate image processing
algorithms to facilitate the task, but they are mainly based on costly manual
work. Achieving full automation with a polished professional outcome (uniform
style, spatial and temporal consistency) remains unfeasible at present.
However, recent advances in Artificial Intelligence (AI) opens the possibility
to attempt advancing to fully automation.

To the best of our knowledge, there are two main strategies to address the
problem. On the one hand, the challenge can be addressed with conditional
generative models, prevalent in image translation tasks nowadays. The first
methods (e.g. CartoonGAN [4], AnimeGAN[5] or the work in [6]) used conditional
generative adversarial networks (GAN) to transform photos of real-world scenes
into cartoon style images. These methods have recently been surpassed by those
based on diffusion models such as Stable Diffusion[7]. Although most of these
methods have been designed for text-to-image generation, they have been
successfully adapted for image-to-image translation using techniques such as
SDEdit [8] or ControlNet [9]. Pure image-to-image translation methods with
diffusion models also exists (e.g. SR3 [10] and Palette [11]), but with a
supervised setup that can only be applied to situations in which a groundtruth
of pairs of images (original, target) such as super-resolution or
colorization. In any case, even if we have a good conditional generative
model, here is still the challenge of applying it to the frames of a video
while preserving temporal consistency. Latest attempts (e.g. [12]) combine
tricks (e.g. tiled inputs) to uniformize the output of the generative model
with video-guided style interpolation techniques (e.g. EbSynth [13]). The
results of these methods have improved greatly, and for some styles the
results are so good that they are already used in independent productions.
However, for most digital art styles these methods impose an undesirable
trade-off between style fidelity and temporal consistency.

On the other hand, a completely different approach consists on applying a 3D
human pose estimation algorithm to track the 3D pose through the different
frames. The 3D pose is then applied to a 3D character model with the desired
aspect (hair, clothes, etc.). Typically this 3D character model is created by
hand with a specialised software. Finally, a proper shader and rendering
configuration are setup to obtain the desired visual style. This is the
approach taken in [14], [15], [16], [17] and [18]. This approach has many
drawbacks. First, it depends on the state-of-art of 3D human pose estimation,
which has advanced a lot but still fails in many situations (low resolution,
truncations, uncommon poses) mainly due to the difficulty to obtain enough
groundtruth data. Second, it requires to design a proper 3D model for each
different character. Third, the 3D pose is based on a given model, typically
SMPL. To apply it to a different model it requires a pose retargetting step.
The joint rotations can be easily retargetted but the SMPL blend shapes not,
resulting in unrealistic body deformations. This is why demos of this kind of
works are usually done with robots, as the don’t have deformable parts. An
advantage of this approach with respect to generative models is that it
provides a good temporal consistency, as inconsistencies from the 3D human
pose estimation step can be properly solved with interpolation such as is done
in [17].

In this work we propose a third approach, mainly based on segmentation and
tracking. Unlike the first mentioned strategy (generative models), the
proposed approach is more predictable and provides much better temporal
consistency. Unlike the second strategy mentioned (3D human pose estimation),
the proposed approach is more robust and much more practical since it does not
require having a custom 3D model for each new character. This is not the first
time that this kind of problem has been addressed this way. In [19]
segmentation was used to stylize a portrait video given an artistic example
image obtaining better results than neural style transfer [20], prevalent at
that time.

##  3 Methodology

###  3.1 Overview

The main input to the method is a target video sequence $T$ of a human
performance. There are not general constraints about the traits of the
depicted person (gender, age, clothes, etc.), pose, orientation, camera
movements, background, lighting, occlusions, truncations (partial body shots)
or resolution. However, certain scenarios are known to potentially impact the
results negatively (this will be addressed in Section 4). While the method is
primarily designed for videos featuring a single person, its theoretical
application with multiple individuals will be discussed in Section 4. When
operating without additional inputs, the method defaults to specific settings
for segmentation, colors, and other attributes. However, it is anticipated
that the user might prefer specifying custom configurations, necessitating
supplementary inputs such as segmentation prompts, color palettes, etc., as
detailed in corresponding subsections below. The output of the method is a set
of PNG images showing the resulting animated character and no background
(where the alpha channel is set to zero).

The method consists of three independent processing steps. The main one is the
segmentation and tracking of the different visual traits (hair, skin, clothes,
etc.) of the person featured in the video frames. The second step focuses on
rectifying segmentation errors through mask purging and simplifying their
polygonal contours. The third and final step involves applying customizable
finishing details, including colors, faux shadows, facial features, and
pixelation. Figure 2 shows the overall workflow of the method. Each component
is elaborated upon in the subsequent subsections.

###  3.2 Segmentation and tracking

In order to segment the different visual traits of the character the Segment
Anything Model (SAM) [1] is used. SAM is a promptable segmentation model which
is becoming a foundation model to solve a variety of tasks. Given an input
image with three channels (height x width x 3) and a prompt, that can include
spatial or text information, it generates a mask (height x weight x 3), with
each pixel assigned an RGB color. Each color represents an object label from
the object label space $L=l_{1},l_{2},...$, inferred from the prompt,
corresponding to the predicted object type at the corresponding pixel location
in the input image. The prompt can be any information indicating what to
segment in an image (e.g. a set of foreground/background points, a rough box
or mask or free-form text). SAM consists on a image encoder (based on MAE [21]
a pre-trained Vision Transformer [22]) that computes an image embedding, a
prompt encoder (being the free-form text based on CLIP [23]) that embeds
prompts, and then the two information sources are combined in a mask decoder
(based on Transformer segmentation models [24]) that predicts segmentation
masks.

Applying SAM to all the frames of the input video is not enough as each time
SAM employs different object labels. To do it consistently we use SAM-Track
[2]. SAM results for the first frame are propagated to the second frame, with
the same labels, with DeAOT [25], a method of hierarchical propagation for
semi-supervised video object segmentation. The hierarchical propagation can
gradually propagate information from past frames to the current frame and
transfer the current frame feature from object-agnostic to object-specific.
SAM is applied to subsequent frames to detect new objects, but the
segmentation of already seen objects is carried out by DeAOT.

We have chosen keyword-based prompts to specify the desired segmentation. We
provide some default prompts for common scenarios (e.g. ’hair’, ’skin’,
’shoes’, ’t-shirt’, ’trousers’) but it’s anticipated that in many scenarios,
the user will need to specify an alternative. These prompts are not processed
with the free-text capabilities of SAM because of their limitations, instead
they are processed with Grounding-DINO [26], an open-set object detector,
which translate them into object boxes that are sent to SAM as box prompts.
One problem we have identified is that certain object combinations lead to
incorrect maks when they are specified in the same prompt. To solve this, we
allow specifying multiple prompts that are processed independently. Figure 3
shows the workflow of the segmentation and tracking stage.

Figure 3: Workflow of the segmentation and tracking stage for the first two
frames (inputs on the left, output on the right).

###  3.3 Contours simplification

One of the desired choices for the target visual style is the possibility to
display very simple geometries, i.e. shape contours with short polygonal
chains, in order to mimic the visual aspect of, e.g., “Another World”. To
achieve it, a series of steps are carried out. For each frame, the result of
the segmentation and tracking step is a mask with each pixel assigned an RGB
color representing an object label. The mask is first divided into submasks,
one for each color/label. Each submask is divided, in turn, into contours with
the algorithm from [27]. Each individual contour is represented with a vector
of the vertex coordinates of the boundary line segments. The number of line
segments of each contour is reduced with the Douglas-Peucker algorithm [3].
The algorithm starts by joining the endpoints of the contour with a straight
line. Vertices closer to the line than a given tolerance value $t$ are
removed. Next, the line is split at the point farthest from the reference
line, creating two new lines. The algorithm works recursively until all
vertices within the tolerance are removed. Tolerance $t$ is one of the
optional parameters of the proposed method and let the user control the level
of simplification. As each contour is simplified independently, undesirable
spaces can appear between the geometries of bordering contours that previously
were perfectly aligned. To avoid this, and before the simplification step, all
contours are dilated convolving the submask with a 4x4 kernel. It is worth
mentioning that he simplification algorithm may fail for certain contour
geometries and values of $t$. In that situations (usually very small contours)
the original contour is kept. See Algorithm 1 and Figure 4.

Algorithm 1 Contour simplification algorithm

1:$M$ a list containing submasks, $\alpha$ a minimum contour area, $t$ a
tolerance value for simplification

2:$R$ a list of contours where $R_{i,j}$ is the simplified contour $j$ of
submask $i$

3:

4:function Simplify($M,\alpha,t$)

5: for $i\leftarrow 0$ to $length(M)-1$ do

6: $C_{i}\leftarrow getContours(M_{i})$

7: for $j\leftarrow 0$ to $length(C_{i})-1$ do

8: if $areaOf(C_{i,j})>\alpha$ then

9: $D_{i,j}\leftarrow dilateContour(C_{i,j})$

10: if $contourCanBeSimplified(D_{i,j},t)$ then

11: $R_{i,j}\leftarrow simplifyContour(D_{i,j},t)$

12: else

13: $R_{i,j}\leftarrow D_{i,j}$

14: end if

15: end if

16: end for

17: end for

18: return $R$

19:end function

(a)

(b)

(c)

(d)

(e)

Figure 4: Graphical overview of the contour simplification stage of one frame:
(a) Result of the segmentation and tracking step (b) Contours detection (c)
Contours dilation (d) Contours simplification (e) Resulting colored frame.

###  3.4 Finishing details

The third and final step of the method entails applying customizable finishing
details, granting users flexibility in shaping the visual style of the
outcome. The most important aspect is the color. The method generates a visual
guide with unique identifiers for each color segment, as shown in Figure 5.
One of the inputs to the method is a color palette, which indicates the RGB
values that will be applied to each segment identifier.

(a)

(b)

(c)

(d)

Figure 5: Coloring process: (a) Input frame (b) Segmentation (c) Visual guide
with unique identifiers for each color segment (d) Final frame.

Another finishing detail, which is optional, is the possibility of adding some
schematic-looking facial features. The features are generated from the 68
facial landmarks that can be obtained using the frontal face detector
described in [28], which employs a Histogram of Oriented Gradients (HOG) [29]
combined with a linear classifier, an image pyramid, and a sliding window
detection scheme. Another option is the possibility of adding a very simple
shadow effect by superimposing the same image with different brightness and
displaced horizontally. Finally, a pixelation effect is achieved by resizing
the image down using bilinear interpolation and scaling it back up with
nearest neighbour interpolation. Figure 6 shows how shadows, facial features
and pixelation are progressively applied to an example frame.

(a)

(b)

(c)

(d)

Figure 6: Optional finishing details: (a) Flat method result (b) Shadows (c)
Facial features (d) Pixelation.

##  4 Experiments and results

In order to evaluate the performance of the method, we collected a set of 25
short videos ranging from 5 to 20 seconds in duration. The videos showcase a
lone human subject displaying various clothing styles and engaging in diverse
actions. The videos were obtained from the UCF101 Human Actions dataset [30],
the Fashion Video Dataset [31], the Kinetics dataset [32] and YouTube. The
assembled datased includes both low and hight resolution videos and both in-
the-wild and controlled-setting videos. The dataset and the results of the
experiments can be found in [33]. Figures 7 and 8 shows some of example
results.

The subjective quality assessment of the results is summarized in Table 1. We
report the MOS (mean opinion score described in ITU-T P.800.1 [34]) of the ACR
(absolute category rating described in ITU-T P.910 recommendation [35] scale
from 1=Bad to 5=Excellent) of 3 aspects: shape quality, color consistency and
temporal consistency, plus and overall score. The videos are grouped into
three categories, in-the-wild low resolution (ITW-LR), in-the-wild high
resolution (ITW-HR) and controlled setting high resolution (C-HR).

subset | # | shape | color | temporal | overall
---|---|---|---|---|---
ITW-LRa | 5 | 3.20 | 3.24 | 3.84 | 3.36
ITW-HRb | 14 | 3.81 | 4.03 | 4.33 | 3.99
C-HRc | 6 | 4.70 | 4.60 | 4.67 | 4.67
all |  | 3.90 | 4.01 | 4.31 | 4.02
Table 1: Subjective quality assessment MOS values of the ACR (1=Bad to
5=Excellent) of shape quality, color consistency, temporal consistency and
overall quality. aITW-LR $\rightarrow$ in-the-wild low resolution; bITW-HR
$\rightarrow$ in-the-wild high resolution; cC-HR $\rightarrow$ controlled
setting high resolution

Strengths: The results show that the method is able to synthetize quality
animations from diverse input videos (an average MOE 4.02 over 5). The
performance significantly improves (4.67) on high resolution videos taken in
controlled conditions (videos 2, 11, 12, 13, 14 and 15). The method exhibits a
very good temporal consistency and can correctly process videos with different
poses and appearances, dynamic shots (e.g. video 5), videos with partial shots
(e.g. the medium close-up shot in video 7) and diverse backgrounds.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

(u)

(v)

(w)

(x)

(y)

(a)

(b)

(c)

(d)

(e)

Figure 7: Example results. Odd rows: input videos (faces have been
anonymized); Even rows: method results.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

(u)

(v)

(w)

(x)

(y)

(a)

(b)

(c)

(d)

(e)

Figure 8: Example results. Odd rows: input videos (faces have been
anonymized); Even rows: method results.

Limitations: The performance decreases in videos with low resolution (all
videos from the UCF101 Human Actions dataset), complex backgrounds (videos 17
and 20), low contrast colors (videos 8, 10, 22 and 23), inappropriate lighting
(videos 14 and 25) and occlusions (video 23). Figure 9 shows four example
frames with errors from four different videos (videos 22, 23, 25 and 8 from
the test dataset).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 9: Example failure cases. Top row: input videos (videos 22, 23, 25 and
8 from the test dataset, faces have been anonymized); bottom row: results.

Computational cost: The experiments have been run on a high-end server with a
quadcore Intel i7-3820 at 3.6 GHz with 64 GB of DDR3 RAM, and 4 NVIDIA Tesla
K40 GPU cards with 12 GB of GDDR5 each. For a single frame, step 1
(segmentation and tracking) takes 0.56 seconds, step 2 (contours
simplification) takes 0.43 seconds and step 3 (finishing details, mainly
facial features) takes 1.32 seconds. Overall the method operates at 0.43 FPS
(1 FPS without facial features).

Discussion: The results obtained are very satisfactory and validate the base
hypothesis, that the latest advances in segmentation and tracking make this
technique an effective method for generating rotoscopic animations from
videos. The limitations detected do not invalidate the initial hypothesis
since it is expected that the method will be used on selected or purposely
created videos. Nor does its computational cost since the method is expected
to be used in editing tasks that do not require real-time execution. The
method is primarily designed for videos featuring a single person, but could
be adapted to videos with multiple people by using full-body masks and
multiple parameterizations for segmentation prompts and color palettes.

##  5 Conclusions

In this work we present Lester, a method for automatically synthetise 2D
retro-style animated characters from videos. Rather than addressing the
challenge with an image-to-image generative model, or with human pose
estimation tools, it is addressed mainly as an object segmentation and
tracking problem. Subjective quality assessment shows that the approach is
feasible and that the method can generate quality animations from videos with
different poses and appearances, dynamic shots, partial shots and diverse
backgrounds. Although it is not the objective, it is worth mentioning that the
performance of the method is less consistent in-the-wild, especially against
low resolution videos. But the method is expected to be used on selected or
custom-made videos in semi-controlled environments.

Due to its capabilities, the method constitutes an effective tool that
significantly reduces the cost of generating rotoscopic animations for video
games or animated films. Due to its characteristics, it is more simple and
deterministic than generative models based pipelines and more practical and
flexible than techniques based on human pose estimation.

section*Acknowledgements This work is partially supported by the Spanish
Ministry of Science and Innovation under contract PID2019-107255GB, and by the
SGR programme 2021-SGR-00478 of the Catalan Government.

## Declarations

### Competing interests

The authors declare that they have no competing interests.

### Data availability

The dataset and the results of the experiments are publicly available at
https://github.com/rtous/lester.

## References

  *   * Kirillov et al. [2023] Kirillov, A., Mintun, E., Ravi, N., Mao, H., Rolland, C., Gustafson, L., Xiao, T., Whitehead, S., Berg, A.C., Lo, W.-Y., et al.: Segment anything. arXiv preprint arXiv:2304.02643 (2023)
  * Cheng et al. [2023] Cheng, Y., Li, L., Xu, Y., Li, X., Yang, Z., Wang, W., Yang, Y.: Segment and track anything. arXiv preprint arXiv:2305.06558 (2023)
  * Douglas and Peucker [2011] Douglas, D.H., Peucker, T.K.: Algorithms for the Reduction of the Number of Points Required to Represent a Digitized Line or its Caricature, pp. 15–28 (2011). https://doi.org/10.1002/9780470669488.ch2
  * Chen et al. [2018] Chen, Y., Lai, Y.-K., Liu, Y.-J.: Cartoongan: Generative adversarial networks for photo cartoonization. In: CVPR, pp. 9465–9474. IEEE Computer Society, ??? (2018)
  * Chen et al. [2020] Chen, J., Liu, G., Chen, X.: AnimeGAN: A Novel Lightweight GAN for Photo Animation, pp. 242–256. Artificial Intelligence Algorithms and Applications, ??? (2020). https://doi.org/10.1007/978-981-15-5577-0_18
  * Liu et al. [2022] Liu, Z., Li, L., Jiang, H., Jin, X., Tu, D., Wang, S., Zha, Z.: Unsupervised coherent video cartoonization with perceptual motion consistency. CoRR abs/2204.00795 (2022) https://doi.org/10.48550/arXiv.2204.00795 2204.00795
  * Rombach et al. [2021] Rombach, R., Blattmann, A., Lorenz, D., Esser, P., Ommer, B.: High-Resolution Image Synthesis with Latent Diffusion Models (2021)
  * Meng et al. [2022] Meng, C., He, Y., Song, Y., Song, J., Wu, J., Zhu, J.-Y., Ermon, S.: SDEdit: Guided image synthesis and editing with stochastic differential equations. In: International Conference on Learning Representations (2022)
  * Zhang et al. [2023] Zhang, L., Rao, A., Agrawala, M.: Adding Conditional Control to Text-to-Image Diffusion Models (2023)
  * Saharia et al. [2021] Saharia, C., Ho, J., Chan, W., Salimans, T., Fleet, D.J., Norouzi, M.: Image super-resolution via iterative refinement. arXiv:2104.07636 (2021)
  * Saharia et al. [2022] Saharia, C., Chan, W., Chang, H., Lee, C., Ho, J., Salimans, T., Fleet, D., Norouzi, M.: Palette: Image-to-image diffusion models. In: ACM SIGGRAPH 2022 Conference Proceedings. SIGGRAPH ’22. Association for Computing Machinery, New York, NY, USA (2022). https://doi.org/10.1145/3528233.3530757 . https://doi.org/10.1145/3528233.3530757
  * Yang et al. [2023] Yang, S., Zhou, Y., Liu, Z., , Loy, C.C.: Rerender a video: Zero-shot text-guided video-to-video translation. In: ACM SIGGRAPH Asia Conference Proceedings (2023)
  * Jamriska [2018] Jamriska, O.: Ebsynth: Fast Example-based Image Synthesis and Style Transfer. GitHub (2018)
  * Liang et al. [2005] Liang, D., Liu, Y., Huang, Q., Zhu, G., Jiang, S., Zhang, Z., Gao, W.: Video2cartoon: generating 3d cartoon from broadcast soccer video. In: MULTIMEDIA ’05 (2005)
  * Ngo and Cai [2008] Ngo, V., Cai, J.: Converting 2d soccer video to 3d cartoon. In: 2008 10th International Conference on Control, Automation, Robotics and Vision, pp. 103–107 (2008). https://doi.org/10.1109/ICARCV.2008.4795500
  * Weng et al. [2019] Weng, C.-Y., Curless, B., Kemelmacher-Shlizerman, I.: Photo wake-up: 3d character animation from a single photo. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) (2019)
  * Tous [2023] Tous, R.: Pictonaut: movie cartoonization using 3d human pose estimation and gans. Multimedia Tools and Applications, 1–15 (2023) https://doi.org/10.1007/s11042-023-14556-1
  * Ruben Tous and Igual [2023] Ruben Tous, J.N., Igual, L.: Human pose completion in partial body camera shots. Journal of Experimental & Theoretical Artificial Intelligence 0(0), 1–11 (2023) https://doi.org/10.1080/0952813X.2023.2241575
  * Fišer et al. [2017] Fišer, J., Jamriška, O., Simons, D., Shechtman, E., Lu, J., Asente, P., Lukáč, M., Sýkora, D.: Example-based synthesis of stylized facial animations. ACM Transactions on Graphics 36(4) (2017)
  * Gatys et al. [2015] Gatys, L.A., Ecker, A.S., Bethge, M.: A neural algorithm of artistic style. CoRR abs/1508.06576 (2015) 1508.06576
  * He et al. [2022] He, K., Chen, X., Xie, S., Li, Y., Dollar, P., Girshick, R.: Masked autoencoders are scalable vision learners. In: 2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pp. 15979–15988 (2022). https://doi.org/10.1109/CVPR52688.2022.01553
  * Dosovitskiy et al. [2021] Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., Uszkoreit, J., Houlsby, N.: An image is worth 16x16 words: Transformers for image recognition at scale. In: 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021. OpenReview.net, ??? (2021)
  * Radford et al. [2021] Radford, A., Kim, J.W., Hallacy, C., Ramesh, A., Goh, G., Agarwal, S., Sastry, G., Askell, A., Mishkin, P., Clark, J., Krueger, G., Sutskever, I.: Learning transferable visual models from natural language supervision. In: Meila, M., Zhang, T. (eds.) Proceedings of the 38th International Conference on Machine Learning. Proceedings of Machine Learning Research, vol. 139, pp. 8748–8763. PMLR, ??? (2021)
  * Carion et al. [2020] Carion, N., Massa, F., Synnaeve, G., Usunier, N., Kirillov, A., Zagoruyko, S.: End-to-end object detection with transformers. In: Computer Vision – ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part I, pp. 213–229. Springer, Berlin, Heidelberg (2020). https://doi.org/10.1007/978-3-030-58452-8_13 . https://doi.org/10.1007/978-3-030-58452-8_13
  * Yang and Yang [2022] Yang, Z., Yang, Y.: Decoupling features in hierarchical propagation for video object segmentation. In: Advances in Neural Information Processing Systems (NeurIPS) (2022)
  * Liu et al. [2023] Liu, S., Zeng, Z., Ren, T., Li, F., Zhang, H., Yang, J., Li, C., Yang, J., Su, H., Zhu, J., et al.: Grounding dino: Marrying dino with grounded pre-training for open-set object detection. arXiv preprint arXiv:2303.05499 (2023)
  * Suzuki and be [1985] Suzuki, S., be, K.: Topological structural analysis of digitized binary images by border following. Computer Vision, Graphics, and Image Processing 30(1), 32–46 (1985) https://doi.org/10.1016/0734-189X(85)90016-7
  * Kazemi and Sullivan [2014] Kazemi, V., Sullivan, J.: One millisecond face alignment with an ensemble of regression trees. In: 2014 IEEE Conference on Computer Vision and Pattern Recognition, pp. 1867–1874 (2014). https://doi.org/10.1109/CVPR.2014.241
  * Dalal and Triggs [2005] Dalal, N., Triggs, B.: Histograms of oriented gradients for human detection. In: 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), vol. 1, pp. 886–8931 (2005). https://doi.org/10.1109/CVPR.2005.177
  * Soomro et al. [2012] Soomro, K., Zamir, A.R., Shah, M.: UCF101: A dataset of 101 human actions classes from videos in the wild. CoRR abs/1212.0402 (2012) 1212.0402
  * Zablotskaia et al. [2019] Zablotskaia, P., Siarohin, A., Zhao, B., Sigal, L.: Dwnet: Dense warp-based network for pose-guided human video generation. In: 30th British Machine Vision Conference 2019, BMVC 2019, Cardiff, UK, September 9-12, 2019, p. 51. BMVA Press, ??? (2019)
  * Kay et al. [2017] Kay, W., Carreira, J., Simonyan, K., Zhang, B., Hillier, C., Vijayanarasimhan, S., Viola, F., Green, T., Back, T., Natsev, P., Suleyman, M., Zisserman, A.: The kinetics human action video dataset. CoRR abs/1705.06950 (2017) 1705.06950
  * [33] Tous, R.: Lester Project Home Page. Accessed: 2023-12-20. https://github.com/rtous/lester
  * ITU [2016] ITU: Mean opinion score (MOS) terminology. ITU-T Recommendation P.800.1 Methods for objective and subjective assessment of speech and video quality (2016). https://www.itu.int/rec/dologin_pub.asp?lang=s&id=T-REC-P.800.1-201607-I!!PDF-E&type=items
  * ITU [2008] ITU: Subjective Video Quality Assessment Methods for Multimedia Applications. ITU-T Recommendation P.910 (2008). https://www.itu.int/rec/T-REC-P.910


# Viscoelastic response of impact process on dense suspensions

Pradipto<EMAIL_ADDRESS>Hisao Hayakawa  Yukawa Institute for
Theoretical Physics, Kyoto University, Kitashirakawaoiwake-cho, Sakyo-ku,
Kyoto 606-8502, Japan

###### Abstract

We numerically study impact processes on dense suspensions using the lattice
Boltzmann method to elucidate the connection between the elastic rebound of an
impactor and relations among the impact speed $u_{0}$, maximum force acting on
the impactor $F_{\rm max}$, and elapsed time $t_{\rm max}$ to reach $F_{\rm
max}$. We find that $t_{\rm max}$ emerges in the early stage of the impact,
while the rebound process takes place in the late stage. We find a crossover
of $F_{\rm max}$ from $u_{0}$ independent regime for low $u_{0}$ to a power
law regime satisfying $F_{\rm max}\propto u_{0}^{\alpha}$ with $\alpha\approx
1.5$ for high $u_{0}$. Similarly, $t_{\rm max}$ satisfies $t_{\rm max}\propto
u_{0}^{\beta}$ with $\beta\approx-0.5$ for high $u_{0}$. Both power-law
relations for $F_{\rm max}$ and $t_{\rm max}$ versus $u_{0}$ for high $u_{0}$
are independent of the system size, but the rebound phenomenon strongly
depends on the depth of the container for suspensions. Thus, we indicate that
the rebound phenomenon is not directly related to the relations among $u_{0}$,
$F_{\rm max}$ and $t_{\rm max}$. We propose a floating + force chain model,
where the rebound process is caused by an elastic term that is proportional to
the number of the connected force chains from the impactor to the bottom
plate. On the other hand, there are no elastic contributions in the relations
for $F_{\rm max}$ and $t_{\rm max}$ against $u_{0}$ because of the absence of
percolated force chains in the early stage. This phenomenology predicts
$F_{\rm max}\propto u_{0}^{3/2}$ and $t_{\rm max}\propto u_{0}^{-1/2}$ for
high $u_{0}$ and also recovers the behavior of the impactor quantitatively
even if there is the rebound process.

††preprint: APS/123-QED

##  I Introduction

A dense suspension can behave as a liquid or a solid depending on the
situation. One of the most interesting behaviors of dense suspension is the
impact-induced hardening in which the suspension is solidified if the speed of
an impactor hitting on a suspension is high enough [1]. An example of this
non-Newtonian behaviors is a running person on the top of a cornstarch
suspension, while a walking person sinks [1]. This process is practically
important for various industrial applications such as protective vests [2]. A
similar hardening process is also observable in fractures on a thin layer of a
suspension under an impact [3]. Moreover, this kind of non-equilibrium solid-
liquid phase transition is interesting even for physicists. The impact-induced
hardening is often regarded as a process related to the discontinuous shear
thickening (DST) under simple shear [2, 4], which attracts much interest among
many researchers recently [5, 6, 7, 8, 9, 10]. Nevertheless, the underlying
mechanism of the impact-induced hardening differs from that of DST as
indicated by Ref. [11]. Indeed, the former is only dominated by the normal
stress, while both the normal and shear stresses play important roles in the
latter case. Thus, impact-induced hardening deserves to be studied on its own.

Let us review some previous studies on the impact-induced hardening.
Waitukaitis and Jaeger conducted an experiment with a rod impactor and
discovered the existence of a dynamically jammed region which is a solid plug
beneath the impactor [12]. They proposed the added-mass model to explain the
solidification induced by the impact. Then, a series of experiments found that
such solidifications take place when the dynamically jammed region is spanned
between the impactor and boundaries [4, 13, 14]. The impact-induced hardening
can be also observed by dropping an impactor into a dense suspension [15, 11].
As a result of the hardening, the free-falling impactor can rebound [15]. It
is obvious that the elastic effect of dense suspensions is responsible for
this rebound phenomenon.

Recently, some papers have discussed the relation between the impact speed
$u_{0}$ and the maximum force acting on the impactor $F_{\rm max}$ or the
elapsed time $t_{\rm max}$ to reach $F_{\rm max}$ in impact processes.
Previous studies[12, 14, 16] showed the existence of power-law relations such
as $F_{\rm max}\propto u_{0}^{\alpha}$ and $t_{\rm max}\propto u_{0}^{\beta}$.
It is noteworthy that similar relations are also found in impact processes for
dry granular materials [17]. The numerical solution of the added-mass model
[14] suggests $\alpha=2$ and $\beta=-1$, though the fitted values in their
experiment are $\alpha=1.5$ and $\beta=-1/2$. Moreover, a closer look at the
data in Ref. [12] suggested that $u_{0}$-independent exponents $\alpha$ and
$\beta$ are not appropriate to fit the data in all ranges of the impact speed.
A recent experiment [16] also suggested $\alpha=1.5$ and $\beta=-1/2$. These
values of the exponents are obtained as the solution of the viscous force
model [16], which is inspired by the existence of a growing dynamically jammed
region below the impactor [12, 18]. Nevertheless, the viscous force model [16]
has two defects in which (i) the model cannot explain the behavior for low
$u_{0}$ regime observed in Ref. [12], and (ii) the model cannot explain the
mechanism of the rebound process since any elastic term is absent. Therefore,
the connection between the rebound of the impactor and the relationships among
$u_{0}$, $F_{\rm max}$, and $t_{\rm max}$ should be clarified to understand
the viscoelastic response of an impactor on dense suspensions.

In this paper, we try to clarify the connection between the relations among
$u_{0}$, $F_{\rm max}$, and $t_{\rm max}$ and the rebound phenomena by
performing simulations of a free-falling impactor onto dense suspensions based
on a coupled model of the lattice Boltzmann method (LBM) and discrete element
method (DEM). In addition, we propose a phenomenology to explain these
processes, including the elastic force as a result of percolated force chains
between the impactor and bottom plate to describe the rebound phenomenon. This
model is essentially reduced to the viscous force model if percolated force
chains are absent.

The outline of the paper is as follows. In Sec. II, we briefly explain the
method and setup of our simulation. In Sec. III, we present the results of our
simulation including the impactor motion and the relationships among $u_{0}$,
$F_{\rm max}$ and $t_{\rm max}$. In Sec. IV we examine our phenomenology with
and without the elastic force between the impactor and bottom plate to explain
the results of the simulation. In Sec. V, we summarize our results and discuss
future perspectives. In Appendix A, we describe the details of our simulation
method. In Appendix B, we present the exact solution and its approximate
treatment of our phenomenology when the elastic force is absent. In Appendix
C, we describe the details of the force chains analysis in the phenomenology.
Finally in Appendix D, we discuss the dependence of our results on the volume
fraction of the suspensions.

##  II  Setup of our simulation

Figure 1: An illustration of an initial setup of our simulation. Figure 2: Successive snapshots of an impactor in a quasi-two-dimensional projection, where the black dashed lines correspond to the maximum penetration (deepest position) of the impactor. Volume fraction $\phi$ | Depth $H$ | Width $W$ | No. of particles $N$
---|---|---|---
0.00 | $3D_{I}$ | $6D_{I}$ | 0
0.10 | $3D_{I}$ | $6D_{I}$ | 409
0.25 | $3D_{I}$ | $6D_{I}$ | 1021
0.40 | $3D_{I}$ | $6D_{I}$ | 1634
0.48 | $3D_{I}$ | $6D_{I}$ | 1960
| $2D_{I}$ | $4D_{I}$ | 617
0.51 | $3D_{I}$ | $6D_{I}$ | 2083
| $6D_{I}$ | $5D_{I}$ | 2893
| $7D_{I}$ | $4D_{I}$ | 2160
| $2D_{I}$ | $4D_{I}$ | 642
0.53 | $2D_{I}$ | $6D_{I}$ | 1443
| $3D_{I}$ | $6D_{I}$ | 2164
| $7D_{I}$ | $4D_{I}$ | 2245
0.56 | $2D_{I}$ | $4D_{I}$ | 677
| $7D_{I}$ | $4D_{I}$ | 2371
Table 1: All variations of simulated volume fractions $\phi$ and box sizes
with the corresponding numbers of suspended particles $N$.

We simulate a binary mixture of suspension consisting of equal number of large
and small particles with bidispersity ratio
$a_{\text{max}}=1.2a_{\text{min}}$, where the radii of the large and small
particles are $a_{\text{max}}$ and $a_{\text{min}}$, respectively. We use the
mixture to avoid the crystallization of suspended particles in high density
regions. These suspended particles have the identical density $\rho_{p}$ which
is equal to the density $\rho_{f}$ of the solvent. The suspension liquid is
confined in a rectangular box surrounded by smooth sidewalls and a smooth
bottom plate. Since we simulate free falling processes of an impactor, there
is no lid above the container. The volume $V$ of the suspension liquid at rest
is expressed as $V=W\times D\times H$, where $H$ is the depth of the
suspension and $W=D$ is the width of the container as shown in Fig. 1. The
volume fraction $\phi$ of the suspension at rest without the impactor is
defined as $\phi=2N\pi(a_{\text{max}}^{3}+a_{\text{min}}^{3})/3V$, where $N$
is the number of suspended particles used in the simulation. The hydrodynamic
interaction among particles is simulated using the LBM. The contact force
between suspended particles is modeled by the DEM [19] with the spring
constant $k_{n}$ between contacting particles. We also introduce the
frictional contact model between suspended particles with Coulomb’s friction
rule, which is important to recover the hardening behavior of dense
suspensions [6, 7, 20, 11]. In this paper, we adopt the friction coefficient
$\mu=1$ for all cases. As it is known [8, 11], the rheological properties of
dense suspensions is insensitive to $\mu$ for $\mu\geq 0.3$ . Details of our
simulation method can be seen in Appendix A.

A spherical impactor with diameter $D_{I}$ and density $\rho_{I}$, is released
from the height $H_{0}$ which corresponds to the impact speed
$u_{0}=\sqrt{2gH_{0}}$ with the gravitational acceleration $g$. In our
simulation $\rho_{I}$ and $D_{I}$ satisfy $\rho_{I}=4\rho_{f}$ and
$D_{I}=6a_{\rm min}$, respectively. We also introduce the time scale
$t_{g}=\sqrt{a_{\text{min}}/2g}$, speed scale $u^{*}=\sqrt{2ga_{\text{min}}}$,
and force scale $F_{g}=\frac{4}{3}\pi\rho_{f}(D_{I}/2)^{3}g$. It should be
noted that there is another important time scale $t_{k}=\sqrt{m_{0}/k_{n}}$
where $m_{0}=\frac{4}{3}\pi\rho_{f}a_{\rm min}^{3}$. Thus, the mass of the
impactor $m_{I}$ is expressed as $m_{I}=\frac{\pi}{6}\rho_{I}D_{I}^{3}$. In
our simulation the ratio $t_{k}/t_{g}=0.045$ is fixed. Note that we also adopt
the DEM for the contact interactions between the impactor and suspended
particles, between suspended particles and container’s walls, and between the
impactor and container’s walls. All variations of simulated volume fractions
$\phi$ and box sizes are summarized in Table. 1. Note that we use three
ensembles for $\phi=0.53$, $W=D=6D_{I}$, and $H=3D_{I}$ and for $\phi=0.53$,
$W=D=6D_{I}$, and $H=2D_{I}$. We only simulate one ensemble for the other
cases.

##  III  Simulation results

###  III.1  Impactor motion

Figure 3: Plots of the time evolutions of the impactor motion (blue solid
lines) for $\phi=0.53,$ $W=D=6D_{I}$, and $H=2D_{I}$ for (a) the velocity
$u^{I}_{z}/u^{*}$ (black dashed line represents $u^{I}_{z}/u^{*}=0$), (b) the
position of the deepest point of the impactor $z(t)/a_{\rm min}$, (c) the
force exerted on the impactor $F^{I}_{z}/F_{g}$, and (d) the force exerted on
the bottom plate $F^{w,b}_{z}/F_{g}$, respectively. Dashed purple lines in
(a), (b), and (c) represent the solution of Eq. (6) and dashed light blue
lines in (a), (b), and (c) represent the solution of Eq. (4). Black dotted
lines highlight $F_{\rm max}$ and $t_{\rm max}$. Figure 4: Plots of the time
evolutions of the impactor velocities $u^{I}_{z}/u^{*}$ for $\phi=0.53$ and
$W=D=6D_{I}$ (left vertical axis) for various $u_{0}$ and $H$. Blue and cyan
solid lines represent the results for $H=2D_{I}$ with $u_{0}/u^{*}=3.2$, and
$H=2D_{I}$ with $u_{0}/u^{*}=2.6$, respectively. Red and yellow solid lines
represent the results for $H=3D_{I}$ with $u_{0}/u^{*}=3.2$, and $H=2D_{I}$
with $u_{0}/u^{*}=2.6$, respectively. Black dashed line represents
$u^{I}_{z}/u^{*}=0$. The dot-dashed lines represent the corresponding forces
exerted on the bottom plate $F^{w,b}_{z}/F_{g}$ (right vertical axis).

((a))

((b))

Figure 5: (a) Plots of maximum forces exerted on the impactor $F_{\rm max}$
scaled by the gravitational force $F_{g}$ against $u_{0}/u^{*}$ for various
container sizes, where the green dashed line represents
$1.58(u_{0}/u^{*})^{1.432}$. (b) Plots of time $t_{\rm max}$ to reach $F_{\rm
max}$ scaled by $t_{g}$ against $u_{0}/u^{*}$, where the green dashed line
represents $0.03(u_{0}/u^{*})^{-0.523}$. The blue solid lines in both figures
represent the solution of the floating model (Eq. (4)).

Figure 2 shows successive snapshots of an impactor in a quasi-two-dimensional
projection of our three-dimensional simulation, where the black dashed lines
correspond to the maximum penetration position of the impactor. We can also
identify a rebound process from the middle and right figures in which the
vertical position of the impactor in the right figure is higher than that in
the middle.

In Fig. 3, we plot the time evolutions of the velocity and vertical position
of the impactor, force acting on the impactor, and force acting on the bottom
plate obtained from our simulation. From Fig.3(a), one can see the existence
of a rebound process i.e. the region for $u_{z}^{I}<0$. As a result, the
vertical position $z(t)$ increases with time in the rebound process as shown
in Fig. 3(b). From Fig. 3(c), one can define $F_{\rm max}$ as its peak value
and $t_{\rm max}$ as the time to reach $F_{\rm max}$. Note that $t_{\rm max}$
coincides with the onset time of the force exerted on the bottom plate
$F^{w,b}_{z}$, while the rebound takes place around and after the peak of
$F^{w,b}_{z}$ (see Fig. 3(d)). This indicates that the rebound of the impactor
takes place when the force from the impactor is transmitted through the force
chains to the bottom plate. In Sec. IV.3, we will confirm this picture.

We also find that $t_{\rm max}$ is located much earlier than time of the
rebound region ($u_{z}^{I}<0$). This suggests that $F_{\rm max}$ and $t_{\rm
max}$ are not related to the rebound motion. This is consistent with the
following observation in which $F_{\rm max}$ and $t_{\rm max}$ are independent
of system size [16] but the rebound motion strongly depends on the system
size. Indeed, as can be seen in Fig. 4, the rebound takes place only for the
suspension in a shallow vessel as in the case of $H=2D_{I}$, while the rebound
cannot be observed for the suspension in a deep vessel ($H=3D_{I}$). At the
early stage for $t/t_{g}<0.05$, such depth dependence does not exist.

###  III.2  Relations among $u_{0}$, $F_{\rm max}$ and $t_{\rm max}$

In Fig. 5, we plot $F_{\rm max}$ exerted on the impactor scaled by the
gravitational force $F_{g}$ against $u_{0}$ for $\phi\geq 0.48$ (see Appendix
D for the results of $\phi\leq 0.40$). Here, the results of our simulation for
$F_{\rm max}$ and $t_{\rm max}$ show the existence of power-law regimes
satisfying

$F_{\rm max}\propto u_{0}^{\alpha},\quad t_{\rm max}\propto u_{0}^{\beta}$ (1)

with $\alpha=1.432\pm 0.0003$ and $\beta=-0.523\pm 0.042$ for $u_{0}>u^{*}$.
One can find that the data for all volume fractions and system sizes are
collapsed on a universal curve for $F_{\rm max}$, while $t_{\rm max}$ does not
have the beautiful data collapse. Our observed exponents agree with those in
the experiment [12, 16] and is smaller than the solution of the added-mass
model [14]. The values of $\alpha$ and $\beta$ also are close to those
obtained by the viscous force model [16]. This is understandable since the
peak of the force exists in the early stage where the elastic force to produce
the rebound does not play any role. Thus, one does not need to take into
account the elastic force to explain the relations among $u_{0}$, $F_{\rm
max}$, and $t_{\rm max}$. Moreover, we have simulated variations of widths and
depths in Fig. 5 to confirm that the relations among $u_{0}$, $F_{\rm max}$,
and $t_{\rm max}$ are independent of the system size. This is in contrast to
the rebound phenomenon which strongly depends on the width and depth of the
simulation box (see Fig. 4). This observation is another evidence that the
relations among $u_{0}$, $F_{\rm max}$, and $t_{\rm max}$ are not related to
the rebound phenomenon.

Our simulation also illustrates that a single power-law is insufficient for
$F_{\rm max}$ versus $u_{0}$ to fit the data in all ranges of the impact
speed. Instead, we find a crossover of the relation between $F_{\rm max}$ and
$u_{0}$ from $u_{0}$ independent regime for low $u_{0}$ to the power-law
region for high $u_{0}$ regime (see Fig. 5(a) ). The corresponding $u_{0}$
independent regime of $t_{\rm max}$ for low $u_{0}$ is also visible in Fig.
5(b), though the data are not clear enough. Even though the authors of Ref.
[12] did not mention such a crossover in their paper, their data suggest the
existence of a subtle crossover in the relation between $u_{0}$ and $F_{\rm
max}$, similar to what we have observed. Furthermore, it is obvious that a set
of single values of $\alpha$ and $\beta$ is no longer valid if the
acceleration due to gravity plays some roles, as will be shown in the next
section. This might be the reason why the viscous force model in Ref. [16]
cannot explain the existence of $u_{0}$ independent regime.

##  IV  PHENOMENOLOGY

###  IV.1 Overview

Judging from the observations in our simulation, we propose the following
simple phenomenology to describe the vertical motion of an impactor:

$m_{I}\frac{d^{2}z_{I}}{dt^{2}}=-m_{I}\tilde{g}+F^{I}_{D},$ (2)

where $z_{I}(t)$ is the vertical position of the center of mass of the
impactor, $\tilde{g}$ is the effective gravitational acceleration defined as
$\tilde{g}=g(\rho_{I}-\rho_{f})/\rho_{f}$, and $F^{I}_{D}$ is the drag force
acting on the impactor. It should be noted that Ref. [11] adopted the
dynamical Hertzian contact model (DHCM) but the predictions of the DHCM,
$6/5<\alpha<4/3$ and $-1/3<\beta<-1/5$, disagree with the simulation and
experimental results. DHCM has also another drawback in which it cannot
recover the $u_{0}$ independent regime observed in our simulation.

Recently, Brassard et al. [16] proposed the viscous force model including a
drag term which is proportional to the depth of the impactor, though their
model ignores the gravity term $m_{I}\tilde{g}$ and the elastic force to
reproduce the rebound process. Although their model cannot explain $u_{0}$
independent regime and the rebound process, the analytic solution of the model
yields $\alpha=1.5$ and $\beta=-0.5$. Our proposed model in the early stage is
essentially the same as that in Ref. [16] with keeping the gravity term (Sec.
IV.2). Of course, we should take into account the elastic force in the late
stage if there are percolated force chains from the impactor to the bottom
plate (Sec. IV.3). The model used in Sec. IV.3 reduces to the model in Sec.
IV.2 because the former contains the number of percolated force chains $n(t)$
which becomes zero in the latter case. Nevertheless, we will explain the
floating model with $n(t)=0$ in the early stage in Sec. IV.2, and introduce
the floating + force chain model with $n(t)\neq 0$ in the late stage in Sec.
IV.3 separately.

###  IV.2 Floating model

Figure 6: Plot of the drag exerted on the impactor $F^{I}_{D,z}$ scaled by
$F_{g}$ and impactor velocity $u^{I}_{z}/u^{*}$ against impactor depth scaled
by diameter of the impactor $z_{I}/D_{I}$. Black dashed line represents the
linear fit of $F^{I}_{D,z}/F_{g}$ for $z_{I}/D_{I}\leq 0.7$

Let us propose a simple phenomenology which we call the floating model to
explain the behavior of the impactor for both $F_{\rm max}$ and $t_{\rm max}$
in the early state. To model the motion of the impactor, we assume that the
impactor is only influenced by the gravity and viscous drag force from the
surrounding suspension in the early stage. This assumption is based on the
observation that the dynamically jammed region is floating without touching
the bottom plate in the early stage of the impact [12, 18]. We also assume
that the drag force is proportional to the impactor velocity because the fluid
drag should be determined by the Stokes flow. Thus, in order to extract the
coefficient, we plot the drag exerted on the impactor $F_{D}^{I}$ divided by
the velocity against $|z|$ in Fig. 6, where $z$ is the deepest position of the
impactor ($z=0$ is the instance of attachment of the impactor on the surface
of the liquid). Here, we confirm that the drag is proportional to the impactor
depth when $|z|/D_{I}\leq 0.7$. If we assume that $F_{D}^{I}$ is proportional
to $|z|\dot{z}_{I}$, $F_{D}^{I}$ in Eq. (2) is identical to the drag force in
the viscous force model [16]. It should be noted that the center of mass of
the impactor $z_{I}$ is related to $z$ as $z_{I}=z+a_{I}$. Then, the total
drag force acting on the impactor is given by

$F^{I}_{D}=3\pi\eta_{\rm eff}\dot{z}_{I}|z|,$ (3)

where we have introduced the effective viscosity $\eta_{\rm eff}$ to
characterize the apparent viscosity of the dynamically jammed region. The
derivation of Eq. (3) can be seen in Appendix B. Then, we can write the
equation of motion of the impactor as

$m_{I}\frac{d^{2}z_{I}}{dt^{2}}=-m_{I}\tilde{g}+3\pi\eta_{\rm
eff}\dot{z}_{I}|z|.$ (4)

Equation (4) can be solved exactly with the aid of the Airy functions (see
Appendix B). The solutions for $z_{I}(t)$, $u_{z}^{I}(t)=-dz_{I}(t)/dt$, and
$F_{z}^{I}(t)=m_{I}d^{2}z_{I}(t)/dt^{2}$ are plotted alongside the simulation
results in Fig. 3. The numerical solutions for $F_{\rm max}$ and $t_{\rm max}$
(the blue solid lines) are presented in Figs. 5(a) and 5(b), respectively. We
use the value of the effective viscosity $\eta_{\rm eff}=4.9\times
10^{4}m_{0}/(a_{\rm min}t_{g})$ as a fitting parameter. This value is about a
hundred times larger than the viscosity of the solvent $\eta_{0}$ and five
times larger than the observed viscosity for DST under simple shear using LBM
simulation [20] (see Appendix D). The enhancement of the viscosity is
reasonable, because the impactor contacts with the dynamically jammed region
which must have larger viscosity than that of the averaged suspension. One can
see that Eq. (4) can recover the crossover from $u_{0}$ independent regime for
low $u_{0}$ to the power law regime for high $u_{0}$ observed in our
simulations. The solution of Eq. (4) yields

$\alpha=\frac{3}{2},\quad\beta=-\frac{1}{2}$ (5)

for high $u_{0}$ (see Appendix B for how to obtain these exponents). $F_{\rm
max}$ should be independent of $u_{0}$ for low $u_{0}$ because the second term
on the right-hand side (r.h.s.) of Eq. (4) is much smaller than the first term
for low $u_{0}$. This is the simple explanation for the crossover observed in
Fig. 5.

###  IV.3 Floating + force chains model

Figure 7: Plots of time evolutions of velocities of the impactors for a
rebound case with $H=2D_{I}$ and a no-rebound case with $H=3D_{I}$. Filled
triangles represent the simulation results of $\phi=0.53,$ $W=D=6D_{I}$, and
$H=2D_{I}$. The green solid line represents the solution of Eq. (4), and the
purple solid line represents the solution of Eq. (6) (the black dashed line
represents $u^{I}_{z}/u^{*}=0$). Here, we also plot the simulation results for
$\phi=0.53$, $W=D=6D_{I}$, and $H=3D_{I}$ (black squares), where rebound does
not take place.

Unfortunately, Eq. (4) cannot explain the rebound of the impactor because of
the absence of elastic force which is the origin of the rebound (see Fig. 3).
This indicates the drawback of the viscous force model which cannot explain
the rebound process. We also note that the solution of Eq. (4) is independent
of system size, which is consistent with the results of Ref. [16] and Fig. 5
in this paper. On the other hand, our simulation in Sec. III.1 indicates that
rebound is related to the force acting on the bottom plate. Since the force
can be transmitted via contacts of suspended particles along the chains, we
can calculate the elastic force along the chains (see Appendix C for
definition and visualizations of the force chains). Then, we include an
elastic term to Eq. (4) caused by connected force chains between the impactor
and bottom plate as

$m_{I}\frac{d^{2}z_{I}}{dt^{2}}=-m_{I}\tilde{g}+3\pi\eta_{\rm
eff}\dot{z}_{I}|z|+n(t)k_{n}z_{I},$ (6)

where $n(t)$ is the number of connected chains from the impactor to the bottom
plate, and $k_{n}$ is the spring constant of the DEM. In other words, the
elastic force (the third term on the r.h.s. of Eq. (6)) is originated from the
contacting elastic force along the force chains of contacting suspended
particles between the impactor and bottom plate. Details of the algorithm to
determine $n(t)$ is written in Appendix C and is illustrated in the Fig. 12
(Multimedia view). It is obvious that Eq. (6) is reduced to Eq. (4) if the
percolated force chains do not exist, i. e. $n(t)=0$ in the early stage. In
this sense, the model in Eq. (6) is more general than the floating model
described by Eq. (4).

In Fig. 7, we plot time evolutions of the impactor velocity from our
simulation alongside with the corresponding results of Eqs. (4) and (6) with
$\eta_{\rm eff}=4.9\times 10^{4}m_{0}/(a_{\rm min}t_{g})$ and $k_{n}=2.5\times
10^{4}m_{0}/(a_{\rm min}t_{g}^{2})$ which is identical to that used in the DEM
simulation. Here, one can see that the rebound of the impactor can be
recovered by the introduction of the third term on the r.h.s. of Eq. (6) for
the shallow vessel case ($H=2D_{I}$) (see Fig. 7). On the other hand, the
floating model (Eq. (4)) is sufficient to recover the impactor velocity
correctly for the deep vessel case ($H=3D_{I}$) where rebound does not take
place. Thus, the phenomenology described by Eq. (6) can describe the
quantitative behavior of the impactor by the introduction of two fitting
parameters $\eta_{\rm eff}$ and $n(t)$, though $n(t)$ is determined by the
observation as shown in Appendix C. Thus, our phenomenology is more accurate
than the linear model in Ref. [15] and the DHCM in Ref. [11]. Note that the
power-law exponents $\alpha$ and $\beta$ in (5) are not affected by this
elastic term since $F_{\rm max}$ and $t_{\rm max}$ emerge in the early stage
of the impact.

##  V  Discussion and Conclusions

We numerically studied the impact processes on dense suspensions using a
coupled model of LBM and DEM to elucidate the connection between the elastic
rebound of the impactor and the relations among $u_{0}$, $F_{\rm max}$, and
$t_{\rm max}$. Then, we have also proposed a simple phenomenology called the
floating+force chain model to explain our simulation results. This model
reduces to the floating model if there are no percolated force chains from the
impactor to the bottom plate. We numerically find the existence of a power-law
regime satisfying $F_{\rm max}\propto u_{0}^{\alpha}$, with $\alpha=1.432\pm
0.0003$ and $t_{\rm max}\propto u_{0}^{\beta}$, with $\beta=-0.523\pm 0.042$,
while the analytic solution of the floating model indicates $\alpha=3/2$ and
$\beta=-1/2$. We have also confirmed the existence of $u_{0}$-independent
regimes of $F_{\rm max}$ and $t_{\rm max}$ for low $u_{0}$. The crossovers of
$F_{\rm max}$ and $t_{\rm max}$ from $u_{0}$-independent regimes to the power
law regimes can be reproduced by the floating model correctly. We conclude
that the relations among $u_{0}$, $F_{\rm max}$, and $t_{\rm max}$ are not
related to the rebound process based on three observations: (i) We found that
$F_{\rm max}$ emerges in the early stage of the impact, while the rebound of
the impactor takes place in the later stage. (ii) We have confirmed that the
relations among $u_{0}$, $F_{\rm max}$, and $t_{\rm max}$ are independent of
the system size, while the rebound strongly depends on the size of the
container. (iii) One can recover the exponents for $F_{\rm max}$ and $t_{\rm
max}$ when the impactor depth is smaller than its diameter (not completely
sink) in the drag term that is proportional to the impactor depth without
considering any elastic force, which agrees with Ref. [16]. In contrast, the
rebound needs an elastic term caused by the connected force chains from the
impactor to the bottom plate.

Our phenomenology, the floating + force chain model, is only valid for dense
suspensions because any percolated force chains do not exist for dilute
suspensions. Nevertheless, $F_{\rm max}$ and $t_{\rm max}$ can exist even for
dilute suspensions. Although the velocity of the impactor can be partially
explained by a model in which $F_{D}^{I}$ is expressed as the Stokes drag
force, the Stokes drag model cannot explain the existence of $F_{\rm max}$ and
$t_{\rm max}$ in dilute situations (see Appendix D). To estimate $\eta_{\rm
eff}$ within our phenomenology we need to take into account the interaction
between the impactor and the dynamically jammed region. However, the
dynamically jammed region is still not well-defined, and thus, this may not be
a well-defined problem. Furthermore, our algorithm to determine the number of
percolating force chains from the impactor to the bottom plate ($n(t)$)
essentially ignores the role of the sidewalls. This needs to be improved if
one examines the impact-induced hardening phenomena in a narrow channel or
using bumpy sidewalls. We only focused on relatively short time behavior after
the impact, while a sinking impactor in dense suspensions shows a distinct
behavior, as it oscillates and exhibits a stop-go cycle near the bottom of the
container [21]. Our simulation will be able to be used to reproduce these
results.

###### Acknowledgements.

One of the authors (P) expresses his gratitude to Alessandro Leonardi for
sharing his lattice Boltzmann code. The authors thank Satoshi Takada for his
critical reading of the manuscript. All numerical calculations were carried
out at the Yukawa Institute for Theoretical Physics (YITP) Computer
Facilities, Kyoto University, Japan. This research is partially supported by
Grants-in-Aid of MEXT, Japan for Scientific Research, Grant Nos. JP16H04025
and JP21H01006.

##  Appendix A  LBM + DEM with free surface

Figure 8: Summary of the LBM + DEM simulations. Figure 9: An illustration of
the bounce-rule in LBM simulation. Filled circles represent fluid nodes, open
circles represent solid nodes, open squares represent boundary nodes, arrows
represent the streaming discrete distribution functions. Figure 10: An
illustration of an impactor in the suspension liquid to explain $z$, $z_{I}$,
and $\theta_{0}$ (the black solid line connects the deepest position of the
impactor at $z$ with the center of mass at $z_{I}$, the red solid line is the
line between the surface of the suspension and the center of mass, the green
dashed line represents the deepest point $z$, the red dashed line is the
surface of the suspension ($z=0$), and $\theta_{0}$ is the angle between the
black and red solid lines.

We employ the LBM involving suspensions and the free surface of the fluid. The
summary of this method can be seen in Fig. 8. Throughout this paper, we have
adopted the perfect density matching between the solvent and suspended
particles, where the densities of particles and solvent satisfy the relation
$\rho_{p}=\rho_{f}$, where $\rho_{p}$ and $\rho_{f}$ are the densities of a
suspended particle and solvent fluid, respectively. The details of the LBM are
explained in Ref. [11]. The suspended particles in LBM are represented as a
group of solid nodes, while the surrounding fluids are represented by fluid
nodes. The hydrodynamic field is calculated from the time evolution of the
discrete distribution function at each fluid node. We select the lattice unit
$\Delta x=0.2a_{\text{min}}$, where it gives sufficient accuracy but still not
computationally expensive as shown in the previous LBM for suspensions
literatures [22, 23, 24]. In addition, to simulate the free surface of the
fluid, it is necessary to introduce interface nodes between the fluid and gas
nodes [25, 26, 27, 11].

Equations of motion and the torque balance of particle $i$ are, respectively,
given by

$m_{i}\frac{d\bm{u}_{i}}{dt}=\bm{F}_{i}^{c}+\bm{F}_{i}^{h}+\bm{F}_{i}^{\rm
lub}+\bm{F}_{i}^{r}+\bm{F}_{i}^{g},$ (7)
$I_{i}\frac{d\bm{\omega}_{i}}{dt}=\bm{T}_{i}^{c}+\bm{T}_{i}^{\rm
lub}+\bm{T}_{i}^{h}.$ (8)

Here, $\bm{u}_{i}$, $\bm{\omega}_{i}$, $m_{i}$, and
$I_{i}=(2/5)m_{i}a_{i}^{2}$ (with $a_{i}$ the radius of particle $i$), are the
translational velocity, angular velocity, mass, and the moment of inertia of
particle $i$, respectively. $\bm{F}_{i}^{g}=-m_{i}g\hat{\bm{z}}$ is the
gravitational force acting on the suspended particle $i$, where $\hat{\bm{z}}$
is the unit vector in the vertical direction.

Note that our LBM accounts for both the short-range lubrication force
$\bm{F}_{i}^{\rm lub}$ and torque $\bm{T}_{i}^{\rm lub}$, as well as the long-
range hydrodynamic force $\bm{F}_{i}^{h}$ and torque $\bm{T}_{i}^{h}$ as in
Ref. [24, 20]. The long-range parts ($\bm{F}_{i}^{h}$ and $\bm{T}_{i}^{h}$)
are calculated using the direct forcing method [27, 11], while the lubrication
force $\bm{F}_{i}^{\rm lub}$ and torque $\bm{T}_{i}^{\rm lub}$ are expressed
by pairwise interactions as $\bm{F}_{i}^{\rm lub}=\sum_{j\neq
i}\bm{F}_{ij}^{\rm lub}$ and $\bm{T}_{i}^{\rm c}=\sum_{j\neq
i}\bm{T}_{ij}^{\rm lub}$, respectively [5, 6, 24, 20]. The explicit
expressions of $\bm{F}_{ij}^{lub}$ and $\bm{T}_{ij}^{lub}$ can be found in
Ref. [20].

We adopt the linear spring-dashpot version of the DEM [19] for the contact
interaction between particles, which involves both the normal and the
tangential contact forces. Note that we omit the dissipative part for the
tangential contact force. For the particle $i$, the contact force
$\bm{F}^{c}_{i}$ and torque $\bm{T}^{c}_{i}$ are, respectively, written as
$\bm{F}^{c}_{i}=\sum_{i\neq
j}(\bm{F}^{\text{nor}}_{ij}+\bm{F}^{\text{tan}}_{ij})$ and
$\bm{T}^{c}_{i}=\sum_{i\neq j}a_{i}\bm{n}_{ij}\times\bm{F}^{\text{tan}}_{ij}$,
where $a_{i}$ is the radius of particle $i$. The normal force is explicitly
expressed as

$\bm{F}^{\text{nor}}_{ij}=(k_{n}\delta_{ij}^{n}-\zeta^{(n)}u_{ij}^{(n)})\bm{n}_{ij},$
(9)

where $k_{n}$ is the spring constant, $\delta_{ij}^{n}$ is the normal overlap,
$\bm{n}_{ij}$ is the normal unit vector between particles, $u_{ij}^{(n)}$ is
the normal velocity difference of the contact point
$u_{ij}^{(n)}=u_{i}^{(n)}-u_{j}^{(n)}$, and $\zeta^{(n)}=\sqrt{m_{0}k_{n}}$ is
the damping constant, where $m_{0}$ is the average mass of the suspended
particles. If the tangential contact force is smaller than a slip criterion,
tangential contact force is represented as

$\tilde{\bm{F}}^{\text{tan}}_{ij}=k_{t}\delta_{ij}^{t}\bm{t}_{ij},$ (10)

where $k_{t}$, assumed to be $0.2k_{n}$, is the tangential spring constant,
$\delta_{ij}^{t}$ is the tangential compression and $\bm{t}_{ij}$ is the
tangential unit vector at the contact point between particles $i$ and $j$. We
adopt the Coulomb friction rules as

$\displaystyle|\bm{F}_{ij}^{\text{tan}}|$
$\displaystyle=\mu|\bm{F}_{ij}^{\text{nor}}|\quad\text{if
}|\tilde{\bm{F}}^{\text{tan}}_{ij}|\geq\mu|\bm{F}^{\text{nor}}_{ij}|\quad\text{(slip)},$
(11) $\displaystyle|\bm{F}_{ij}^{\text{tan}}|$
$\displaystyle=|\tilde{\bm{F}}^{\text{tan}}_{ij}|\quad\text{if
}|\tilde{\bm{F}}^{\text{tan}}_{ij}|\leq\mu|\bm{F}^{\text{nor}}_{ij}|\quad\text{(stick)},$
(12)

whereas $\delta_{ij}^{t}$ is updated each time with relative tangential
velocity[19].

Finally, $\bm{F}_{i}^{r}$ is the electrostatic repulsive force, also expressed
by pairwise interactions as $\bm{F}_{i}^{\rm r}=\sum_{j\neq
i}\bm{F}_{ij}^{r}$. The explicit expression of $\bm{F}^{r}_{ij}$ is expressed
by the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory[28, 29, 30] for the
double layer electostatic force as

$\bm{F}^{r}_{ij}=F_{0}\exp(-h/\lambda)\bm{n}_{ij},$ (13)

where
$F_{0}=k_{B}T\lambda_{B}\hat{Z}^{2}(e^{a_{\text{min}}/\lambda}/(1+a_{\text{min}}/\lambda))^{2}/h^{2}$
with the charge number $\hat{Z}$, the Bjerrum length $\lambda_{B}$ and the
Debye-Hückel length $\lambda$. Note that $\lambda_{B}$ can be expressed as
$\lambda_{B}=e^{2}/(4\pi\epsilon_{0}\epsilon_{r}k_{B}T)$ where $e$,
$\epsilon_{0}$, $\epsilon_{r}$, and $k_{B}$ are the elementary charge, the
vacuum permittivity, the dielectric constant, and the Boltzmann constant,
respectively[30]. Here, we adopt the Debye length
$\lambda=0.02a_{\text{min}}$. Our simulation ignores the Brownian force. Thus,
the electrostatic repulsion force is important to prevent the suspended
particles from clustering [20, 6].

The impactor is a solid spherical object with the density
$\rho_{I}=4\rho_{f}$. The force and torque acting on the impactor are,
respectively, given by

$\bm{F}^{I}=\bm{F}^{I,h}+\bm{F}^{I,\rm lub}+\bm{F}^{I,c}+\bm{F}^{I,g},$ (14)
$\bm{T}^{I}=\bm{T}^{I,h}+\bm{T}^{I,c}+\bm{T}^{I,\rm lub}.$ (15)

$\bm{F}^{I,g}=-m_{I}g\bm{\hat{z}}$ is the gravitational force acting on the
impactor. The contact force $\bm{F}^{I,c}$ and torque $\bm{T}^{I,c}$, which
arise from the interactions with the suspended particles, are also calculated
by the DEM. The lubrication force $\bm{F}^{I,\rm lub}$ and torque
$\bm{T}^{I,\rm lub}$ are also calculated in a similar manner as used in
suspended particles.

The long-range hydrodynamic force $\bm{F}^{I,h}$ and torque $\bm{T}^{I,h}$ are
calculated using the bounce-back rule which satisfies the no-slip boundary
condition between the fluid and the surface of the impactor [22, 23]. In the
bounce-back rule the LBM discrete distribution function that streams from
fluid nodes to the boundary nodes is reflected. Then, the hydrodynamic force
on each node is calculated from the momentum transferred in this reflection
process. In our implementation, the bounce-back rule is implemented by
treating the surface of the impactor as boundary nodes. An illustration of
this bounce-back rule can be seen in Fig. 9.

Figure 11: An illustration of the terminology in force chains. Lines represent
the links, circles represent the nodes.

##  Appendix B  Derivation and analytical solution of the floating model

###  B.1 Derivation of $F_{D}^{I}$

((a))

((b))

((c))

((d))

Figure 12: Illustrations of the algorithm to determine $n(t)$ from force
chains network (Multimedia view). (a) An initial network pf force chains. (b)
The remaining chains after the lateral chains are removed. (c) All connected
components after removing all edges that do not touch the impactor or bottom
plate, where blue and red lines represent the corresponding connected
components that self-loops and chains between the impactor and bottom plate,
respectively. (d) Percolated force chains from the impactor to the bottom
plate.

The linear relationship between $|z|$ and the drag force in Fig. 6 may be
understood by the following simple model. For $|z|<D_{I}$, the surface of the
impactor is partially surrounded by the liquid and some parts of the surface
are still in the air. Here, we assume that Stokes drag law can be used for the
region surrounded by the liquid. Thus, Stokes’ drag force consists of two
parts, the pressure drag $F_{D,p}^{I}$ and friction drag $F_{D,f}^{I}$ as
$F^{I}_{D}=F^{I}_{D,p}+F^{I}_{D,f}$ [31],

$\displaystyle F^{I}_{D,p}$ $\displaystyle=3\pi\eta_{\rm
eff}a_{I}\dot{z}_{I}\int_{0}^{\theta_{0}}\cos^{2}\theta\sin\theta d\theta,$
(16) $\displaystyle F^{I}_{D,f}$ $\displaystyle=3\pi\eta_{\rm
eff}a_{I}\dot{z}_{I}\int_{0}^{\theta_{0}}\sin^{3}\theta d\theta$
$\displaystyle=3\pi\eta_{\rm
eff}a_{I}\dot{z}_{I}(1-\cos\theta_{0})-F^{I}_{D,p},$ (17)

where $a_{I}$ is the radius of the impactor satisfying $a_{I}=D_{I}/2$,
$\eta_{\rm eff}$ is the effective viscosity of the surrounding fluid, and
$\theta_{0}$ is the separation angle between moving direction ($\theta=0$) and
the line from the impactor center to the surface of the liquid (see Fig. 10).
Note that the deepest position of the impactor satisfies the relation
$|z|=a(1-\cos\theta_{0})$. Thus, one can reach Eq. (3).

###  B.2 Analytical solution of Eq. (4)

The dimensionless form of Eq. (4) with the aid of $a_{I}$ and
$t_{I}=\sqrt{a_{I}/\tilde{g}}$ is written as

$\ddot{z}^{*}_{I}=-1+\eta^{*}\dot{z}^{*}_{I}|z^{*}_{I}-1|,$ (18)

where $z^{*}_{I}=z_{I}/a_{I}$, $t^{*}=t/t_{I}$,
$\dot{z_{I}}^{*}=dz^{*}_{I}(t)/dt^{*}=-u_{z}^{I}/u^{*}$,
$\ddot{z_{I}}^{*}=d^{2}z^{*}_{I}(t)/dt^{*2}$, and $\eta^{*}=3\pi\eta_{\rm
eff}a_{I}\sqrt{a_{I}/\tilde{g}}/m_{I}$. Then, Eq. (18) can be solved exactly
in terms of the Airy functions as

$z^{*}_{I}(t)=\frac{\kappa[-\text{Ai}^{\prime}(\Phi)\text{Bi}^{\prime}(\Theta)+\text{Ai}^{\prime}(\Theta)\text{Bi}^{\prime}(\Phi)]}{\gamma[\text{Bi}(\Phi)\text{Ai}^{\prime}(\Theta)-\text{Ai}(\Phi)\text{Bi}^{\prime}(\Theta)]},$
(19)

where $\gamma=(\eta^{*})^{2/3}$, $\kappa=2^{2/3}$,
$\Theta=u^{*}_{0}\sqrt[3]{\eta^{*}/2}$, and
$\Phi=(u_{0}^{*}+t^{*})\sqrt[3]{\eta^{*}/2}$, where
$u_{0}^{*}=u_{0}t_{I}/a_{I}$. Here, $\text{Ai}(x)$ is the Airy function of the
first kind, which is defined as
$\text{Ai}(x)=\int_{0}^{\infty}\cos(t^{3}/3+xt)dt/\pi$, and
$\text{Ai}^{\prime}(x)$ is its derivative. $\text{Bi}(x)$ is the Airy function
of the second kind, which is defined as
$\text{Bi}(x)=\int_{0}^{\infty}[\exp(-t^{3}/3+xt)+\sin(-t^{3}/3+xt)]dt/\pi$,
and $\text{Bi}^{\prime}(x)$ is its derivative. One can differentiate Eq. (19)
two times to get the expression for $\ddot{z}_{I}$ as

$\displaystyle\ddot{z}^{*}_{I}=$
$\displaystyle\frac{\begin{matrix}\bigg{[}\kappa\gamma^{\frac{3}{2}}\bigg{(}u^{*}_{0}-t^{*}\bigg{)}\bigg{(}\text{Ai}^{\prime}(\Theta)\text{Bi}(\Phi)-\text{Ai}(\Phi)\text{Bi}^{\prime}(\Theta)\bigg{)}^{2}\bigg{(}\text{Ai}^{\prime}(\Theta)\text{Bi}^{\prime}(\Phi)-\text{Ai}^{\prime}(\Phi)\text{Bi}^{\prime}(\Theta)\bigg{)}-\gamma\Lambda\bigg{]}\end{matrix}}{\bigg{(}\text{Ai}^{\prime}(\Theta)\text{Bi}(\Phi)-\text{Ai}(\Phi)\text{Bi}^{\prime}(\Theta)\bigg{)}^{3}},$
$\displaystyle\Lambda=$
$\displaystyle\text{Ai}^{\prime}(\Theta)^{3}\text{Bi}(\Phi)^{3}+2\text{Ai}^{\prime}(\Theta)^{3}\text{Bi}^{\prime}(\Phi)^{3}-3\text{Ai}(\Phi)\text{Ai}^{\prime}(\Theta)^{2}\text{Bi}^{\prime}(\Phi)^{2}\text{Bi}^{\prime}(\Theta)-6\text{Ai}^{\prime}(\Phi)\text{Ai}^{\prime}(\Theta)^{2}\text{Bi}^{\prime}(\Phi)^{2}\text{Bi}^{\prime}(\Theta)$
$\displaystyle+3\text{Ai}(\Phi)^{2}\text{Ai}^{\prime}(\Theta)\text{Bi}(\Phi)\text{Bi}^{\prime}(\Theta)^{2}+6\text{Ai}^{\prime}(\Phi)^{2}\text{Ai}^{\prime}(\Theta)\text{Bi}^{\prime}(\Phi)\text{Bi}^{\prime}(\Theta)^{2}+\text{Ai}(\Phi)^{3}\text{Bi}^{\prime}(\Theta)^{3}-2\text{Ai}^{\prime}(\Phi)^{3}\text{Bi}^{\prime}(\Theta)^{3}.$
(20)

To obtain the expression of $F_{\rm max}$ and $t_{\rm max}$, we adopt the
short time expansion for Eq. (20) since $F_{\rm max}$ appears in the region
$t/t_{g}\ll 1$, Thus, up to third order, one can obtain

$\displaystyle\ddot{z}^{*}_{I}=1$
$\displaystyle-\eta^{*}u_{0}^{*2}t^{*}-\frac{3\eta^{*}u^{*}_{0}t^{*2}}{2}$
$\displaystyle+\bigg{(}\frac{2u^{*3}_{0}(\eta^{*})^{2}}{3}-\frac{\eta^{*}}{2}\bigg{)}t^{*3}+\mathcal{O}\bigg{(}\left[\frac{t^{*}}{u^{*}_{0}}\right]^{4}\bigg{)}$
(21)

Then, we differentiate Eq. (21) to obtain $\dddot{z}^{*}_{I}$ as

$\displaystyle\dddot{z}^{*}_{I}=$
$\displaystyle-\eta^{*}u_{0}^{*2}-3\eta^{*}u^{*}_{0}t^{*}$
$\displaystyle+\bigg{(}\frac{2u^{*3}_{0}(\eta^{*})^{2}}{3}-\frac{\eta^{*}}{2}\bigg{)}t^{*2}+\mathcal{O}\bigg{(}\left[\frac{t^{*}}{u^{*}_{0}}\right]^{3}\bigg{)}$
(22)

Then, for $\dddot{z}^{*}_{I}=0$, one can solve the quadratic equation in Eq.
(22) for $t_{\rm max}$ as

$\frac{t_{\rm
max}}{t_{I}}=\frac{3u^{*}_{0}+\sqrt{3u_{0}^{*2}+8\eta^{*}u_{0}^{*5}}}{4\eta^{*}u_{0}^{*3}-3}.$
(23)

For $\eta_{\rm eff}u_{0}\gg 1$, Eq. (23) reduces to

$\frac{t_{\rm max}}{t_{I}}=\frac{u_{0}^{*-\frac{1}{2}}}{\sqrt{2}}.$ (24)

Thus, we confirm the exponent $\beta=-1/2$. To obtain $F_{\rm max}$, we plug
Eq. (24) into Eq. (21) and take the limit $\eta_{\rm eff}u_{0}\gg 1$. Thus, we
obtain

$\frac{F_{\rm
max}}{m_{I}\tilde{g}}=u_{0}^{*\frac{3}{2}}\sqrt{\frac{2\eta^{*}}{9}}.$ (25)

Thus, we confirm the exponent $\alpha=3/2$ for large $u_{0}$. From Eq. (23),
$t_{\rm max}$ diverges at $u_{0,c}=\sqrt[3]{3/4\eta^{*}}$. This result
suggests the limitation of the short time approximation.

##  Appendix C  Determination of $n(t)$ in the floating + force chain model

In this Appendix, we explain the algorithm to determine the connected force
chains from the impactor to the bottom plate used in the floating + force
chain model in Sec. IV.3. First of all, let us explain how we draw the force
chains. Note that force chains are defined as a collection of nodes and links
representing the contacting suspended particles (see Fig. 11). Thus, for each
pair of contacting suspended particles, we draw a network in which a node
represents the center of a contacting pair of particles and a link is a
straight line connecting a pair of adjacent nodes. The initial force chains
can be seen in Fig. 12(a).

The algorithm to determine $n(t)$ is as follows. Since we are only interested
in the force propagation in the vertical direction, we remove links that
expand in the lateral directions (dangling chains). For this purpose, we
remove all links in which the height difference $|z_{i}-z_{j}|$ for a
contacting pair of particles $i$ and $j$ is less than the smallest radius of
the suspended particles $a_{\rm min}$. The corresponding network after the
removal of lateral chains can be seen in Fig. 12(b).

Our goal is to determine connected networks from the impactor to the bottom
plate. Thus, we remove all links to reach the edges of the force chains which
do not touch the bottom plate nor the impactor. Once such links are removed,
the leftover chains create new links at the edges of the force chain. We
repeat these labeling and removal processes until there are no edges of
dangling chains except for the edges which touch the bottom plate or the
impactor. Then, we label each connected component (the blue and red connected
components in Fig. 12(c)).

Note that the connected components do not need to be percolated from the
impactor to the bottom plate to survive in our algorithm at this stage due to
the existence of connected edges which form a self-loop (blue connected
components in Fig. 12(c)). Therefore, we need to examine whether each
connected component is percolated or not. Then, we remove non-percolated
connected components (blue component) while keeping the percolated connected
component (red component) as shown in Fig. 12(d). Finally, we evaluate $n(t)$
by the number of links that touch the bottom plate. The above processes are
illustrated in Fig. 12 (Multimedia view). The obtained $n(t)$ for $\phi=0.53,$
$W=D=6D_{I}$, and $H=2D_{I}$ with $u_{0}=2.6u^{*}$ against time is plotted in
Fig. 13.

Figure 13: A plot of the number of connected force chains from the impactor to
the bottom plate $n(t)$ against time for $\phi=0.53,$ $W=D=6D_{I}$, and
$H=2D_{I}$ with $u_{0}=2.6u^{*}$.

##  Appendix D  Dependence on volume fraction of the suspensions.

Figure 14: A phase diagram showing whether the impactor has $F_{\rm max}$ as a
function of the volume fraction $\phi$ and the impact speed $u_{0}$ for
$W=D=6D_{I}$ and $H=3D_{I}$. Red squares represent set of parameters where the
relations among $u_{0}$, $F_{\rm max}$, and $t_{\rm max}$ can be explained by
Eq. (4). Green squares are points where $F_{\rm max}$ exist but Eq. (4) fails.
Blue squares are where $F_{\rm max}$ does not even exist.

((a))

((b))

Figure 15: Plots of forces exerted on the impactor against time for various
volume fractions for $W=D=6D_{I}$ and $H=3D_{I}$ for (a) $u_{0}/u^{*}=5.84$
(Dashed lines represent the solutions of Eq. (26)) and (b) $u_{0}/u^{*}=0.93$.
Figure 16: Plots of the deepest points of the impactor scaled by the diameter
of the impactor for $u_{0}/u^{*}=5.84$ for $W=D=6D_{I}$ and $H=3D_{I}$.
Figure 17: Plots of the velocities of the impactor for low volume fractions
and for $W=D=6D_{I}$ and $H=3D_{I}$ alongside with the solutions of Eq. (26).

In this Appendix, we have examined whether the relations among $u_{0}$,
$F_{\rm max}$, and $t_{\rm max}$ only exist in the impact process in dense
suspensions, though the rebound only exists in dense suspensions. We summarize
the dependence on the volume fraction in the phase diagram in Fig. 14. Our
simulation indicates that $F_{\rm max}$ only exists in all range of $u_{0}$
when $\phi\geq 0.48$, while $F_{\rm max}$ does not exist for dilute
suspensions except for very high $u_{0}$. This result clarifies the role of
suspensions in which our analysis in the main text is only valid for dense
suspensions.

In Fig. 15, we plot the force exerted on the impactor against time to clarify
the difference between dense and dilute cases. For high $u_{0}$ (Fig. 15(a)),
$t_{\rm max}$ for the dilute case emerges earlier than that in the denser
case. In addition, dilute cases have smaller $F_{\rm max}$. Such differences
occur since the origin of $F_{\rm max}$ in a dilute case is different from
that in the dense case. In dense situations, the dominant contribution is from
the contact force between the impactor and the suspended particles, while for
the dilute situations, the dominant contribution is from the hydrodynamic
force exerted on the impactor [11]. For low $u_{0}$ (Fig. 15(b)), one can see
that $F_{\rm max}$ only exists in dense situation. Since the acceleration due
to the gravity is dominant for low $u_{0}$, the sufficient drag resistance to
compete with the gravity forces only exists for dense suspensions.

In Fig. 16, we plot the deepest point of the impactor $|z|$ (see Fig. 10)
scaled by the impactor diameter $D_{I}$ against time. Here, one can see that
the impactor sinks right after the impact in dilute suspensions, while the
impactor can keep its position near the surface for dense suspensions. The
behavior in which the impactor stays for a while near the surface of the
suspension is a characteristic of dense suspensions under impact. Thus, the
floating model (Eq. (4)) cannot be used for dilute situations since the
floating model assumes that the impactor is partially surrounded by fluid. We
also summarize the region where Eq. (4) is applicable in Fig. 14. When the
impactor is completely sink, the second term on the r.h.s. of Eq. (4) should
be replaced by the Stokesian drag as

$m_{I}\frac{d^{2}z_{I}}{dt^{2}}=-m_{I}\tilde{g}+3\pi\eta_{\rm
eff}a_{I}\dot{z}_{I}.$ (26)

In Fig. 17, we plot the impactor velocity against time alongside the solutions
of Eq. (26) for the dilute cases. Note that Eq. (4) cannot describe even the
behavior of the impactor velocity in dilute cases.

As expected, the apparent viscosity $\eta_{\rm eff}/\eta_{0}$, where
$\eta_{0}$ is the solvent viscosity, becomes larger as the volume fraction
increases. Although the agreement between the solution of Eq. (26) and our
simulation is remarkable for $u_{z}^{I}$ (see Fig. 17), Eq. (26) cannot
capture $F_{\rm max}$ in dilute suspensions (see the dashed lines in Fig.
15(a)). This is because there is no competition between time increasing and
time decreasing contributions in Eq. (26). Although one may extend the studies
on the impact process on water (without suspended particles) alone to dilute
suspensions [32], such a problem is beyond the scope of this paper.

## References

  * Brown and Jaeger [2014] E. Brown and H. M. Jaeger, “Shear thickening in concentrated suspensions: phenomenology, mechanisms and relations to jamming,” Reports on Progress in Physics 77, 046602 (2014).
  * Lee, Wetzel, and Wagner [2003] Y. S. Lee, E. D. Wetzel, and N. J. Wagner, “The ballistic impact characteristics of kevlar woven fabrics impregnated with a colloidal shear thickening fluid,” J. Mater. Sci. 38, 2825–2833 (2003).
  * Roche  _et al._ [2013] M. Roche, E. Myftiu, M. C. Johnston, P. Kim, and H. A. Stone, “Dynamic fracture of nonglassy suspensions,” Phys. Rev. Lett. 110, 148304 (2013).
  * Allen  _et al._ [2018] B. Allen, B. Sokol, S. Mukhopadhyay, R. Maharjan, and E. Brown, “System-spanning dynamically jammed region in response to impact of cornstarch and water suspensions,” Phys. Rev. E. 97, 052603 (2018).
  * Seto  _et al._ [2013] R. Seto, R. Mari, J. F. Morris, and M. M. Denn, “Discontinuous shear thickening of frictional hard-sphere suspensions,” Phys. Rev. Lett. 111, 218301 (2013).
  * Mari and Seto [2014] R. Mari and R. Seto, “Shear thickening, frictionless and frictional rheologies in non-brownian suspensions,” J. Rheol 58, 1693–1724 (2014).
  * Townsend and Wilson [2017] A. K. Townsend and H. J. Wilson, “Frictional shear thickening in suspensions: The effect of rigid asperities,” Phys. Fluids 29, 121607 (2017).
  * Sivadasan  _et al._ [2019] V. Sivadasan, E. Lorenz, A. G. Hoekstra, and D. Bonn, “Shear thickening of dense suspensions: The role of friction,” Phys. Fluids 31, 103103 (2019).
  * Gado and Morris [2020] E. D. Gado and J. F. Morris, “Preface: Physics of dense suspensions,” J. Rheol 64, 223 (2020).
  * Jamali, Gado, and Morris [2020] S. Jamali, E. D. Gado, and J. F. Morris, “Rheology discussions: Physics of dense suspensions,” J. Rheol 64, 1501 (2020).
  * Pradipto and Hayakawa [2021] Pradipto and H. Hayakawa, “Impact-induced hardening in dense frictional suspensions,” Phys. Rev. Fluids 6, 033301 (2021).
  * Waitukaitis and Jaeger [2012] S. R. Waitukaitis and H. M. Jaeger, “Impact-activated solidification of dense suspensions via dynamic jamming fronts,” Nature 487, 205–209 (2012).
  * Maharjan  _et al._ [2018] R. Maharjan, S. Mukhopadhyay, B. Allen, T. Storz, and E. Brown, “Constitutive relation for the system-spanning dynamically jammed region in response to impact of cornstarch and water suspensions,” Phys. Rev. E. 97, 052602 (2018).
  * Mukhopadhyay, Allen, and Brown [2018] S. Mukhopadhyay, B. Allen, and E. Brown, “Testing constitutive relations by running and walking on cornstarch and water suspensions,” Phys. Rev. E. 97, 052604 (2018).
  * Egawa and Katsuragi [2019] K. Egawa and H. Katsuragi, “Bouncing of a projectile impacting a dense potato-starch suspension layer,” Phys. Fluids 31, 053304 (2019).
  * Brassard  _et al._ [2021] M. Brassard, N. Causley, N. Krizou, J. A. Dijksman, and A. H. Clark, “Viscous-like forces control the impact response of shear-thickening dense suspensions,” J. Fluid. Mech. 923, A38 (2021).
  * Krizou and Clark [2020] N. Krizou and A. H. Clark, “Power-law scaling of early-stage forces during granular impact,” Phys. Rev. Lett. 124, 178002 (2020).
  * Han, Peters, and Jaeger [2016] E. Han, I. R. Peters, and H. M. Jaeger, “High-speed ultrasound imaging in dense suspensions reveals impact-activated solidification due to dynamic shear jamming,” Nat. Commun. 7, 12243 (2016).
  * Luding [2008] S. Luding, “Cohesive, frictional powders: contact models for tension,” Granul. Matter. 10, 235–246 (2008).
  * Pradipto and Hayakawa [2020] Pradipto and H. Hayakawa, “Simulation of dense non-brownian suspensions with the lattice boltzmann method: shear jammed and fragile states,” Soft Matter 16, 945–959 (2020).
  * von Kann  _et al._ [2011] S. von Kann, J. H. Snoeijer, D. Lohse, and D. van der Meer, “Nonmonotonic settling of a sphere in a cornstarch suspension,” Phys. Rev. E. 84, 060401(R) (2011).
  * Ladd [1994a] A. J. C. Ladd, “Numerical simulations of particulate suspensions via a discretized boltzmann equation. part 1. theoretical foundation,” J. Fluid Mech 271, 285–309 (1994a).
  * Ladd [1994b] A. J. C. Ladd, “Numerical simulations of particulate suspensions via a discretized boltzmann equation. part 2. numerical simulations,” J. Fluid Mech 271, 311–339 (1994b).
  * Nguyen and Ladd [2002] N. Q. Nguyen and A. J. C. Ladd, “Lubrication corrections for lattice-boltzmann simulations of particle suspensions,” Phys. Rev. E. 66, 046708 (2002).
  * Švec  _et al._ [2012] O. Švec, J. Skocek, H. Stang, M. R. Geiker, and N. Roussel, “Free surface flow of a suspension of rigid particles in a non-newtonian fluid: A lattice boltzmann approach,” J. Non-Newton. Fluid. 179-180, 32–42 (2012).
  * Leonardi  _et al._ [2014] A. Leonardi, F. K. Wittel, M. Mendoza, and H. J. Hermann, “Coupled dem-lbm method for the free-surface simulation of heterogeneous suspensions,” Comp. Part. Mech. 1, 3–13 (2014).
  * Leonardi  _et al._ [2015] A. Leonardi, M. Cabrera, F. K. Wittel, R. Kaitna, M. Mendoza, W. Wu, and H. J. Herrmann, “Granular-front formation in free-surface flow of concentrated suspensions,” Phys. Rev. E. 92, 052204 (2015).
  * Derjaguin and Landau [1941] B. Derjaguin and L. D. Landau, “Theory of the stability of strongly charged lyophobic sols and of the adhesion of strongly charged particles in solutions of electrolytes,” Acta Physicochim. U.R.S.S 14, 633–662 (1941).
  * Verwey and Overbeek [1948] E. J. W. Verwey and J. T. G. Overbeek,  _Theory of the Stability of Lyophobic Colloids: The Interaction of Sol Particles Having an Electric Double Layer_ (Elsevier, 1948).
  * Israelachvili [2011] J. Israelachvili,  _Intermolecular and Surface Forces_, 3rd ed. (Academic Press, 2011).
  * Batchelor [2000] G. K. Batchelor,  _An Introduction to Fluid Dynamics_ (Cambridge University Press, 2000).
  * Moghisi and Squire [1981] M. Moghisi and P. T. Squire, “An experimental investigation of the initial force of impact on a sphere striking a liquid surface,” J. Fluid. Mech. 108, 133 (1981).