††thanks: 1/A, Pázmány Péter Walkway, 1117 Budapest, HUNGARY

e-mail<EMAIL_ADDRESS>64/13, Volodymyrska Street, City of
Kyiv, Ukraine, 01601 e-mail<EMAIL_ADDRESS>
# Pomeron-pomeron scattering

István Szanyi  Eötvös Loránd University  Volodymyr Svintozelskyi  Taras
Shevchenko National University of Kyiv

###### Abstract

Central exclusive diffractive (CED) production of meson resonances potentially
is a factory producing new particles, in particular a glueball. The produced
resonances lie in trajectories with vacuum quantum numbers, essentially on the
pomeron trajectory. A tower of resonance recurrences, the production cross
section and the resonances widths are predicted. A new feature is the form of
the non-linear pomeron trajectory, producing resonances (glueballs) with
increasing widths. At LHC energies, in the nearly forward direction the
$t$-channel both in elastic, single or double diffraction dissociation as well
as in CED is dominated by pomeron exchange (the role of secondary trajectories
is negligible, however a small contribution from the odderon may be present).

Regge trajectory, pomeron, glueball, CED, LHC.

##  1 Introduction

Central exclusive diffractive (CED) production continues attracting attention
of both theorists and experimentalists, see e.g. Ref. [1] and references
therein. Interest in this subject is triggered by LHC’s high energies, where
even the sub energies at an equal partition is sufficient to neglect the
contribution from secondary Regge trajectories and consequently CED can be
considered as a gluon factory to produce exotic particles such as glueballs.

Fig. 1.: Regge-pole factorization.

Below we study CED shown in Fig. 1 with topology 4. Its knowledge is essential
in studies with diffractive excited protons, topologies 5 and 6.

In single-diffraction dissociation or single dissociation (SD) one of the
incoming protons dissociates (topology 2 in Fig. 1), in double-diffraction
dissociation or double dissociation (DD) both protons dissociate (topology 3)
and in central dissociation (CD) or double-Pomeron exchange (DPE) neither
proton dissociates (topology 4). These processes are tabulated below as,

$\displaystyle{\rm SD}$ $\displaystyle pp\rightarrow Xp$ $\displaystyle{\rm
or}$ $\displaystyle pp\rightarrow pY$ $\displaystyle{\rm DD}$ $\displaystyle
pp\rightarrow XY$ $\displaystyle{\rm CD}\;{\rm(DPE)}$ $\displaystyle
pp\rightarrow pZp,$

where $X$ and $Y$ represent diffractive dissociated protons while $Z$ denotes
a central system, consisting of meson/glueball resonances.

##  2 Pomeron/glueball trajectory

Regge trajectories $\alpha(s)$ connect the scattering region, $s<0$, with that
of particle spectroscopy, $s>0$. In this way they realize crossing symmetry
and anticipate duality i.e. dynamics of two kinematically disconnected regions
is intimately related: the trajectory at $s<0$ should "know" its behavior in
the cross channel and vice versa. Most of the familiar meson and baryon
trajectories follow the above regularity: with their parameters fitted in the
scattering region they fit masses and spins of relevant resonances, see e.g.
Ref. [2]. The behavior of trajectories both in the scattering and particle
region is close to linear, which is an approximation to reality. Resonances on
real and linear trajectories imply unrealistic infinitely narrow resonances.
Analyticity and unitarity also require that trajectories be non-linear complex
functions [3, 4]. Constraints on the threshold- and asymptotic behavior of
Regge trajectories were derived from dual amplitudes with Mandelstam
analyticity [4]. Accordingly, near the threshold (see also Refs. [5, 7, 6])

$\Im\alpha(s)_{s\rightarrow s_{0}}\sim(s-s_{0})^{\Re\alpha(s_{0})+1/2},$ (1)

while asymptotically the trajectories are constrained by [4]

$\bigg{|}\frac{\alpha(s)}{\sqrt{s}\ln{s}}\bigg{|}_{s\rightarrow\infty}\leq{\rm
const}.$ (2)

The above asymptotic constrain can be still lowered to a logarithm by imposing
(see Ref. [8] and earlier references) wide-angle power behavior for the
amplitude.

The above constrains are restrictive but still leave much room for model
building. In Refs. [10, 9] the imaginary part of the trajectories
(resonances’ widths) was recovered form the nearly linear real part of the
trajectory by means of dispersion relations and fits to the data.

While the parameters of meson and baryon trajectories can be determined both
from the scattering data and from the particles spectra, this is not true for
the pomeron (and odderon) trajectory, known from fits to scattering data only
(negative values of its argument). An obvious task is to extrapolate the
pomeron trajectory from negative to positive values to predict glueball states
at $J=2,4,...$ non has been found. Given the nearly linear form of the pomeron
trajectory, known from the fits to the (exponential) diffraction cone, little
room is left for variations in the region of particles ($s>0$.) The non-
observability of any glueball state in the expected values of spin and mass
may have two explanations: 1. glueballs appear as hybrid states mixed with
quarks that makes their identification difficult; 2. their production cross
section is low and their widths is large. To resolve these problems one needs
a reliable model to predict cross sections and decay width of the expected
glueballs, in which the pomeron trajectory plays a crucial role.

Models for the pomeron/glueball trajectories were proposed and discussed in
quite a number of paper [11, 12, 13, 14]. They range from simple
phenomenological (also linear) models to quite sophisticated ones, involving
QCD, lattice calculations, extra dimensions, etc. The basic problem of
production cross section and the decay width of produced glueballs in the
cited papers remains open. Close to the spirit of the present approach are
papers [12, 13, 14], where the pomeron/glueball trajectory, including
threshold singularities is manifestly non-linear and the real part terminates.

Fig. 2.: Pomeron-pomeron total cross section in CED, calculated in Ref. [15].

We continue the lines of research initiated in Refs. [15, 1] in which an
analytic pomeron trajectory was used to calculate the pomeron-pomeron cross
section in central exclusive production measurable in proton proton scattering
e.g. at the LHC. The basic idea in that approach is the use of a non-linear
complex Regge trajectory for the pomeron satisfying the requirements of the
analytic $S$-matrix theory and fitting the data. Fits imply high-energy
elastic proton-proton scattering data. For the scattering amplitude the simple
and efficient Donnachie-Landshoff model [16] was used. The essential
difference with respect to many similar studies lies in the non-linear
behavior of the trajectories. They affect crucially the predicted properties
of the resonances. Our previous papers [15, 1] contains more than that: the
fitted trajectories are used to calculate pomeron-pomeron scattering cross
sections in central exclusive diffraction at the LHC. Fig. 2 shows the result
of those calculations.

Papers [15, 1] contain detailed analyses and fits of both the pomeron and
non-leading (also complex!) Regge trajectories, the emphases being on the
pomeron/gluon one. In the present study we revise the basic object, namely the
model of the pomeron trajectory, postponing other details (secondary reggeons,
CED etc.) to a forthcoming study.

###  2.1 Scattering amplitude, cross sections, resonances

In Ref. [1] the resonances contribution to pomeron-pomeron (PP) cross section
was calculated from the imaginary part of the amplitude by use of the optical
theorem

$\sigma_{t}^{PP}(M^{2})\;\;=\;\;{\Im m\;A}(M^{2},t=0)\;\;=$
$\;\;=a\sum_{i=f,P}\sum_{J}\frac{[f_{i}(0)]^{J+2}\;\Im
m\;\alpha_{i}(M^{2})}{(J-\Re e\;\alpha_{i}(M^{2}))^{2}+(\Im
m\;\alpha_{i}(M^{2}))^{2}},$ (3)

In this Section we concentrate on the pomeron. In this case Eq. (3) reduces to

$\sigma_{t}^{PP}(M^{2})\;=\;a\sum_{J}\frac{k^{J+2}\;\Im
m\;\alpha_{P}(M^{2})}{(J-\Re e\;\alpha_{P}(M^{2}))^{2}+(\Im
m\;\alpha_{P}(M^{2}))^{2}},$ (4)

where $k=f_{P}(0)$, and, for simplicity here we set $k=1$.

We start by comparing the resulting glueball spectra in two ways: first we
plot the real and imaginary parts of the trajectory (Chew-Frautchi plot) and
calculate the resonances’ widths by using the relation (see: e.g. Eq. (18) in
Ref. [9])

$\Gamma(s=M^{2})=\frac{2\Im\alpha(s)}{|\alpha^{\prime}(s)|},$ (5)

where $\alpha^{\prime}(s)=dRe\alpha(\sqrt{s})/d\sqrt{s}$.

###  2.2 Analytic Regge trajectories

In previous studies [15, 18, 1] the following two types of trajectories were
considered:

$\alpha(s)=\alpha_{0}+\alpha_{1}s+\alpha_{2}(\sqrt{s_{0}-s}-{\sqrt{s_{0}}}),$
(6)

and

$\alpha(s)=\alpha_{0}+\alpha_{2}(\sqrt{s_{0}-s}-{\sqrt{s_{0}}})+\alpha_{3}(\sqrt{s_{1}-s}{-\sqrt{s_{1}}}),$
(7)

In trajectory (Eq. 7) the second, heavy threshold was introduced to mimic the
nearly linear rise of the trajectory for $s<s_{1}$, avoiding indefinite rise
as in Eq. (6), thus securing the asymptotic square-root upper bound Eq. (2).
As realized in Refs. [1, 15], these trajectories result in "narrowing"
resonances (here, glueball) whose widths decrease as their masses increase.
Below we show that this deficiency is remedied in a trajectory that satisfies
the constraint of the analytic $S$ matrix theory, namely threshold behavior
and asymptotic boundedness, yet produces fading resonances (glueballs), whose
widths are rising with mass.

The trajectory is:

$\alpha(s)=\frac{a+bs}{1+c(\sqrt{s_{0}-s}-\sqrt{s_{0}})},$ (8)

where $s_{0}=4m_{\pi}^{2}$ and $a,b,c$ are adjustable parameters, to be fitted
to scattering ($s<0$) data with the obvious constrains: $\alpha(0)\approx
1.08$ and $\alpha^{\prime}(0)\approx{0.3}$. Trajectory Eq. (8) has square-root
asymptotic behavior, in accord with the requirements of the analytic
$S$-matrix theory.

With the parameters fitted in the scattering region, we continue trajectory
Eq. (8) to positive values of $s$. When approaching the branch cut at
$s=s_{0}$ one has to chose the right Riemann sheet. For $s>s_{0}$ trajectory
(8) may be rewritten as

$\alpha(s)=\frac{a+bs}{1-c(i\sqrt{s-s_{0}}+\sqrt{s_{0}})},$ (9)

with the sign "minus" in front of $c$, according to the definition of the
physical sheet.

For $s>>s_{0},\ \ |\alpha(s)|\rightarrow\frac{b}{c}\sqrt{|s|}$. For $s>s_{0}$
(on the upper edge of the cut), $\Im\alpha(s)>0.$

The intercept is $\alpha(0)=a$ and the slope at $s=0$ is

$\alpha^{\prime}(0)=b+\frac{ac}{2\sqrt{s_{0}}}.$ (10)

To anticipate subsequent fits and discussions, note that the presence of the
light threshold $s_{0}=4m_{\pi}^{2}$ (required by unitarity and the observed
"break" in the data [17]) results in the increasing, compared with the
"standard" value of $\approx 0.25$ GeV-2, intercept.

###  2.3 Simple Regge-pole fits to high-energy elastic scattering data

High-energy elastic proton-proton and proton-antiproton scattering, including
ISR and LHC energies was successfully fitted with non-linear pomeron
trajectories Eqs. (6) and (7) in number of paper, see Ref. [17] and references
therein. Since here we are interested in the parametrization of the pomeron
trajectory, dominating the LHC energy region, we concentrate on the LHC data,
where secondary trajectories can be completely ignored in the near forward
direction.

While at lower energies, e.g. at the ISR, the diffraction cone shows almost
perfect exponential behavior corresponding to a linear pomeron trajectory in a
wide span of $0<-t<1.3$ GeV2, violated only by the "break" near $t\approx-0.1$
GeV2, at the LHC it is almost immediately followed by another structure,
namely by the dip at $t\approx-0.6$ GeV2. The dynamic of the dip (diffraction
minimum) has been treated fully and successfully [18], however those details
are irrelevant to the behavior of the pomeron trajectory in the resonance
(positive $s$) region and expected glueballs there, that depend largely on the
imaginary part of the trajectory and basically on the threshold singularity in
Eq. (8).

In Fig. 3 we show a fit to the low-$|t|$ elastic proton-proton differential
cross section data [19] at 13 TeV with a simple model:

$A_{P}(s,t)=a_{P}e^{b_{P}t}e^{-i\pi\alpha_{P}(t)/2}(s/s_{0P})^{\alpha_{P}(t)},$
(11)

where $\alpha_{P}(t)$ is given by Eq. (8) (changing variable $s$ to variable
$t$). We used the norm:

$\frac{d\sigma}{dt}=\frac{\pi}{s^{2}}|A_{P}(s,t)|^{2}.$ (12)

Fig. 3.: Fitted $pp$ differential cross section at 13 TeV using the amplitude
Eq. (11) and trajectory Eq. (8).

Fig. 4 shows the normalized form of the differential cross section (used by
TOTEM [19]) illustrating the low-$|t|$ "break" phenomenon [17] related to the
non-linear square root term in the pomeron trajectory. However, it should be
also noted that the "break" may result from the two-pion threshold both in the
trajectory and the non-exponential residue, as discussed in Ref. [17].

Fig. 4.: Normalized form of the fitted $pp$ differential cross section at 13
TeV using the amplitude Eq. (11) and trajectory Eq. (8).

###  2.4 Extrapolating the pomeron trajectory to the resonance region, $s>0$

Fitting to the measured $pp$ scattering data, the values of the pomeron
trajectory parameters became known. Changing back variable $t$ to variable $s$
(crossing symmetry), now we can extrapolate the pomeron trajectory to the
resonance region, $s>0$. Figs. 5 and 6 show, respectively, the real and
imaginary parts of the trajectories (during the calculations, the trajectory
parameter values are taken from the fit shown in Fig. 3). Fig. 5 shows the
glueball spectra lying on the pomeron trajectory. Such glueballs have even
integer spins ($J\equiv\Re\alpha_{P}(s)$ = 2, 4, 6…) and mass square
$M^{2}=s$.

In Fig. 8 and Fig. 7 we can see, respectively, the resonance width and the
pomeron component of pomeron-pomeron total cross section.

Fig. 5.: Real part of the pomeron trajectory Eq. (8) as function of $s$. Fig.
6.: Imaginary part of the pomeron trajectory Eq. (8) as function of $s$.  Fig.
7.: Resonance width Eq. (5) calculated with the trajectory Eq. (8). Fig. 8.:
Pomeron component of the pomeron-pomeron total cross section Eq. (4) (setting
$a$ = 1 GeV-2 and $J\in(2,4,6,8,10,12)$) calculated with the trajectory Eq.
(8). The ratios of neighboring resonances’ widths are also shown.

##  3 Summary

Using a simple pomeron pole model fit to the 13 TeV $pp$ low-$|t|$
differential cross section data we have extrapolated the pomeron trajectory
from negative to positive values to predict glueball states at $J=$ 2, 4, 6,
8, 10 and 12. We have predicted also the production cross section and the
decay widths of the expected glueballs. Applying the pomeron trajectory Eq.
(8) we have obtained such resonances (glueballs) whose widths increase as
their masses increase.

Acknowledgment. We thank the organizers of the New Trends in High-Energy
Physics 2019 for their support and the inspiring discussions provided by the
Conference. We thank also László Jenkovszky for his guidance during the
preparation of the manuscript. The work of I. Szanyi was supported by the
"Márton Áron Szakkollégium" program.

## References

  * [1] R. Fiore, L. Jenkovszky and R. Schicker, Exclusive diffractive resonance production in proton-proton collisions at high energies, Eur. Phys. J. C 78, 468 (2018).
  * [2] P. D. B. Collins, An Introduction to Regge Theory and High Energy Physics, Cambridge University Press, (1977).
  * [3] A. Degasperis and E. Predazzi, Dynamical calculation of Regge trajectories, Nuovo Cimento A 65, 764 (1970).
  * [4] A.I. Bugrij, G. Cohen-Tannoudji, L. L. Jenkovszky and N. A. Kobylinsky, Dual Amplitudes with Mandelstam Analyticity, Fortschritte der Physik 21, 427 (1973).
  * [5] A.O. Barut and D.E. Zwanziger, Complex Angular Momentum in Relativistic S-Matrix Theory, Phys. Rev. 127, 974 (1962).
  * [6] V.N. Gribov and Ya.I. Pomeranchuk, Properties of the elastic scattering amplitude at high energies, Nucl. Phys. 38, 516 (1962).
  * [7] R. Oehme, Complex angular momentum in elementary particle scattering, in Proceedings of the Scottish Universities Summer School Edinburgh 1963, edited by R. G. Moorhouse p.129 (1964).
  * [8] L.L. Jenkovszky, High-energy elastic hadron scattering, La Rivista del Nuovo Cimento 10, 1 (1987).
  * [9] R. Fiore, L. L. Jenkovszky, F. Paccanoni and A. Prokudin, Baryonic Regge trajectories with analyticity constraints, Phys. Rev. D 70, 054003 (2004), arXiv:hep-ph/0404021.
  * [10] R. Fiore, L.L. Jenkovszky, V. Magas, F. Paccanoni and A. Papa, Analytic model of Regge trajectories, Eur.Phys.J. A 10, 217 (2001), arXiv:hep-ph/0011035.
  * [11] M. N. Sergeenko, Gluonium states and the Pomeron trajectory, arXiv:0807.0911.
  * [12] M.M. Brisudova, L. Burakovsky, T. Goldman and A. Szczepaniak, Nonlinear Regge trajectories and glueballs, arXiv:nucl-th/0303012.
  * [13] T. Goldman and A. Szczepaniak, Nonlinear Regge trajectories and glueballs, arXiv:hep-ph/0009126.
  * [14] M.M. Brisudova, L. Burakovsky and T. Goldman, Spectroscopy "windows" of quark-antiquark mesons and glueballs with effective Regge trajectories, arXiv:hep-ph/0008212.
  * [15] R. Fiore, L. Jenkovszky and R. Schicker, Resonance production in PomeronPomeron collisions at the LHC, Eur. Phys. J. C 76, 38 (2016).
  * [16] A. Donnachie and P.V. Landshoff, Hard diffraction: Production of high pT jets, W or Z, and Drell-Yan pairs, Nucl. Phys. B 303, 634 (1988).
  * [17] L. Jenkovszky, I. Szanyi and CI. Tan, Shape of proton and the pion cloud, Eur. Phys. J. A 54, 116 (2018).
  * [18] I. Szanyi, N. Bence and L. Jenkovszky, New physics from TOTEM’s recent measurements of elastic and total cross sections, J. Phys. G: Nucl. Part. Phys. 46, 055002 (2019).
  * [19] The TOTEM collaboration, First determination of the $\rho$ parameter at $\sqrt{s}=13TeV$ – probing the existence of a colourless three-gluon bound state, arXiv:1812.04732.


# Drag of Two Cylindrical Intruders in a Two-Dimensional Granular Environment

Takumi Kubota1 Haruto Ishikawa1 and Satoshi Takada1,2 Corresponding author.
E-mail<EMAIL_ADDRESS>of Mechanical Systems
Engineering1Department of Mechanical Systems Engineering Tokyo University of
Agriculture and Technology Tokyo University of Agriculture and Technology
Koganei Koganei Tokyo 184–8588 Tokyo 184–8588 Japan
2Institute of Engineering Japan
2Institute of Engineering Tokyo University of Agriculture and Technology Tokyo
University of Agriculture and Technology Koganei Koganei Tokyo 184–8588 Tokyo
184–8588 Japan Japan

###### Abstract

The drag of two cylindrical intruders in a two-dimensional granular
environment is numerically studied by the discrete element method. We find the
yield force, below which the intruders cannot move because of interactions
with the surrounding particles. Above the yield force, on the other hand, the
intruders can move at a constant speed. We investigate the relationship
between the drag force and the steady speed of the intruders, where the speed
becomes higher as the distance between the intruders decreases. We confirm
that the origin of the yield is the Coulombic friction between the particles
and the bottom plate by changing the value of the friction coefficient. We
also find that the yield force is almost proportional to the friction
coefficient, which means that the number of particles determining the yield
force is almost constant. On the other hand, the two-dimensional elasticity is
applicable to determine the stress fields around the intruders. We confirm
that fields asymmetric with respect to the drag direction are reproduced by
using the information of the stresses on the surfaces of the intruders by
introducing bipolar coordinates.

##  1 Introduction

An understanding of the resistance to intruders into a flow is important in
many applications [(1), (2), (3)]. The drag in a fluid is one of the principal
problems that has been studied for a long time [(2)], where the drag is
proportional to the speed or the square of the speed if the flow speed is low
or high, respectively. It is known that the drag is characterized by the
Reynolds number, which is the ratio of the flow speed to the viscosity of the
fluid.

On the other hand, the drag in granular environments changes depending on the
protocol of the experiments or simulations [(11), (26), (12), (13), (14),
(20), (8), (9), (10), (5), (6), (16), (17), (18), (19), (4), (21), (22), (7),
(23), (24), (25), (15)] using spheres [(4), (5), (6), (7)], cylinders [(11),
(8), (9), (10)], or other shapes of objects [(12), (13), (14)]. One of the
simplest setups was numerically prepared by Hilton and Tordesillas [(4)]. They
pulled a spherical intruder in the direction perpendicular to gravity. In
their setup, the drag is given by the combination of a constant term and a
term proportional to the drag speed. Another important setup was studied by
Takehara and co-workers [(8), (9)]. They reported that the drag of an intruder
in a two-dimensional environment is the sum of a term independent of the drag
speed and one depending on the square of the speed. This was also numerically
verified by Takada and Hayakawa in a similar setup [(10)] where the intruder
was pulled in a direction perpendicular to gravity.

We also know that intruders cannot move owing to the existence of surrounding
particles if the drag force is sufficiently small. The threshold is related to
the friction between particles or between particles and the bottom [(10)].
Although there are few studies on the force below the threshold, it is
important to understand this regime because there are fracture and other
phenomena. Once this regime is understood, we can expect to be able to predict
the minimum drag force required to move an object, which will yield a
reduction of energy loss. In addition, knowledge of this phenomena may be
applied to the application of fracture mechanics or plastic mechanics to
granular systems. Referring to those studies, we previously investigated the
drag of the intruder using a similar setup to check how the stress fields
around the intruder. We found good agreement between the stress fields
obtained by simulations and those based on the theory of elasticity [(15)]. We
note that the drag problems have also been studied in the context of planetary
science [(27)].

Most previous studies focus on the drag of a single intruder. However, it is
important to investigate the drag of multiple intruders [(7), (22), (21),
(28)] when we consider more realistic systems such as the movement of a fork
hoe. In this case, the existence of one intruder must affect the other;
therefore, the drag should be different from that in the one-intruder cases.
To clarify this, we study the drag of multiple intruders in a two-dimensional
granular environment. As a first step, we consider two identical cylindrical
objects vertically aligned to the drag direction. In the following, the drag
law and the stress around the intruders is investigated by the discrete
element method (DEM) [(29)].

This paper is organized as follows. In the next section, we explain the model
and the setup of our simulations. Section 3 is the main part of this paper, in
which we show the relationship between the drag force and the steady speed of
the intruders. We also reproduce the stress fields around the intruders from
the information on the surface of the intruders. In Sects. 4 and 5, we discuss
and conclude our results. In Appendix A, we explain how the system size is
determined in the simulations. In Appendix B, we examine the drag law when we
increase the number of intruders. In Appendix C, we present the general
framework to obtain the stress fields in terms of the bipolar coordinates.

##  2 Model and Setup of Our Simulation

Figure 1: (Color online) Typical snapshot of the simulation. Two intruders are
placed with interval $\Delta$. They are pulled with force $F_{\rm drag}$ in
the $x$-direction. Here, gravity acts perpendicular to the $xy$ plane.

In this section, we explain the model and the setup of our simulations.
Initially, two identical intruders are placed together at
$(x,y)=(0,\pm\Delta/2)$, as shown Fig. 1. Here, the intruders are
characterized by diameter $D$ and mass $M$, and the distance between them is
fixed. Then, smaller particles are placed around the intruders at random.
These particles are characterized by diameter $d_{i}$ and mass $m_{i}$. To
avoid crystallization, the surrounding particles are prepared with the number
ratio $1$$:$$1$ and the size ratio $d$$:$$1.4d$. We note that the mass ratio
is given by $m$:$1.4^{2}m$. Here, the mass density of the intruders is assumed
to be the same as that of the surrounding particles. The equations of motion
of the intruders and the particles are, respectively, described by [(29)]

$\begin{cases}\displaystyle
M\frac{d^{2}\bm{R}_{i}}{dt^{2}}\cdot\hat{\bm{e}}_{x}=\sum_{j}\bm{F}_{ij}\cdot\hat{\bm{e}}_{x}+F_{\rm
drag},\\\ \displaystyle
m_{i}\frac{d^{2}\bm{r}_{i}}{dt^{2}}=\sum_{j}\bm{F}_{ij}-\mu_{\rm
b}m_{i}g\hat{\bm{v}}_{i},\end{cases}$ (1)

where $\bm{R}_{i}$ and $\bm{r}_{i}$ are the position vectors of the $i$-th
intruder and the $i$-th particle, respectively, $\bm{F}_{ij}$ is the contact
force between the $i$-th and the $j$-th intruders or particles,
$F_{\rm{drag}}$ is the drag force acting on the intruders, $\hat{\bm{e}}_{x}$
is the unit vector parallel to the $x$-direction, $\mu_{\rm b}$ is a
coefficient of Coulombic friction acting on the surrounding particles and the
bottom plate, $g$ is the acceleration due to gravity, and $\hat{\bm{v}}_{i}$
is the unit vector parallel to the velocity of the $i$-th particle [(30)].
Here, we ignore the motions of the intruders in the $y$-direction. Hereafter,
for simplicity, we do not distinguish $\bm{R}_{i}$ from $\bm{r}_{j}$ when we
consider the positions of the particles. The contact force $\bm{F}_{ij}$ is
described by the sum of the normal force $F^{n}_{ij}$ and the tangential force
$F^{t}_{ij}$, i.e., $\bm{F}_{ij}=F^{n}_{ij}\bm{n}+F^{t}_{ij}\bm{t}$, where
$\bm{n}$ and $\bm{t}$ are unit vectors in the normal and tangential directions
between the two contacting particles, respectively. The force in the normal
direction $F^{n}_{ij}$ is given by

$F^{n}_{ij}=\left[k_{n}\delta-\eta_{n}\left(\bm{v}_{ij}\cdot\hat{\bm{r}}_{ij}\right)\right]\Theta\left(d_{ij}-r_{ij}\right),$
(2)

where $k_{n}$ and $\eta_{n}$ are the spring constant and the magnitude of the
dashpot in the normal direction, respectively, $\delta$ is the overlap of the
two contacting particles, $\bm{v}_{ij}$ is the relative velocity described by
$\bm{v}_{ij}=\bm{v}_{i}-\bm{v}_{j}$, $d_{ij}=(d_{i}+d_{j})/2$ is the sum of
the radii of the two contacting particles,
$r_{ij}=|\bm{r}_{ij}|=|\bm{r}_{i}-\bm{r}_{j}|$ is the relative displacement,
$\hat{\bm{r}}_{ij}=\bm{r}_{ij}/r_{ij}$, and $\Theta(x)$ is the step function.
Here, it is well known that the viscoelastic force becomes linear with a
logarithmic correction, which means that the linear model is better than the
Hertzian model [(31), (32)]. Next, the force in the tangential direction
$F^{t}_{ij}$ is given by

$F^{t}_{ij}=\min\left(\mu
F^{n}_{ij},\left|-k_{t}\bm{\xi}-\eta_{t}\bm{v}^{t}_{ij}\right|\right)\Theta\left(d_{ij}-r_{ij}\right),$
(3)

where $\mu$ is the coefficient of Coulombic friction acting between the
particles, $k_{t}$ and $\eta_{t}$ are the spring constant and the magnitude of
the dashpot in the tangential direction, and
$\bm{v}^{t}_{ij}=\bm{v}_{ij}-(\bm{v}_{ij}\cdot\hat{\bm{r}}_{ij})\hat{\bm{r}}_{ij}$
and $\bm{\xi}=\int^{t}\bm{v}^{t}_{ij}(t^{\prime})dt^{\prime}$ are the
tangential relative velocity and the tangential displacement from the
beginning of the contact, respectively [(30)]. In the simulation, we set
$D=10d$, number of the surrounding particles $N=4000$, packing fraction
$\varphi=0.80$, $k_{t}=k_{n}$, $\eta_{t}=\eta_{n}$, $\mu=\mu_{\rm b}=0.20$,
and $g=1.0\times 10^{-4}k_{n}d/m$. For our choice of parameters, the
restitution coefficients in the normal and tangential directions correspond to
$e_{n}=e_{t}=4.3\times 10^{-2}$. Here, the determination of the system sizes
($L_{x}$ and $L_{y}$) is discussed in Appendix A. We apply the periodic
boundary condition in the $x$-direction, and we put the smaller particles at
$y=\pm L_{y}/2$ as the bumpy boundary condition in the $y$-direction [(10)].
In this paper, we mainly focus on the case where there exist two intruders in
the system. The results for three intruders are briefly discussed in Appendix
B.

##  3 Results

In this section, we explain the results obtained from the simulation in two
parts. In the first subsection, we show the drag law when the drag force is
sufficiently larger than the yield force, which will be defined there. In the
next subsection, we focus on the stress fields when the drag force is smaller
than but close to the yield force. We try to reproduce the stress fields
around the intruders in terms of the Airy stress function known in elasticity.

Figure 2: (Color online) Steady speed of the intruders against the drag force
for various distances $\Delta$ and the packing fraction $\varphi$ for $D=10d$.
Here, we have introduced $F_{\rm drag}^{*}\equiv F_{\rm drag}/(k_{n}d)$ and
$V^{*}\equiv V/(d\sqrt{k_{n}/m})$. The solid line represents Eq. (4) for the
single-intruder case.

###  3.1 Drag law

First, we observe the steady speed $V$ when we change the drag force $F_{\rm
drag}$, the distance of the intruders $\Delta$, and the packing fraction
$\varphi$. Note that the rotational motion of the intruders is not observed.
As shown in Fig. 2, there is a yield force $F_{\rm Y}$, and the value of
$F_{\rm Y}$ is found to be almost unchanged for the distance $\Delta$,
although it increases as the packing fraction $\varphi$ becomes larger. We
note that the velocity dependence is almost the same. Now, we focus on the
dependence of the drag on the distance $\Delta$.

For $F_{\rm drag}<F_{\rm Y}$, the intruders cannot move because the drag force
and the frictional forces of the surrounding particles are balanced. From the
simulation results, the value of the yield force becomes $F_{\rm Y}\simeq
1.8\times 10^{-2}k_{n}d$ for $\mu_{\rm b}=0.20$. Above the yield force $F_{\rm
Y}$, on the other hand, the intruders can move with a steady speed $V$, where
$V$ is an increasing function depending on $F_{\rm drag}$. Note that the rate
of increase of $V$ depends on the distance $\Delta$. Here, in the case of a
single intruder, we previously reported that the relationship between the drag
force $F_{\rm drag}$ and the steady speed $V$ is given by [(15)]

$F_{\rm drag}=F_{\rm
Y}+\frac{2}{\pi}\varphi\frac{Mm}{M+m}\frac{D+d}{d^{2}}V^{2}\left(\frac{14}{9}+\frac{4}{3}e_{n}+\frac{2}{9}e_{t}\right),$
(4)

which is derived under the approximation that the surrounding particles are
monodisperse and the collisions are characterized by the restitution
coefficients $e_{n}$ and $e_{t}$ (see, e.g., Ref. 33). When there exist two
intruders, the steady speed becomes higher as the interval between the
intruders decreases, which was also reported in Ref. 7. This is because the
area pushed by the upper intruder has a part that is common with that of the
lower intruder. This fact is also understood from Fig. 3, which shows the
displacements of particles with the time interval $\Delta t=50\sqrt{m/k_{n}}$
when we change the interval between the intruders, $\Delta$. It is clear that
the particles between the intruders are pushed by both intruders. Note that
the drag might change when the distance between the intruders becomes smaller
as $\Delta-D\gtrsim d$ because we keep its distance unchanged. For the
sufficiently small interval, the surrounding particles cannot pass between the
intruders Dhiman20 .

Figure 3: (Color online) Displacements of particles with the time interval
$\Delta t=50\sqrt{m/k_{n}}$ in the steady state for $\varphi=0.80$, $D=10d$,
and $F_{\rm drag}=2.5\times 10^{-2}k_{n}d$. Panels (a) and (b) correspond to
$\Delta=15d$ and $20d$, respectively. The color represents the angle of the
displacements. The line width is propotional to the displacement. Figure 4:
(Color online) Yield force of the intruders against frictional coefficient
$\mu_{\rm b}$. The dashed line represents $F_{\rm Y}=9.2\times 10^{2}\mu_{\rm
b}mg$.

Next, we observe the yield force $F_{\rm Y}$ of the intruders for various
bottom friction coefficients $\mu_{\rm b}$ when we fix $\Delta=15d$. The yield
force $F_{\rm Y}$ is related to the frictional forces of the surrounding
particles interacting with the intruders. This means that the yield force
should be given by

$F_{\rm Y}=\mu_{\rm b}M_{\rm Y}g,$ (5)

where $M_{\rm Y}$ is the total mass of all particles interacting with the
intruders. Figure 4 indicates that the above interpretation is correct because
the ratio of the yield force to the frictional force is almost constant.

###  3.2 Stress fields around intruders

For $F_{\rm drag}<F_{\rm Y}$, let us investigate the stress fields around the
intruders. Here, the stress obtained from the simulation is given by

$\sigma_{\alpha\beta}=\frac{1}{L_{x}L_{y}}\sum_{i}\left(m_{i}v_{i,\alpha}v_{i,\beta}+\frac{1}{2}\sum_{j}r_{ij,\alpha}F_{ij,\beta}\right),$
(6)

where we have used the coarse-grained method for smoothing the data with the
characteristic length $w=2d$Zhang10 .

We introduce the Airy stress function for comparison with the theory of
elasticity Timoshenko ; Daniels17 . For this purpose, let us expand
$\sigma_{\varrho\varrho}$, $\sigma_{\varrho\vartheta}$, and
$\sigma_{\vartheta\vartheta}$ in the Fourier series as

$\begin{Bmatrix}\sigma_{\varrho\varrho}\\\ \sigma_{\varrho\vartheta}\\\
\sigma_{\vartheta\vartheta}\end{Bmatrix}=\frac{1}{2}\begin{Bmatrix}\widetilde{\sigma}_{\varrho\varrho}^{{\rm
C}(0)}\\\ \widetilde{\sigma}_{\varrho\vartheta}^{{\rm C}(0)}\\\
\widetilde{\sigma}_{\vartheta\vartheta}^{{\rm
C}(0)}\end{Bmatrix}+\sum_{n=1}^{8}\left[\begin{Bmatrix}\widetilde{\sigma}_{\varrho\varrho}^{{\rm
C}(n)}\\\ \widetilde{\sigma}_{\varrho\vartheta}^{{\rm C}(n)}\\\
\widetilde{\sigma}_{\vartheta\vartheta}^{{\rm
C}(n)}\end{Bmatrix}\cos(n\vartheta)+\begin{Bmatrix}\widetilde{\sigma}_{\varrho\varrho}^{{\rm
S}(n)}\\\ \widetilde{\sigma}_{\varrho\vartheta}^{{\rm S}(n)}\\\
\widetilde{\sigma}_{\vartheta\vartheta}^{{\rm
S}(n)}\end{Bmatrix}\sin(n\vartheta)\right].$ (7)

In the Fourier series, we consider terms up to the eighth order. As shown in
Fig. 5, the Fourier series up to the eighth order well reproduce those
obtained from the simulations without any fitting parameter. We mainly focus
on the upper intruder throughout this paper because the stress fields are
symmetric with respect to the $x$-axis.

Figure 5: (Color online) Profiles of the stresses on the surfaces of the upper
intruders for $D=10d$, $\Delta=20d$, and $F_{\rm drag}=1.6\times
10^{-2}k_{n}d$, where we have introduced
$\sigma_{\alpha\beta}^{*}=\sigma_{\alpha\beta}/(F_{\rm drag}/(\pi d))$ with
$(\alpha,\beta)=(\varrho,\varrho)$, $(\varrho,\vartheta)$, and
$(\vartheta,\vartheta)$. The marks represent the stresses obtained from the
simulations, and the lines represent the Fourier series expansions of the
stresses.

Now, let us construct the stress field around the intruders in terms of the
surface stress. As discussed in Appendix C, it is convenient to adopt the
bipolar coordinates as

$\begin{cases}\displaystyle
x=\frac{\sqrt{\Delta^{2}-D^{2}}}{2}\frac{\sin\vartheta}{\cosh\varrho-\cos\vartheta}\\\
\displaystyle
y=\frac{\sqrt{\Delta^{2}-D^{2}}}{2}\frac{\sinh\varrho}{\cosh\varrho-\cos\vartheta}\end{cases}.$
(8)

Using the same procedure as in Ref. 15, the stress field can be constructed in
terms of the Airy stress function $\chi$ as Jeffery21

$\displaystyle\sigma_{\varrho\varrho}$
$\displaystyle=\frac{2}{\sqrt{\Delta^{2}-D^{2}}}\left[(\cosh\varrho-\cos\vartheta)\frac{\partial^{2}}{\partial\vartheta^{2}}-\sinh\varrho\frac{\partial}{\partial\varrho}\right.$
$\displaystyle\hskip
60.00009pt\left.-\sin\vartheta\frac{\partial}{\partial\vartheta}+\cosh\varrho\right]\left(\frac{\chi}{J}\right)$
(9a) $\displaystyle\sigma_{\vartheta\vartheta}$
$\displaystyle=\frac{2}{\sqrt{\Delta^{2}-D^{2}}}\left[(\cosh\varrho-\cos\vartheta)\frac{\partial^{2}}{\partial\varrho^{2}}-\sinh\varrho\frac{\partial}{\partial\varrho}\right.$
$\displaystyle\hskip
60.00009pt\left.-\sin\vartheta\frac{\partial}{\partial\vartheta}+\cos\vartheta\right]\left(\frac{\chi}{J}\right)$
(9b) $\displaystyle\sigma_{\varrho\vartheta}$
$\displaystyle=-\frac{2}{\sqrt{\Delta^{2}-D^{2}}}(\cosh\varrho-\cos\vartheta)\frac{\partial^{2}}{\partial\varrho\partial\vartheta}\left(\frac{\chi}{J}\right),$
(9c)

where

$J=\frac{\sqrt{\Delta^{2}-D^{2}}}{2(\cosh\varrho-\cos\vartheta)},$ (10)

and

$\displaystyle\frac{\chi}{J}$
$\displaystyle=(E\cos\vartheta+F\sin\vartheta+G\cosh\varrho+H\sinh\varrho)\vartheta$
$\displaystyle\hskip 10.00002pt+(B_{0}\cosh\varrho+D_{0}\sinh\varrho)\varrho$
$\displaystyle\hskip
10.00002pt+[A_{1}\cosh(2\varrho)+B_{1}+C_{1}\sinh(2\varrho)+D_{1}\varrho]\cos\vartheta$
$\displaystyle\hskip
10.00002pt+[A_{1}^{\prime}\cosh(2\varrho)+C_{1}^{\prime}\sinh(2\varrho)+D_{1}^{\prime}\varrho]\sin\vartheta$
$\displaystyle\hskip
10.00002pt+\sum_{n=2}^{8}\left\\{A_{n}\cosh[(n+1)\varrho]+B_{n}\cosh[(n-1)\varrho]\right.$
$\displaystyle\hskip
20.00003pt\left.+C_{n}\sinh[(n+1)\varrho]+D_{n}\sinh[(n-1)\varrho]\right\\}\cos(n\vartheta)$
$\displaystyle\hskip
10.00002pt+\sum_{n=2}^{8}\left\\{A_{n}^{\prime}\cosh[(n+1)\varrho]+B_{n}^{\prime}\cosh[(n-1)\varrho]\right.$
$\displaystyle\hskip
20.00003pt\left.+C_{n}^{\prime}\sinh[(n+1)\varrho]+D_{n}^{\prime}\sinh[(n-1)\varrho]\right\\}\sin(n\vartheta),$
(11)

with constants $E$, $F$, $G$, $H$, $A_{n}$, $B_{n}$, $C_{n}$, $D_{n}$,
$A_{n}^{\prime}$, $B_{n}^{\prime}$, $C_{n}^{\prime}$, and $D_{n}^{\prime}$
($n\geq 0$). Note that higher order terms ($n\geq 9$) are dropped in Eq. (11)
because they do not affect the results. Here, we consider the functional form
of each stress component, $H$, $B_{0}$, $B_{0}^{\prime}$ $C_{n}$, $D_{n}$,
$C_{n}^{\prime}$, and $D_{n}^{\prime}$ to be zero.

To determine each Fourier coefficient, we introduce $35$-dimensional vectors
$\bm{\widetilde{\Sigma}}_{\rm{surf}}=(\widetilde{\sigma}_{\varrho\varrho}^{{\rm
C}(0)},\ \cdots,\ \widetilde{\sigma}_{\varrho\varrho}^{{\rm C}(8)},\
\widetilde{\sigma}_{\varrho\varrho}^{{\rm S}(1)},$ $\cdots,\
\widetilde{\sigma}_{\varrho\varrho}^{{\rm S}(8)},\
\widetilde{\sigma}_{\varrho\vartheta}^{{\rm C}(0)},\ \cdots,\
\widetilde{\sigma}_{\varrho\vartheta}^{{\rm C}(8)},\
\widetilde{\sigma}_{\varrho\vartheta}^{{\rm S}(1)},\ \cdots,\
\widetilde{\sigma}_{\varrho\vartheta}^{{\rm S}(8)},\
\widetilde{\sigma}_{\varrho\varrho}^{{\rm
C}(0)}+\widetilde{\sigma}_{\vartheta\vartheta}^{{\rm C}(0)})^{T}$ and
$\bm{X}=(x_{1},x_{2},\cdots,x_{35})^{T}$. Here, the components of $\bm{X}$ are
given by

$\displaystyle x_{1}=E.\quad x_{2}=F,\quad x_{3}=G,\quad x_{4}=D_{0},$
$\displaystyle x_{5}=A_{1},\ \cdots,\ x_{12}=A_{8},\ x_{13}=A^{\prime}_{1},\
\cdots,\ x_{20}=A^{\prime}_{8},$ $\displaystyle x_{21}=B_{1},\ \cdots,\
x_{28}=B_{8},\ x_{29}=B^{\prime}_{2},\ \cdots,\ x_{35}=B^{\prime}_{8}.$ (12)

Here, the Fourier components are functions of the vector $\bm{X}$, where their
expressions are summarized in Table 1.

Table 1: Coefficients in the Fourier series of $\sigma_{\varrho\varrho}$,
$\sigma_{\varrho\vartheta}$, and $\sigma_{\vartheta\vartheta}$ in Eqs. (9)
with Eq. (11).
$a_{\varrho\varrho}^{(0)}=-3F+D_{0}+2B_{1}-(D_{0}-2A_{1})\cosh(2\varrho)$
$a_{\varrho\varrho}^{(1)}=(2F+A_{1}+3B_{2})\cosh(\varrho)-(A_{1}-3A_{2})\cosh(3\varrho)$
---
$a_{\varrho\varrho}^{(2)}=-\frac{1}{2}F-B_{2}-2(B_{2}-3B_{3})\cosh(2\varrho)-3(A_{2}-2A_{3})\cosh(4\varrho)$
$a_{\varrho\varrho}^{(3)}=(B_{2}-3B_{3})\cosh(\varrho)+(A_{2}-2A_{3}-5B_{3}+10B_{4})\cosh(3\varrho)-2(3A_{3}-5A_{4})\cosh(5\varrho)$
$a_{\varrho\varrho}^{(4)}=3(B_{3}-2B_{4})\cosh(2\varrho)+(3A_{3}-5A_{4}-9B_{4}+15B_{5})\cosh(4\varrho)-5(2A_{4}-3A_{5})\cosh(6\varrho)$
$a_{\varrho\varrho}^{(5)}=2(3B_{4}-5B_{5})\cosh(3\varrho)+(6A_{4}-9A_{5}-14B_{5}+21B_{6})\cosh(5\varrho)-3(5A_{5}-7A_{6})\cosh(7\varrho)$
$a_{\varrho\varrho}^{(6)}=5(2B_{5}-3B_{6})\cosh(4\varrho)+2(5A_{5}-7A_{6}-10B_{6}+14B_{7})\cosh(6\varrho)-7(3A_{6}-4A_{7})\cosh(8\varrho)$
$a_{\varrho\varrho}^{(7)}=3(5B_{6}-7B_{7})\cosh(5\varrho)+(15A_{6}-20A_{7}-27B_{7}+36B_{8})\cosh(7\varrho)-4(7A_{7}-9A_{8})\cosh(9\varrho)$
$a_{\varrho\varrho}^{(8)}=7(3B_{7}-4B_{8})\cosh(6\varrho)+(21A_{7}-27A_{8}-35B_{8})\cosh(8\varrho)-36A_{8}\cosh(10\varrho)$
$b_{\varrho\varrho}^{(1)}=2(E+G)+(-2E-G+A^{\prime}_{1}+3B^{\prime}_{2})\cosh(\varrho)-(A^{\prime}_{1}-3A^{\prime}_{2})\cosh(3\varrho)$
$b_{\varrho\varrho}^{(2)}=-\frac{1}{2}E-G-B^{\prime}_{2}-2(B^{\prime}_{2}-3B^{\prime}_{3})\cosh(2\varrho)-3(A^{\prime}_{2}-2A^{\prime}_{3})\cosh(4\varrho)$
$b_{\varrho\varrho}^{(3)}=\frac{2}{3}(E+G)+(B^{\prime}_{2}-3B^{\prime}_{3})\cosh(\varrho)+(A^{\prime}_{2}-2A^{\prime}_{3}-5B^{\prime}_{3}+10B^{\prime}_{4})\cosh(3\varrho)-2(3A^{\prime}_{3}-5A^{\prime}_{4})\cosh(5\varrho)$
$b_{\varrho\varrho}^{(4)}=-\frac{1}{2}(E+G)+3(B^{\prime}_{3}-2B^{\prime}_{4})\cosh(2\varrho)+(3A^{\prime}_{3}-5A^{\prime}_{4}-9B^{\prime}_{4}+15B^{\prime}_{5})\cosh(4\varrho)-5(2A^{\prime}_{4}-3A^{\prime}_{5})\cosh(6\varrho)$
$b_{\varrho\varrho}^{(5)}=\frac{2}{5}(E+G)+2(3B^{\prime}_{4}-5B^{\prime}_{5})\cosh(3\varrho)+(6A^{\prime}_{4}-9A^{\prime}_{5}-14B^{\prime}_{5}+21B^{\prime}_{6})\cosh(5\varrho)-3(5A^{\prime}_{5}-7A^{\prime}_{6})\cosh(7\varrho)$
$b_{\varrho\varrho}^{(6)}=-\frac{1}{3}(E+G)+5(2B^{\prime}_{5}-3B^{\prime}_{6})\cosh(4\varrho)+2(5A^{\prime}_{5}-7A^{\prime}_{6}-10B^{\prime}_{6}+14B^{\prime}_{7})\cosh(6\varrho)-7(3A^{\prime}_{6}-4A^{\prime}_{7})\cosh(8\varrho)$
$b_{\varrho\varrho}^{(7)}=\frac{2}{7}(E+G)+3(5B^{\prime}_{6}-7B^{\prime}_{7})\cosh(5\varrho)+(15A^{\prime}_{6}-20A^{\prime}_{7}-27B^{\prime}_{7}+36B^{\prime}_{8})\cosh(7\varrho)-4(7A^{\prime}_{7}-9A^{\prime}_{8})\cosh(9\varrho)$
$b_{\varrho\varrho}^{(8)}=-\frac{1}{4}(E+G)+7(3B^{\prime}_{7}-4B^{\prime}_{8})\cosh(6\varrho)+(21A^{\prime}_{7}-27A^{\prime}_{8}-35B^{\prime}_{8})\cosh(8\varrho)-36A^{\prime}_{8}\cosh(10\varrho)$
$a_{\varrho\vartheta}^{(0)}=-(G-2A^{\prime}_{1})\sinh(2\varrho)$
$a_{\varrho\vartheta}^{(1)}=(G-A^{\prime}_{1}+B^{\prime}_{2})\sinh(\varrho)-(A^{\prime}_{1}-3A^{\prime}_{2})\sinh(3\varrho)$
$a_{\varrho\vartheta}^{(2)}=(A^{\prime}_{1}-3A^{\prime}_{2}-B^{\prime}_{2}+3B^{\prime}_{3})\sinh(2\varrho)-3(A^{\prime}_{2}-2A^{\prime}_{3})\sinh(4\varrho)$
$a_{\varrho\vartheta}^{(3)}=(B^{\prime}_{2}-3B^{\prime}_{3})\sinh(\varrho)+3(A^{\prime}_{2}-2A^{\prime}_{3}-B^{\prime}_{3}+2B^{\prime}_{4})\sinh(3\varrho)-2(3A^{\prime}_{3}-5A^{\prime}_{4})\sinh(5\varrho)$
$a_{\varrho\vartheta}^{(4)}=3(B^{\prime}_{3}-2B^{\prime}_{4})\sinh(2\varrho)+2(3A^{\prime}_{3}-5A^{\prime}_{4}-3B^{\prime}_{4}+5B^{\prime}_{5})\sinh(4\varrho)-5(2A^{\prime}_{4}-3A^{\prime}_{5})\sinh(6\varrho)$
$a_{\varrho\vartheta}^{(5)}=2(3B^{\prime}_{4}-5B^{\prime}_{5})\sinh(3\varrho)+5(2A^{\prime}_{4}-3A^{\prime}_{5}-2B^{\prime}_{5}+3B^{\prime}_{6})\sinh(5\varrho)-3(5A^{\prime}_{5}-7A^{\prime}_{6})\sinh(7\varrho)$
$a_{\varrho\vartheta}^{(6)}=5(2B^{\prime}_{5}-3B^{\prime}_{6})\sinh(4\varrho)+3(5A^{\prime}_{5}-7A^{\prime}_{6}-5B^{\prime}_{6}+7B^{\prime}_{7})\sinh(6\varrho)-7(3A^{\prime}_{6}-4A^{\prime}_{7})\sinh(8\varrho)$
$a_{\varrho\vartheta}^{(7)}=3(5B^{\prime}_{6}-7B^{\prime}_{7})\sinh(5\varrho)+7(3A^{\prime}_{6}-4A^{\prime}_{7}-3B^{\prime}_{7}+4B^{\prime}_{8})\sinh(7\varrho)-4(7A^{\prime}_{7}-9A^{\prime}_{8})\sinh(9\varrho)$
$a_{\varrho\vartheta}^{(8)}=7(3B^{\prime}_{7}-4B^{\prime}_{8})\sinh(6\varrho)+4(7A^{\prime}_{7}-9A^{\prime}_{8}-7B^{\prime}_{8})\sinh(8\varrho)-36A^{\prime}_{8}\sinh(10\varrho)$
$b_{\varrho\vartheta}^{(1)}=(A_{1}-B_{2})\sinh(\varrho)+(A_{1}-3A_{2})\sinh(3\varrho)$
$b_{\varrho\vartheta}^{(2)}=-(A_{1}-3A_{2}-B_{2}+3B_{3})\sinh(2\varrho)+3(A_{2}-2A_{3})\sinh(4\varrho)$
$b_{\varrho\vartheta}^{(3)}=-(B_{2}-3B_{3})\sinh(\varrho)-3(A_{2}-2A_{3}-B_{3}+2B_{4})\sinh(3\varrho)+2(3A_{3}-5A_{4})\sinh(5\varrho)$
$b_{\varrho\vartheta}^{(4)}=-3(B_{3}-2B_{4})\sinh(2\varrho)-2(3A_{3}-5A_{4}-3B_{4}+5B_{5})\sinh(4\varrho)+5(2A_{4}-3A_{5})\sinh(6\varrho)$
$b_{\varrho\vartheta}^{(5)}=-2(3B_{4}-5B_{5})\sinh(3\varrho)-5(2A_{4}-3A_{5}-2B_{5}+3B_{6})\sinh(5\varrho)+3(5A_{5}-7A_{6})\sinh(7\varrho)$
$b_{\varrho\vartheta}^{(6)}=-5(2B_{5}-3B_{6})\sinh(4\varrho)-3(5A_{5}-7A_{6}-5B_{6}+7B_{7})\sinh(6\varrho)+7(3A_{6}-4A_{7})\sinh(8\varrho)$
$b_{\varrho\vartheta}^{(7)}=-3(5B_{6}-7B_{7})\sinh(5\varrho)-7(3A_{6}-4A_{7}-3B_{7}+4B_{8})\sinh(7\varrho)+4(7A_{7}-9A_{8})\sinh(9\varrho)$
$b_{\varrho\vartheta}^{(8)}=-7(3B_{7}-4B_{8})\sinh(6\varrho)-4(7A_{7}-9A_{8}-7B_{8})\sinh(8\varrho)+36A_{8}\sinh(10\varrho)$
$a_{\vartheta\vartheta}^{(0)}=-F+2B_{1}+3D_{0}+(D_{0}-2A_{1})\cosh(2\varrho)$
$a_{\vartheta\vartheta}^{(1)}=(-2D_{0}+3A_{1}+B_{2})\cosh(\varrho)+(A_{1}-3A_{2})\cosh(3\varrho)$
$a_{\vartheta\vartheta}^{(2)}=\frac{1}{2}(F+2B_{2})-2(A_{1}-3A_{2})\cosh(2\varrho)+3(A_{2}-2A_{3})\cosh(4\varrho)$
$a_{\vartheta\vartheta}^{(3)}=-(B_{2}-3B_{3})\cosh(\varrho)-(5A_{2}-10A_{3}-B_{3}+2B_{4})\cosh(3\varrho)+2(3A_{3}-5A_{4})\cosh(5\varrho)$
$a_{\vartheta\vartheta}^{(4)}=-3(B_{3}-2B_{4})\cosh(2\varrho)-(9A_{3}-15A_{4}-3B_{4}+5B_{5})\cosh(4\varrho)+5(2A_{4}-3A_{5})\cosh(6\varrho)$
$a_{\vartheta\vartheta}^{(5)}=-2(3B_{4}-5B_{5})\cosh(3\varrho)-(14A_{4}-21A_{5}-6B_{5}+9B_{6})\cosh(5\varrho)+3(5A_{5}-7A_{6})\cosh(7\varrho)$
$a_{\vartheta\vartheta}^{(6)}=-5(2B_{5}-3B_{6})\cosh(4\varrho)-2(10A_{5}-14A_{6}-5B_{6}+7B_{7})\cosh(6\varrho)+7(3A_{6}-4A_{7})\cosh(8\varrho)$
$a_{\vartheta\vartheta}^{(7)}=-3(5B_{6}-7B_{7})\cosh(5\varrho)-(27A_{6}-36A_{7}-15B_{7}+20B_{8})\cosh(7\varrho)+4(7A_{7}-9A_{8})\cosh(9\varrho)$
$a_{\vartheta\vartheta}^{(8)}=-7(3B_{7}-4B_{8})\cosh(6\varrho)-(35A_{7}-45A_{8}-21B_{8})\cosh(8\varrho)+36A_{8}\cosh(10\varrho)$
$b_{\vartheta\vartheta}^{(1)}=2(E+G)-(G-3A^{\prime}_{1}-B^{\prime}_{2})\cosh(\varrho)+(A^{\prime}_{1}-3A^{\prime}_{2})\cosh(3\varrho)$
$b_{\vartheta\vartheta}^{(2)}=-\frac{1}{2}(3E+2G-2B^{\prime}_{2})-2(A^{\prime}_{1}-3A^{\prime}_{2})\cosh(2\varrho)+3(A^{\prime}_{2}-2A^{\prime}_{3})\cosh(2\varrho)$
$b_{\vartheta\vartheta}^{(3)}=\frac{2}{3}(E+G)-(B^{\prime}_{2}-3B^{\prime}_{3})\cosh(\varrho)-(5A^{\prime}_{2}-10A^{\prime}_{3}-B^{\prime}_{3}+2B^{\prime}_{4})\cosh(3\varrho)+2(3A^{\prime}_{3}-5A^{\prime}_{4})\cosh(5\varrho)$
$b_{\vartheta\vartheta}^{(4)}=-\frac{1}{2}(E+G)-3(B^{\prime}_{3}-2B^{\prime}_{4})\cosh(2\varrho)-(9A^{\prime}_{3}-15A^{\prime}_{4}-3B^{\prime}_{4}+5B^{\prime}_{5})\cosh(4\varrho)+5(2A^{\prime}_{4}-3A^{\prime}_{5})\cosh(6\varrho)$
$b_{\vartheta\vartheta}^{(5)}=\frac{2}{5}(E+G)-2(3B^{\prime}_{4}-5B^{\prime}_{5})\cosh(3\varrho)-(14A^{\prime}_{4}-21A^{\prime}_{5}-6B^{\prime}_{5}+9B^{\prime}_{6})\cosh(5\varrho)+3(5A^{\prime}_{5}-7A^{\prime}_{6})\cosh(7\varrho)$
$b_{\vartheta\vartheta}^{(6)}=-\frac{1}{3}(E+G)-5(2B^{\prime}_{5}-3B^{\prime}_{6})\cosh(4\varrho)-2(10A^{\prime}_{5}-14A^{\prime}_{6}-5B^{\prime}_{6}+7B^{\prime}_{7})\cosh(6\varrho)+7(3A^{\prime}_{6}-4A^{\prime}_{7})\cosh(8\varrho)$
$b_{\vartheta\vartheta}^{(7)}=\frac{2}{7}(E+G)-3(5B^{\prime}_{6}-7B^{\prime}_{7})\cosh(5\varrho)-(27A^{\prime}_{6}-36A^{\prime}_{7}-15B^{\prime}_{7}+20B^{\prime}_{8})\cosh(7\varrho)+4(7A^{\prime}_{7}-9A^{\prime}_{8})\cosh(9\varrho)$
$b_{\vartheta\vartheta}^{(8)}=-\frac{1}{4}(E+G)-7(3B^{\prime}_{7}-4B^{\prime}_{8})\cosh(6\varrho)-(35A^{\prime}_{7}-45A^{\prime}_{8}-21B^{\prime}_{8})\cosh(8\varrho)+36A^{\prime}_{8}\cosh(10\varrho)$

When we introduce a coefficient matrix $C$, we can write simultaneous
equations to determine the vector $\bm{X}$ as

$C\bm{X}=\widetilde{\bm{\Sigma}}_{\rm surf}.$ (13)

This means that we can evaluate $\bm{X}=C^{-1}\widetilde{\bm{\Sigma}}_{\rm
surf}$ from the information of the surfaces of the intruders. Figure 6 shows
the absolute value of each component $|x_{i}|$. This numerically justifies the
truncation of higher order terms in Eq. (11) because the values become
exponentially smaller for larger order. Figure 6 shows that only $x_{1}=E$ and
$x_{21}=B_{1}$ have large absolute values, while the rest are much smaller and
almost zero. This means that $E$ and $B_{1}$ play dominant roles in the stress
fields. We will revisit this in the next section.

Figure 6: (Color online) Magnitude of each component $|x_{i}|$
($i=1,\cdots,35$) for $\varphi=0.80$, $D=10d$, $\Delta=20d$, and $F_{\rm
drag}=1.6\times 10^{-2}k_{n}d$. Here, the open (filled) circles represent
positive (negative) $x_{i}$.

Now, let us check whether the above estimation is valid. Figure 7 shows the
stress fields $\sigma_{xx}$ and $\sigma_{xy}$ estimated from the Airy stress
function and those obtained from the simulations with $\varphi=0.80$, $D=10d$,
$\Delta=20d$, and $F_{\rm drag}=1.6\times 10^{-2}k_{n}d$. This shows that our
estimation qualitatively reproduces the anisotropic behaviors, e.g., the
positive and negative peak angles of $\sigma_{xy}$ in front of the intruders.

Figure 7: (Color online) Spatial distributions of (a) $\sigma_{xx}$ and (b)
$\sigma_{xy}$ obtained from the simulations and from the Airy stress function
(11) with the coefficients in Fig. 6 for $\varphi=0.80$, $D=10d$,
$\Delta=20d$, and $F_{\rm drag}=1.6\times 10^{-2}k_{n}d$. The solid and dashed
lines in panel (a) represent the contours of $5.0\times 10^{-2}$ and
$2.0\times 10^{-3}$, respectively. The solid, dashed, dotted, and dot–dashed
lines in panel (b) correspond to the contours of $3.0\times 10^{-2}$,
$1.0\times 10^{-4}$, $-1.0\times 10^{-4}$, and $-3.0\times 10^{-2}$,
respectively.

##  4 Discussion

Let us discuss our results. First, we consider the magnitude of $M_{\rm Y}$.
As shown in Fig 4, $F_{\rm Y}$ is almost constant ($F_{\rm Y}=9.2\times
10^{2}\mu_{\rm b}mg$). Now, let us estimate this yield force on the bases of a
simple argument. From Fig. 3, it might be a good estimation that the intruders
push the surrounding particles inside a square region whose size is
$\Delta+D=30d$. In this case, the number of particles $N_{\rm c}$ present
inside this region is approximately given by

$N_{\rm
c}\simeq\frac{(\Delta+D)^{2}}{\displaystyle\frac{\pi}{8}(1+1.4^{2})d^{2}}\varphi,$
(14)

where we assume that the numbers of particles of diameters $d$ and $1.4d$ are
the same. Under this assumption, the frictional force from the surrounding
particles becomes

$F_{\rm Y}\simeq\frac{N_{\rm c}}{2}\mu_{\rm
b}(1+1.4^{2})mg\simeq\frac{4}{\pi}(\Delta+D)^{2}\varphi\mu_{\rm b}mg.$ (15)

For our choice of parameters, this approximately becomes $F_{\rm Y}\simeq
9\times 10^{2}\mu_{\rm b}mg$, which gives a good estimation of the simulation
result in Fig. 4.

Next, let us discuss the stress fields. Figure 6 shows that two coefficients
$E$ and $B_{1}$ are dominant over other coefficients. Then, let us assume that

$E=-\frac{F_{\rm drag}}{\pi d},\ B_{1}=\frac{F_{\rm drag}}{\pi d},$ (16)

as suggested in the previous section. In this case, the Airy function becomes

$\chi=\frac{F_{\rm drag}}{\pi}J(1-\vartheta)\cos\vartheta.$ (17)

Let us consider this meaning. The angle $\vartheta$ is related to the angle
$\theta$, which is defined as the angle between the vector $(x,y-\Delta/2)$
and the $x$-axis. The relationship is approximately given by $\theta\sim
1-\vartheta\sim\pi/2-\vartheta$. We also introduce the distance $r$ from the
center of the upper (bottom) intruder to the point $(x,y)$. With the aid of
Eq. (29), this is approximately given by

$r=\sqrt{x^{2}+\left(y\mp\frac{\Delta}{2}\right)^{2}}=\pm\frac{\sqrt{\Delta^{2}-D^{2}}}{2\sinh\varrho}.$
(18)

Substituting Eq. (18) into Eq. (10), we get

$J=\frac{r\sinh\varrho}{\cosh\varrho-\cos\vartheta}.$ (19)

From Eq. (19), Eq. (11) can be approximately estimated as

$\chi\sim\frac{F_{\rm
drag}}{\pi}r\theta\sin\theta\frac{\sinh\varrho}{\cosh\varrho-\cos\vartheta}.$
(20)

In the vicinity of the upper intruder, $\varrho\to\infty$ is satisfied. This
means that $\sinh\varrho/(\cosh\varrho-\cos\vartheta)\to 1$, and then, we get

$\chi\to\frac{F_{\rm drag}}{\pi}r\theta\sin\theta,$ (21)

which is consistent with the case when the concentration force is applied to a
semi-infinite plane Timoshenko .

Figure 8: (Color online) Spatial distribution of $\sigma_{xy}$ when we adopt
the Airy stress function in Eq. (17). The set of parameters is the same as
that in Fig. 7. The solid, dashed, dotted, and dot–dashed lines correspond to
the contours of $3.0\times 10^{-2}$, $1.0\times 10^{-4}$, $-1.0\times
10^{-4}$, and $-3.0\times 10^{-2}$, respectively.

Figure 8 shows the stress field $\sigma_{xy}$ obtained from the Airy stress
function (17) when we choose the same set of parameters as that in Fig. 7.
This figure indicates that the above simple treatment works as the zeroth
order approximation, while the magnitude is underestimated compared with that
obtained from the simulation.

In this study, we do not consider the motions of the intruders in the
$y$-direction for simplicity of discussion. In some previous studies Dhiman20
; Escalante17 ; Cruz16 , forces perpendicular to the drag direction were
investigated, and it was reported that the force between the intruders is
basically attractive, but becomes repulsive owing to particle clogging if the
distance between the intruders is too narrow ($\Delta-D\gtrsim d$). In our
system, let us check the profile of the stress on the surfaces of the
intruders (see Fig. 9), where the stress $\sigma_{yy}$ outside the intruders
is larger than that inside. Integrating the stress around the surface, we can
understand that the force acting on the intruder in the $y$-direction becomes
attractive. This is consistent with the results of previous studies Dhiman20 ;
Escalante17 ; Cruz16 . More details should be studied in the future.

Figure 9: (Color online) Stress field $\sigma_{yy}$ on the surface of the
upper intruder for $\varphi=0.80$, $D=10d$, $\Delta=20d$, and $F_{\rm
drag}=1.6\times 10^{-2}k_{n}d$.

##  5 Summary

We investigated the drag of two cylindrical objects in a two-dimensional
granular environment by the discrete element method. A yield force was found
to exist, and this yield force was insensitive to the distance between the
intruders. When the drag force was above the yield force, we found that the
steady speed of the intruders increased as the distance between intruders
decreased, because particles between them were influenced by both of them. On
the other hand, the stress fields around the intruders were investigated below
the yield force, where the Airy stress function in the bipolar coordinates
well reproduced the stress fields obtained from the simulations.

One of the authors (ST) thanks Hisao Hayakawa for fruitful discussions. This
work is financially supported by Grants-in-Aid of MEXT for Scientific Research
(Grant Nos. JP20K14428 and JP21H01006).

##  Appendix A Determination of System Size

In this Appendix, we explain how the system size ($L_{x}$ and $L_{y}$) is
determined in the simulations. Because we put the smaller particles at $y=\pm
L_{y}/2$, the system size in the $x$-direction must be an integer multiple of
the diameter $d$. We also note that the aspect ratio of the system should be
close to a certain value $\nu$. Under these conditions, let us determine
$L_{x}$ and $L_{y}$ for a given packing fraction $\varphi$.

The areas occupied by the intruders, the bulk particles, and the wall
particles are, respectively, given by

$\displaystyle S_{\rm intruder}$ $\displaystyle=\frac{\pi}{4}D^{2}\cdot N_{\rm
intruder},$ (22a) $\displaystyle S_{\rm bulk}$
$\displaystyle=\frac{1}{2}\frac{\pi}{4}(N-N_{\rm wall})(1+1.4^{2})d^{2},$
(22b) $\displaystyle S_{\rm wall}$
$\displaystyle=\frac{1}{2}\frac{\pi}{4}d^{2}N_{\rm wall},$ (22c)

where $N_{\rm intruder}$ and $N_{\rm wall}$ are the numbers of intruders and
wall particles, respectively. Here, the factor $1/2$ in Eq. (22c) comes from
the fact that half of the particles are located inside the system. Summing up
Eqs. (22a)–(22c), the total area occupied by the intruders, the bulk
particles, and the wall particles is given by

$S=\frac{\pi}{4}N_{\rm
intruder}D^{2}+\frac{\pi}{8}(1+1.4^{2})Nd^{2}-\frac{\pi}{8}1.4^{2}N_{\rm
wall}d^{2}.$ (23)

If the system size is given by $L\times\nu L$, the size $L$ is determined by
Eq. (23) as

$\varphi\nu L^{2}=\frac{\pi}{4}N_{\rm
intruder}D^{2}+\frac{\pi}{8}(1+1.4^{2})Nd^{2}-\frac{\pi}{8}1.4^{2}N_{\rm
wall}d^{2}.$ (24)

Because the wall is made of the smaller particles, $N_{\rm wall}$ should be
given by

$N_{\rm wall}=2\left\lfloor\frac{L}{d}\right\rfloor,$ (25)

where $\lfloor x\rfloor$ is the floor function that returns an integer $n$ for
$x-1<n\leq x$. In this case, the system area is unchanged when we choose

$L_{x}=\lfloor L\rfloor,\quad L_{y}=\frac{\nu L^{2}}{\lfloor L\rfloor}.$ (26)

Rigorously, the aspect ratio is larger than $\nu$, but it is very close to
$\nu$. Substituting Eq. (25) into Eq. (24), we get

$L^{2}+\frac{\pi}{4\varphi\nu}1.4^{2}d^{2}\left\lfloor\frac{L}{d}\right\rfloor-\frac{\pi}{4\varphi\nu}N_{\rm
intruder}D^{2}-\frac{\pi}{8\varphi\nu}(1+1.4^{2})d^{2}N=0.$ (27)

We can numerically determine $L$ from this equation. For our choice of
parameters ($N_{\rm intruder}=2$, $N=4000$, $\varphi=0.80$, and $\nu=3/2$), we
obtain $L/d\simeq 62.65\cdots$, and then $L_{x}=62d$ and $L_{y}=94.96d$.

##  Appendix B Case of three intruders

In this Appendix, we investigate how the drag changes when there exist three
intruders in the system.

First, we discuss the drag law when the drag force is larger than the yield
force. Figure 10 shows the drag for two and three intruders. This suggests
that the yield force is almost unchanged because the areas of interactions
with the intruders are limited to near them. On the other hand, the difference
becomes larger as the steady speed of the intruders increases. This can be
understood as follows. As the speed increases, the number of particles that
interact with the intruders also increases. Then, there exist areas where
particles interact with two adjacent intruders. This area becomes more movable
than other areas. Therefore, the resistance from the surrounding particles
decreases. Here, a more systematic study is needed but is beyond the scope of
this paper.

Figure 10: (Color online) Steady speed of the intruders against the drag force
in the cases of two (open squares) and three intruders (open circles) for
$\varphi=0.80$, $D=10d$, and $\Delta=20d$.

##  Appendix C Derivation of Stress Fields in terms of Bipolar Coordinates

In this Appendix, we revisit the general solutions of the stress fields around
intruders when we adopt bipolar coordinates Jeffery21 . The bipolar
coordinates $(\varrho,\vartheta)$ are transformed from the $(x,y)$–plane via
the following conformal mapping:

$x+iy=ia\coth\left(\frac{\varrho+i\vartheta}{2}\right),$ (28)

with an imaginary number $i$ and a constant $a=\sqrt{\Delta^{2}-D^{2}}/2$.
Here, the explicit relationship between $(\varrho,\vartheta)$ and $(x,y)$ is
given by

$\begin{cases}\displaystyle
x=\frac{a\sin\vartheta}{\cosh\varrho-\cos\vartheta}\\\ \displaystyle
y=\frac{a\sinh\varrho}{\cosh\varrho-\cos\vartheta}\end{cases}.$ (29)

We note that $\varrho={\rm const.}$ and $\vartheta={\rm const.}$ correspond to
the following circles in the $(x,y)$–plane, respectively:

$\displaystyle x^{2}+\left(y-a\coth\varrho\right)^{2}$
$\displaystyle=\frac{a^{2}}{\sinh^{2}\varrho},$ (30a)
$\displaystyle\left(x-a\cot\vartheta\right)^{2}+y^{2}$
$\displaystyle=\frac{a^{2}}{\sin^{2}\vartheta}.$ (30b)

For later convenience, we introduce a quantity $J$ as

$J=\sqrt{\left(\frac{\partial
x}{\partial\varrho}\right)^{2}+\left(\frac{\partial
x}{\partial\vartheta}\right)^{2}}=\frac{a}{\cosh\varrho-\cos\vartheta},$ (31)

which corresponds to the growth rate of the conformal mapping (28).

Now, let us consider the Airy stress function $\chi$. In the $(x,y)$–plane,
the stress function should satisfy

$\Delta\Delta\chi(x,y)=0$ (32)

with the Laplacian $\Delta\equiv(\partial^{2}/\partial
x^{2}+\partial^{2}/\partial y^{2})$. By using the relationship (29), Eq. (32)
in the $(\varrho,\vartheta)$–plane becomes

$\displaystyle\left(\frac{\partial^{4}}{\partial\varrho^{4}}+2\frac{\partial^{4}}{\partial\varrho^{2}\partial\vartheta^{2}}+\frac{\partial^{4}}{\partial\vartheta^{4}}-2\frac{\partial^{2}}{\partial\varrho^{2}}+2\frac{\partial^{2}}{\partial\vartheta^{2}}+1\right)\left(\frac{\chi}{J}\right)=0.$
(33)

The general solution of this equation is given by

$\displaystyle\frac{\chi}{J}$
$\displaystyle=(E\cos\vartheta+F\sin\vartheta+G\cosh\varrho+H\sinh\varrho)\vartheta$
$\displaystyle\hskip 10.00002pt+(B_{0}\cosh\varrho+D_{0}\sinh\varrho)\varrho$
$\displaystyle\hskip
10.00002pt+[A_{1}\cosh(2\varrho)+B_{1}+C_{1}\sinh(2\varrho)+D_{1}\varrho]\cos\vartheta$
$\displaystyle\hskip
10.00002pt+[A_{1}^{\prime}\cosh(2\varrho)+C_{1}^{\prime}\sinh(2\varrho)+D_{1}^{\prime}\varrho]\sin\vartheta$
$\displaystyle\hskip
10.00002pt+\sum_{n=2}^{\infty}\left\\{A_{n}\cosh[(n+1)\varrho]+B_{n}\cosh[(n-1)\varrho]\right.$
$\displaystyle\hskip
20.00003pt\left.+C_{n}\sinh[(n+1)\varrho]+D_{n}\sinh[(n-1)\varrho]\right\\}\cos(n\vartheta)$
$\displaystyle\hskip
10.00002pt+\sum_{n=2}^{\infty}\left\\{A_{n}^{\prime}\cosh[(n+1)\varrho]+B_{n}^{\prime}\cosh[(n-1)\varrho]\right.$
$\displaystyle\hskip
20.00003pt\left.+C_{n}^{\prime}\sinh[(n+1)\varrho]+D_{n}^{\prime}\sinh[(n-1)\varrho]\right\\}\sin(n\vartheta),$
(34)

with constants $E$, $F$, $G$, $H$, $A_{n}$, $B_{n}$, $C_{n}$, $D_{n}$,
$A_{n}^{\prime}$, $B_{n}^{\prime}$, $C_{n}^{\prime}$, and $D_{n}^{\prime}$
($n\geq 0$). Using Eq. (34), the stress can be calculated as Jeffery21

$\displaystyle a\sigma_{\varrho\varrho}$
$\displaystyle=\left[(\cosh\varrho-\cos\vartheta)\frac{\partial^{2}}{\partial\vartheta^{2}}-\sinh\varrho\frac{\partial}{\partial\varrho}\right.$
$\displaystyle\hskip
20.00003pt\left.-\sin\vartheta\frac{\partial}{\partial\vartheta}+\cosh\varrho\right]\left(\frac{\chi}{J}\right)$
(35a) $\displaystyle a\sigma_{\vartheta\vartheta}$
$\displaystyle=\left[(\cosh\varrho-\cos\vartheta)\frac{\partial^{2}}{\partial\varrho^{2}}-\sinh\varrho\frac{\partial}{\partial\varrho}\right.$
$\displaystyle\hskip
20.00003pt\left.-\sin\vartheta\frac{\partial}{\partial\vartheta}+\cos\vartheta\right]\left(\frac{\chi}{J}\right)$
(35b) $\displaystyle a\sigma_{\varrho\vartheta}$
$\displaystyle=-(\cosh\varrho-\cos\vartheta)\frac{\partial^{2}}{\partial\varrho\partial\vartheta}\left(\frac{\chi}{J}\right).$
(35c)

These are the general solutions of the stress fields. Once we get the stresses
$\sigma_{\varrho\varrho}$, $\sigma_{\varrho\vartheta}$, and
$\sigma_{\vartheta\vartheta}$, the stresses on Cartesian coordinates
$\sigma_{xx}$, $\sigma_{xy}$, and $\sigma_{yy}$ can be calculated by

$\displaystyle\sigma_{xx}$
$\displaystyle=\sigma_{\varrho\varrho}\cos^{2}\gamma+\sigma_{\vartheta\vartheta}\sin^{2}\gamma-2\sigma_{\varrho\vartheta}\sin\gamma\cos\gamma,$
(36a) $\displaystyle\sigma_{xy}$
$\displaystyle=(\sigma_{\varrho\varrho}-\sigma_{\vartheta\vartheta})\sin\gamma\cos\gamma+\sigma_{\varrho\vartheta}(\cos^{2}\gamma-\sin^{2}\gamma),$
(36b) $\displaystyle\sigma_{yy}$
$\displaystyle=\sigma_{\varrho\varrho}\sin^{2}\gamma+\sigma_{\vartheta\vartheta}\cos^{2}\gamma+2\sigma_{\varrho\vartheta}\sin\gamma\cos\gamma,$
(36c)

(see, e.g., Ref. 35) where $\gamma$ satisfies Jeffery21

$\displaystyle\cos\gamma$
$\displaystyle=-\frac{\sinh\varrho\sin\vartheta}{\cosh\varrho-\cos\vartheta},$
(37a) $\displaystyle\sin\gamma$
$\displaystyle=-\frac{\cosh\varrho\cos\vartheta-1}{\cosh\varrho-\cos\vartheta}.$
(37b)

## References

  * (1) S. H. Lamb, Hydrodynamics (Dover, New York, 1945).
  * (2) G. K. Batchelor, An Introduction to Fluid Dynamics (Cambridge Univ. Press, Cambridge, 1967).
  * (3) S. Takada and H. Hayakawa, Phys. Rev. Res. 2, 033468 (2020).
  * (4) J. E. Hilton and A. Tordesillas, Phys. Rev. E 88, 062203 (2013).
  * (5) S. Kumar, K. A. Reddy, S. Takada, and H. Hayakawa, arXiv:1712.09057.
  * (6) S. Takada and H. Hayakawa, Granul. Matter 22, 6 (2020).
  * (7) M. Dhiman, S. Kumar, K. A. Reddy, and R. Gupta, J. Fluid Mech. 886, A23 (2020).
  * (8) Y. Takehara, S. Fujimoto, and K. Okumura, EPL 92, 44003 (2010).
  * (9) Y. Takehara and K. Okumura, Phys. Rev. Lett. 112, 148001 (2014).
  * (10) S. Takada and H. Hayakawa, J. Eng. Mech. 143, C4016004 (2017).
  * (11) C. R. Wassgren and J. A. Cordova, Phys. Fluids 15, 3318 (2003).
  * (12) K. A. Reddy, Y. Forterre, and O. Pouliquen, Phys. Rev. Lett. 106, 108301 (2011).
  * (13) F. Guillard, Y. Forterre, and O. Pouliquen, Phys. Fluids 26, 043301 (2014).
  * (14) A. Pal and A. Kudrolli, Phys. Rev. Fluids 6, 124302 (2021).
  * (15) T. Kubota, H. Ishikawa, and S. Takada, EPJ Web Conf. 249, 03033 (2021).
  * (16) X. Cheng, G. Varas, D. Citron, H. M. Jaeger, and S. R. Nagel, Phys. Rev. Lett. 99, 188001 (2007).
  * (17) X. Cheng, L. Gordillo, W. W. Zhang, H. M. Jaeger, and S. R. Nagel, Phys. Rev. E 89, 042201 (2014).
  * (18) T. G. Sano and H. Hayakawa, Phys. Rev. E 86, 041308 (2012).
  * (19) T. G. Sano and H. Hayakawa, Prog. Theor. Exp. Phys. 2013, 103J02 (2013).
  * (20) R. Jewel, A. Panaitescu, and A. Kudrolli, Phys. Rev. Fluids 3, 084303 (2018).
  * (21) R. A. López de la Cruz and G. A. Caballero-Robledo, J. Fluid Mech. 800, 248 (2016).
  * (22) M. F. Acevedo-Escalante and G. A. Caballero-Robledo, EPJ Web Conf. 140, 03048 (2017).
  * (23) T. Hossain and P. Rognon, Phys. Rev. E 102, 022914 (2020).
  * (24) T. Hossain and P. Rognon, Granul, Matter 22, 72 (2020).
  * (25) T. Hossain and P. Rognon, Phys. Rev. Fluids 5, 054306 (2020).
  * (26) R. Bharadwaj and C. Wassgren, Phys. Fluids 18, 043301 (2006).
  * (27) H. Katsuragi, Physics of Soft Impact and Cratering (Springer, Berlin, 2016).
  * (28) M. Espinosa, V. L. Díaz, A. Serrano, and E. Altshuler, arXiv:2110.15311.
  * (29) P. A. Cundall and O. D. L. Strack, Géotechnique 29, 47 (1979).
  * (30) S. Luding, Granul. Matter 10, 235 (2008).
  * (31) H. Kuninaka and H. Hayakawa, J. Phys. Soc. Jpn. 70, 2220 (2001).
  * (32) H. Hayakawa and H. Kuninaka, Chem. Eng. Sci. 57, 239 (2002).
  * (33) V. Garzó, Granular Gaseous Flows —A Kinetic Theory Approach to Granular Gaseous Flows— (Springer, Berlin, 2019).
  * (34) J. Zhang, R. P. Behringer, and I. Goldhirsch, Prog. Theor. Phys. Suppl. 184, 16 (2010).
  * (35) S. P. Timoshenko and J. N. Goodier, Theory of Elasticity (McGraw-Hill, New York, 1970) 3rd ed.
  * (36) K. E. Daniels, J. E. Kollmer, and J. G. Puckett, Rev. Sci. Instrum. 88, 051808 (2017).
  * (37) G. B. Jeffery, Phil. Trans. Roy. Soc. A 221, 265 (1921).