# Dynamics of a droplet in shear flow by smoothed particle hydrodynamics

Kuiliang Wang  Hong Liang  Chong Zhao  Xin Bian  State Key Laboratory of Fluid
Power and Mechatronic Systems, Department of Engineering Mechanics, Zhejiang
University, Hanghzhou 310027, China.  Department of Physics, Hangzhou Dianzi
University, Hangzhou 310018, China.  Hangzhou Shiguangji Intelligient
Electronics Technology Co., Ltd, Hangzhou, 310018, China.

###### Abstract

We employ a multi-phase smoothed particle hydrodynamics (SPH) method to study
droplet dynamics in shear flow. With an extensive range of Reynolds number,
capillary number, wall confinement, and density/viscosity ratio between the
droplet and the matrix fluid, we are able to investigate systematically the
droplet dynamics such as deformation and breakup. We conduct the majority of
the simulations in two dimensions due to economical computations, while
perform a few representative simulations in three dimensions to corroborate
the former. Comparison between current results and those in literature
indicates that the SPH method adopted has an excellent accuracy and is capable
of simulating scenarios with large density or/and viscosity ratios. We
generate slices of phase diagram in five dimensions, scopes of which are
unprecedented. Based on the phase diagram, critical capillary numbers can be
identified on the boundary of different states. As a realistic application, we
perform simulations with actual parameters of water droplet in air flow to
predict the critical conditions of breakup, which is crucial in the context of
atomization.

###### keywords:

droplet; multiphase flow; SPH;

††journal: Elsevier

##  1 Introduction

The deformation and breakup of droplets in shear flow are ubiquitous in
engineering applications. On microfluidic chips, droplets are utilized for
microbial cultivation and material transport [1, 2], and a thorough
understanding of their dynamics in confined flows may improve the efficiency
of production and transportation. In other environmental and industrial
applications such as protection against harmful aerosols, ink-jet printing and
atomization in nozzles [3, 4, 5, 6, 7], liquid droplets are typically in gas
flows. Accordingly, a decent knowledge on their dynamics with a high
density/viscosity ratio against the matrix fluid is significant. To this end,
a comprehensive investigation on the dynamics of a droplet in shear flow,
which involves a wide range of Reynolds number, capillary number, confinements
of the wall, viscosity/density ratio between the two phases, is called for.

Since pioneering works by Taylor on droplet deformation in shear and
extensional flows [8, 9], enormous theoretical and experimental studies have
been conducted. A series of works by the group of Mason [10, 11, 12] further
studied the deformation and burst of droplets, and even depicted the
streamlines inside and around the droplets. Chaffey and Brenner [13] extended
a previous analytical approximation to a second order form, which is crucial
for the non-elliptic deformation of a highly viscous droplet under large shear
rate. Barthes-Biesel and Acrivos [14] expressed the solution of creeping-flow
equations in powers of deformation parameters and applied a linear stability
theory to determine the critical values for the droplet breakup. Hinch and
Acrivos [15] investigated theoretically the stability of a long slender
droplet, which is largely deformed in shear flow. However, early analytical
works rarely considered effects of finite Reynolds number or wall
confinements. In addition, numerous experimental studies have been conducted
on the droplet deformation and breakup [16, 17, 18, 19], where not only the
effects of viscosity ratio between the droplet and the matrix fluid [20, 21],
but also wall confinements [22, 23] have been taken into account.

With advance in computational science, numerical simulation has become a
popular approach to study droplet dynamics in the past decades. Boundary
integral method was among the first to be applied to study deformation of
droplets in stationary and transient states [24], non-Newtonian droplets [25],
and migration of a droplet in shear flow [26]. Moreover, Li et al. [27]
employed a volume-of-fluid (VOF) method and Galerkin projection technique to
simulate the process of droplet breakup. In the work of Amani et al. [28], a
conservative level-set (CLS) method built on a conservative finite-volume
approximation is applied to study the effect of viscosity ratio and wall
confinement on the critical capillary number. In addition, lattice Boltzmann
method (LBM) has been widely employed to study deformation, breakup and
coalescence of droplets [29, 30, 31, 32, 33]; to model viscoelastic droplet
[34] and surfactant-laden droplet [35]. We note that an interface tracing
technique such as VOF, CLS, a phase-field formulation, or immersed boundary
method is often necessary by a flow solver based on Eulerian meshes.

As a Lagrangian method, smoothed particle hydrodynamics (SPH) method has some
advantages in simulating multiphase flows. Since different phases are
identified by different types of particles, the interface automatically
emerges without an auxillary tracing technique, even for a very large
deformation. Moreover, inertia and wall effects can be taken into account
straightforward, in contrast to theoretical analysis or the boundary integral
method. Since its inception in astrophysics, SPH method has been largely
developed and widely applied in various flow problems [36, 37]. Morris [38]
considered the surface tension based on a continuous surface force model and
simulated an oscillating two-dimensional rod in SPH. Hu et al. [39] proposed a
multi-phase model that handles both macroscopic and mesoscopic flows in SPH,
where a droplet in shear flow was selected as a benchmark to validate the
method. Other improvements and modifications have also been proposed for SPH
in the context of multiphase problems [40, 41, 42, 43]. Furthermore, a droplet
or matrix flow with special properties can also be considered. For example,
Moinfar et al. [44] studied the drop deformation under simple shear flow of
Giesekus fluids and Vahabi [45] investigated the effect of thixotropy on
deformation of a droplet under shear flow. Saghatchi et al. [46] studied the
dynamics of a 2D double emulsion in shear flow with electric field based on an
incompressible SPH method. There are also studies on colliding and coalescence
process of droplets by SPH [47, 48]. Simulation of bubbles in liquid is
similar, but can encounter special challenges [49], due to the reverse
density/viscosity ratio as that of droplet in gas.

Previously, simulations of multiphase flows by SPH method often investigated
specific circumstances. Therefore, the objective of this paper is two fold:
firstly, to simulate an extensive range of parameters to examine the SPH
method for multiphase flows; secondly, to fill gaps of unexplored range of
parameters and systematically investigate their influence on the droplet
dynamics. The rest of the paper is arranged as follows: in Sec. 2, we
introduce the multiphase SPH method and a specific surface tension model. We
present validations and extensive numerical results in Sec. 3. We summarize
this work after discussions in Sec. 4.

##  2 Method

###  2.1 Governing equations and surface tension model

We consider isothermal Navier-Stokes equations with a surface tension for
multiphase flow in Lagrangian frame

$\displaystyle\frac{d\rho}{dt}=-\rho\nabla\cdot\mathbf{v},$ (1)
$\displaystyle\frac{d\mathbf{v}}{dt}=\frac{1}{\rho}\left(-\nabla
p+\mathbf{F}_{b}+\mathbf{F}_{v}+\mathbf{F}_{s}\right),$

where $\rho$, $\mathbf{v}$ and $p$ are density, velocity and pressure
respectively. $\mathbf{F}_{b}$ is the body force, which is not considered in
this study. $\mathbf{F}_{v}$, $\mathbf{F}_{s}$ denote viscous force and
surface tension at the interface between two phases, respectively.

Following previous studies of quasi-incompressible flow modeling [38], an
artificial equation of state relating pressure to density can be written as

$p=c_{s}^{2}\left(\rho-\rho_{\mathrm{ref}}\right),$ (2)

where $c_{s}$ is an artificial sound speed and $\rho_{\mathrm{ref}}$ is a
reference density. Theoretically, subtracting the reference density has no
influence on the gradient of pressure, but it can reduce the numerical error
of SPH discretizations for the gradient operator.

For a Newtonian flow, the viscous force $\mathbf{F}_{v}$ simplifies to

$\mathbf{F}_{v}=\mu\nabla^{2}\mathbf{v},$ (3)

where $\mu$ is the dynamic viscosity. We assume surface tension to be uniform
along the interface and do not consider Marangoni force. Therefore, the
surface tension acts on the normal direction of the interface. Moreover, its
magnitude depends on the local curvature as

$\mathbf{F}_{s}=\sigma\kappa\mathbf{\hat{n}}\delta_{s},$ (4)

where $\sigma$, $\kappa$, $\mathbf{\hat{n}}$ are surface tension coefficient,
curvature and unit normal vector to the concave side, respectively;
$\delta_{s}$ is a surface delta function and its discrete form shall be
described later.

To describe the surface tension at the interface between two fluids, a
continuous surface tension model is adopted. As a matter of fact, surface
tension my be written as the divergence of a tensor $\mathbf{T}$ [50, 51]

$\sigma\kappa\mathbf{\hat{n}}\delta_{s}=\nabla\cdot\mathbf{T},$ (5)

where

$\mathbf{T}=\sigma\left(\mathbf{I}-\mathbf{\hat{n}}\otimes\mathbf{\hat{n}}\right)\delta_{s}.$
(6)

To represent a multiphase flow, we define a color function $c$ and set a
unique value for each phase, that is, $c^{\mathrm{I}}=0$ and
$c^{\mathrm{II}}=1$ for the two phases, respectively. Apparently, the color
function has a jump from $0$ to $1$ at the interface between phase
$\mathrm{I}$ and $\mathrm{II}$. Therefore, the unit normal vector can be
represented by the normalized gradient of the color function as

$\mathbf{\hat{n}}=\frac{\nabla c}{\left|\nabla c\right|},$ (7)

and the surface delta function is replaced by the scaled gradient as

$\delta_{s}=\left|\mathbf{n}\right|=\frac{\left|\nabla
c\right|}{\left|c^{\mathrm{I}}-c^{\mathrm{II}}\right|}.$ (8)

###  2.2 SPH method

In SPH, fluid is represented by moving particles carrying flow properties such
as density, velocity and pressure. We largely follow the work of Hu and Adams
[39] and provide a brief derivation here. Density of a particle is calculated
by interpolating the mass of neighboring particles as

$\rho_{i}=m_{i}\sum_{j}W_{ij},$ (9)

where mass $m_{i}$ is constant for every particle. $W_{ij}$ denotes a weight
function for interpolation

$W_{ij}=W\left(\mathbf{r}_{ij},h\right),$ (10)

where $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$ is a relative position
vector from particle $j$ to $i$ and $h$ is the smoothing length. We further
define

$V_{i}=\frac{1}{\sum_{j}W_{ij}},$ (11)

to be an equivalent volume of particle $i$ so that $V_{i}=m_{i}/\rho_{i}$.

The pressure gradient can be computed as

$-\left(\frac{1}{\rho}\nabla
p\right)_{i}=-\sum_{j}\left(V_{i}^{2}p_{i}+V_{j}^{2}p_{j}\right)\frac{\partial
W}{\partial r_{ij}}\mathbf{e}_{ij},$ (12)

where $p_{i}$ and $p_{j}$ are obtained by Eq. (2). The viscous force can be
calculated as

$\left(\mu\nabla^{2}\mathbf{v}\right)_{i}=\sum_{j}\frac{2\mu_{i}\mu_{j}}{\mu_{i}+\mu_{j}}\left(V_{i}^{2}+V_{j}^{2}\right)\frac{\mathbf{v}_{ij}}{r_{ij}}\frac{\partial
W}{\partial r_{ij}},$ (13)

where $\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}$ is the relative velocity
of particle $i$ and $j$ and $r_{ij}=\left|\mathbf{r}_{ij}\right|$ is the
distance between them.

As suggested by Morris [38] and Hu et al. [39], a part of pressure
contribution $\sigma\frac{d-1}{d}\delta_{s}$ is removed to avoid attractive
force and improve the stability of the interactions between SPH particles.
Therefore, we employ

$\mathbf{T}^{\prime}=\sigma\left(\frac{1}{d}\mathbf{I}-\mathbf{\hat{n}}\otimes\mathbf{\hat{n}}\right)\delta_{s}$
(14)

to replace Eq. (6), where $d$ is the spatial dimension. Combining Eq. (8), Eq.
(7) and Eq. (14), we obtain

$\mathbf{T}^{\prime}=\frac{\sigma}{\left|c^{\mathrm{I}}-c^{\mathrm{II}}\right|\left|\nabla
c\right|}\left(\frac{\left|\nabla c\right|^{2}}{d}\mathbf{I}-\nabla
c\otimes\nabla c\right).$ (15)

The gradient of color function between phase $\mathrm{I}$ and phase
$\mathrm{II}$ can be calculated in SPH as

$\nabla
c_{i}=\frac{1}{V_{i}}\sum_{j}V_{j}^{2}\left(c_{j}-c_{i}\right)\frac{\partial
W}{\partial r_{ij}}\mathbf{e}_{ij},$ (16)

where $c_{i}$ (or $c_{j}$) is initially assigned to be $c^{\mathrm{I}}$ or
$c^{\mathrm{II}}$ according to which phase particle $i$ (or $j$) consititutes.
Substitute Eq. (16) into Eq. (15) to obtain stress tensor

$\mathbf{T}^{\prime}_{i}=\frac{\sigma}{\left|\nabla
c_{i}\right|\left|c^{\mathrm{I}}-c^{\mathrm{II}}\right|}\left(\frac{\left|\nabla
c_{i}\right|^{2}}{d}\mathbf{I}-\nabla c_{i}\otimes\nabla c_{i}\right).$ (17)

Finally, the surface force term is calculated by the stress tensor using the
SPH expression for divergence

$\left(\sigma\kappa\mathbf{\hat{n}}\delta_{s}\right)_{i}=\sum_{j}\frac{\partial
W}{\partial
r_{ij}}\mathbf{e}_{ij}\cdot\left(V_{i}^{2}\mathbf{T}^{\prime}_{i}+V_{j}^{2}\mathbf{T}^{\prime}_{j}\right).$
(18)

It is simple to see that the discrete version of $\delta_{s}$ in SPH is

$\left(\delta_{s}\right)_{i}=\frac{1}{V_{i}\left|c^{\mathrm{I}}-c^{\mathrm{II}}\right|}\left|\sum_{j}V_{j}^{2}\left(c_{j}-c_{i}\right)\frac{\partial
W}{\partial r_{ij}}\mathbf{e}_{ij}\right|,$ (19)

which has a finite support to remove the singularity and distributes the
surface tension onto a thin layer of two fluids across the interface.

###  2.3 Computational settings

The quintic kernel is adopted as weight function

$W=\phi\begin{cases}(3-R)^{5}-6(2-R)^{5}+15(1-R)^{5}&0\leq R<1;\\\
(3-R)^{5}-6(2-R)^{5}&1\leq R<2;\\\ (3-R)^{5}&2\leq R<3;\\\ 0&R\geq
3,\end{cases}$ (20)

where $R=r/h$ and $h$ is the smoothing length. $\phi$ is a normalization
coefficient which equals $1/120$, $7/(478\pi)$ and $1/(120\pi)$ in one, two
and three dimensions, respectively. We set $h=1.2\Delta x$ with $\Delta x$ as
the initial spacing distance between particles. This means that the support
domain of the kernel function is truncated at $3.6\Delta x$, namely the cutoff
$r_{c}=3.6\Delta x$. According to our tests, a smoothing length of $1.2\Delta
x$ is almost optimal for an excellent accuracy while avoiding the pairing
instability. A detailed discussion on this issue is referred to Price [52].

Since we adopt a weakly compressible formulation, the sound speed $c_{s}$
should be large enough to restrict the density fluctuations. Based on a scale
analysis, Morris et al. [38, 53] suggested that $c^{2}_{s}$ should be
comparable to the largest of

$\frac{U^{2}}{\Delta},\ \frac{\mu U}{\rho_{0}L\Delta},\ \frac{FL}{\Delta},\
\frac{\sigma\kappa}{\rho_{0}\Delta},$ (21)

where $\Delta$ is the density variation and $U$, $L$, $F$, $\kappa$ and
$\sigma$ are typical velocity, length, body force, curvature and surface
tension coefficient, respectively. Accordingly, for multiphase flows the sound
speed may be different for each phase. In all simulations, we set identical
$\Delta\leq 0.5\%$ for each phase and calculate $c_{s}$ accordingly.

At every time step, the minimal relative density is recorded among all
particles, that is,

$\rho_{min}=min\left\\{min\left\\{\frac{\rho_{i}}{\rho_{0}^{\mathrm{I}}}\right\\},min\left\\{\frac{\rho_{j}}{\rho_{0}^{\mathrm{II}}}\right\\}\right\\},$
(22)

where particle $i$ belongs to phase $\mathrm{I}$ and particle $j$ belongs to
phase $\mathrm{II}$; $\rho_{0}^{\mathrm{I}}$, $\rho_{0}^{\mathrm{II}}$ are
initial densities for the two phases, respectively. Thereafter,
$\rho_{\mathrm{ref}}^{\mathrm{I}}=0.99\rho_{min}\rho_{0}^{\mathrm{I}}$,
$\rho_{\mathrm{ref}}^{\mathrm{II}}=0.99\rho_{min}\rho_{0}^{\mathrm{II}}$ are
subtracted as reference density for each phase in Eq. (2) to compute the
particle pressure. This operation is performed to reduce numerical errors in
calculating the pressure gradient while still keeping repulsive forces between
particles.

The explicit velocity-Verlet method is adopted for time integration and a time
step is chosen appropriately for stability [38].

##  3 Numerical Results

We consider a shear flow generated by two parallel walls with opposite
velocity of magnitude $U$. Periodic boundaries apply in the $x$ direction. The
computational domain is with length $L$ and height $H$. A circular droplet
with radius $R_{0}$ is initially located at the center of the computational
domain, as shown in Fig. 1.

Figure 1: Schematic representation of a droplet with initial radius $R_{0}$ in
a matrix fluid between two parallel walls with distance $H$. The blue dashed
lines represent periodic boundaries with a distance $L$. The continuous phase
has viscosity $\mu_{c}$ while the dispersed phase has viscosity
$\mu_{d}=\lambda\mu_{c}$.

No-slip boundary condition is applied at the wall-fluid interfaces using the
method proposed by Morris [53].

Five dimensionless parameters that determine the deformation of the droplet
are Reynolds number $Re=\rho_{c}\dot{\gamma}R_{0}^{2}/\mu_{c}$, Capillary
number $Ca=\dot{\gamma}R_{0}\mu_{c}/\sigma$, confinement ratio $R_{0}/H$,
viscosity ratio $\lambda=\mu_{d}/\mu_{c}$ and density ratio
$\alpha=\rho_{d}/\rho_{c}$, where $\dot{\gamma}=2U/H$ is the shear rate,
$\sigma$ is the surface tension coefficient, $\rho_{d}$ and $\mu_{d}$ are
density and viscosity of the dispersed fluid phase inside the droplet while
and $\rho_{c}$ and $\mu_{c}$ are for the continuous fluid phase, respectively.

In Sec. 3.1, we study the deformation for an intact droplet while considering
the effects due to the five dimensionless numbers. In Sec. 3.2, we examine the
breakup of the droplet. In Sec. 3.3, we summarize the droplet dynamics for
both intact shape and breakup in phase diagrams. In Sec. 3.4, we demonstrate
the deformation and breakup with physical parameters of a water droplet in air
flow as an industrial application.

###  3.1 Droplet deformation

Figure 2: Trace of interface and parameters for the measurement of droplet
deformation.

When the shear is mild, the droplet remains intact and deforms to arrive at a
stable shape eventually. The degree of droplet deformation can be quantified
by the Taylor deformation parameter $D=(A-B)/(A+B)$, where $A$ is the greatest
length and $B$ is the breadth of the droplet as shown in Fig. 2. To validate
our method, we first compare our results of transient deformations with that
of Sheth and Pozrikidis using immersed boundary method within the finite
difference method [54]. We follow their work to set $L=H=4R_{0}=1$,
$\rho_{d}=\rho_{c}=1$, $\mu_{d}=\mu_{c}=0.5$ and adjust shear rate and surface
tension. The two walls slide with velocities $\pm\frac{1}{2}\dot{\gamma}H$ to
generate a clockwise rotation of the droplet. Two resolutions are considered
for particles initially placed on squared lattice: $\Delta x=2R_{0}/25$ and
$R_{0}/25$, corresponding to the droplet containing $N=484$ and $1976$
particles, respectively.

(a) $\dot{\gamma}t=0$

(b) $\dot{\gamma}t=4$

(c) $\dot{\gamma}t=20$

Figure 3: Particle distribution for $H=4R_{0}$, $\alpha=\lambda=1$,
$Re=0.125$, $Ca=0.45$, $L=4R_{0}$ at (a) initial configuration. (b) maximum
elongation. (c) steady state. (d) Taylor deformation parameter as a function
of time.

We present particle distributions and $D$ as functions of time in Fig. 3 for a
typical simulation with $Re=0.125$ and $Ca=0.45$. We note that the deformation
of the droplet may oscillate in time and its maximum elongation does not
necessarily takes place at the steady state of a very long time.

(a) $Re=0.125$

(b) $Re=1.25$

(c) $Re=6.25$

(d) $Re=12.5$

Figure 4: Transient deformations of a droplet compared with results of Sheth
and Pozrikidis [54] with $\alpha=\lambda=1$, $H=4R_{0}$, $L=4R_{0}$ and (a)
$Re=0.25$, $Ca=0.1,\,0.45$. (b) $Re=2.5$, $Ca=0.1,\,0.4$. (c) $Re=12.5$,
$Ca=0.1,\,0.2$. (d) $Re=25$, $Ca=0.025,\,0.2$

We further focus on the transient deformations in short time in Fig. 4 so that
we can compare our results with those of Sheth and Pozrikidis [54]. It can be
readily seen that our results with low resolution $\Delta x=0.02$ or $N=484$
already reproduce the reference very well for different Reynolds numbers
and/or capillary numbers. As the reference is within a rather short time
period, some interesting phenomenon such as oscillation of the Taylor
deformation parameter $D$ is not captured, as indicated for $Re=12.5$ and
$Ca=0.025$ on Fig. 4(d).

Figure 5: Effects of $Re$ and $Ca$ on 2D droplet deformation when
$\alpha=\lambda=1$, $H=8R_{0}$, $L=8R_{0}$. (a) Taylor deformation parameters.
(b) Droplet orientation. (c) Steady shapes under different $Ca$ when $Re=0.1$.
The droplet already breaks up at $Ca\gtrsim 0.4$ for $Re=1$, results of which
are omitted here and presented in a later section.

To validate our method for vanishing Reynolds numbers, we calculate the
stationary deformation and orientation of the droplet with respect to $Ca$. We
follow Zhou and Pozrikidis [55] to set $L=H=8R_{0}=2$, $\rho_{d}=\rho_{c}=1$,
$\mu_{d}=\mu_{c}=0.5$ and adjust shear rate and surface tension accordingly.
The deformation parameter $D$ and orientation $\theta$ (defined on Fig. 2) as
functions of $Ca$ (up to $Ca=1$) for $Re=0.01$ are shown in Fig. 5. Results
for $Re=0.1$ and $1$ are also given for comparison, where droplet breakup
already takes places at $Ca\gtrsim 0.4$ for $Re=1$. The difference between the
results of $Re=0.1$ and $Re=0.01$ is insignificant and they both resemble the
results of boundary integral method for Stokes flow [55]. We can readily
conclude that $Re=0.1$ is small enough to approximate the Stokes flow and
present the steady shapes accordingly on Fig. 5.

Figure 6:  A typical evolution of deformation for an initially circular 2D
droplet in shear flow: $\alpha=\lambda=1$, $Re=0.1$, $Ca=0.4$, $H=16R_{0}$,
$L=16R_{0}$. (a) Droplet deformation over time. (b) Streamlines: the color
represents the magnitude of velocity and a red line indicates the droplet
interface.

We further present the contours and streamlines for a typical evolution of
droplet deformation at vanishing Reynolds number in Fig. 6.

Figure 7: Effects of confinement ratio and $Ca$ on 2D droplet deformation when
$\alpha=\lambda=1$, $Re=0.1$, $L=16R_{0}$. (a) Taylor deformation parameters.
(b) Droplet orientation.

We commence to investigate the effects of confinement and set $L=16R_{0}$ to
minimize the periodic artifacts. We first restrict out attention to $Re=0.1$,
$\alpha=1$ and $\lambda=1$. Four ratios of confinement are considered:
$H=2.4R_{0}$, $4R_{0}$, $8R_{0}$ and $16R_{0}$. The deformation parameter as a
function of $Ca$ is shown in Fig. 7. As we can see, a smaller distance of the
two walls enhances the elongation of droplet and makes its long axis align
more horizontally. As we relax the confinement, the relation between $D$ and
$Ca$ becomes linear and the difference between $H=16R_{0}$ and $H=8R_{0}$ is
already negligible.

Figure 8: Effects of viscosity ratio $\lambda$ and $Ca$ on 2D droplet
deformation when $\alpha=1$, $Re=0.1$, $H=16R_{0}$, $L=16R_{0}$. (a) Taylor
deformation parameters. (b) Droplet inclination.

Furthermore, we simulate cases where the droplet and the matrix flow are two
fluids with different physical properties. We first consider two fluids of the
same density but with different viscosities. We choose a computational domain
of $16R_{0}\times 16R_{0}$ and set $Re=0.1$, $\alpha=1$ and $\lambda$ ranges
from 0.1 to 10. Initial space $\Delta x$ among nearest particles is
$2R_{0}/25$ so a droplet contains 484 particles. The deformation parameter as
a function of $Ca$ is shown in Fig. 8. As we can see, the deformation
increases as $\lambda$ increases from $0.1$ to $10$. In this range of
$\lambda$, a droplet with lower viscosity has a smooth inside circulation and
fast reaction which can reduce the elongation [16, 20].

Figure 9: Effects of density ratio $\alpha$ and $Ca$ on 2D droplet deformation
when $\lambda=1$, $Re=0.1$, $H=16R_{0}$, $L=16R_{0}$. (a) Taylor deformation
parameters. (b) Droplet inclination.

The other case is that fluids inside and outside the droplet have the same
viscosity but different densities. The sound speed is chosen according to the
ratio of initial density to balance the pressure

$\frac{c_{s}^{c}}{c_{s}^{d}}=\frac{\rho_{ref}^{d}}{\rho_{ref}^{c}}=\frac{\rho_{d}}{\rho_{c}}=\alpha,$
(23)

where $c_{s}^{c}$, ${c_{s}^{d}}$ and $\rho_{ref}^{c}$, $\rho_{ref}^{d}$ are
sound speeds and reference densities used for fluids outside and inside the
droplet. As shown in Fig. 9, the difference between deformations of droplet
under density ratio $0.1-10$ is very small except obvious lower inclination at
small $Ca$ when $\alpha=0.1$. In this small Reynolds number regime ($Re=0.1$),
the density ratio has negligible influence and only the capillary number
determines the droplet deformation.

Figure 10: Effects of box length $L$ and width $W$ on the Talor deformation
parameter $D$ of 3D droplet deformation when $\alpha=\lambda=1$, $Re=0.1$,
$Ca=0.2$, $H=4R_{0}$, compared to the analytical prediction of Shapira and
Haber (1990)

Figure 11:  Effects of confinement ratio and $Ca$ on 3D droplet deformation in
shear flow when $\alpha=1$ and $Re=0.1$, $L=24R_{0}$, $W=8R_{0}$, compared
with predictions of Shapira and Haber [56] and experiment data of Sibillo et
al. [57]

.

In 3D simulations, the width of simulation box $W$ is an additional
computational parameter compared to the 2D simulations. To compare with
analytical predictions or experiment data, the length and width of simulation
box are numerical and should be large enough. One set of parameters of
$Re=0.1$, $Ca=0.2$, $\alpha=\lambda=1$, $H=4R_{0}$ are selected and different
length $L$ and width $W$ of simulation box are examined. According to our
simulations, the deformation basically decreases with the increase of $L$
and/or width $W$. We compare the Taylor deformation parameter $D$ in steady
states of our simulations with the analytical prediction of Shapira and Haber
[56]. The differences between our results and analytical prediction under
different $L$ and $W$ are plotted in Fig. 10. It can be seen that when $L$ is
larger than $24R_{0}$ and $W$ is larger than $8R_{0}$, the results has little
change with the increase of $L$ and/or $W$. Fig. 11 shows the steady
deformation of 3D droplets in shear flow when $L=24R_{0}$, $W=8R_{0}$,
$Re=0.1$ and $\alpha=\lambda=1$ with different $Ca$ and confinement in $H$
direction, compared with theoretical predictions of Shapira and Haber [56] and
experiment data of Sibillo et al. [57]. Our results agree well with both
anlaytical and experiment references at $Ca=0.1$ and $0.2$, whereas are closer
to the experimental data at $Ca=0.3$. The deformation increases with the
confinement ratio $R_{0}/H$, which has the same trend as for 2D cases.

###  3.2 Droplet breakup

Figure 12:  Breakup type A: evolution of an initially circular 2D droplet
breakup in shear flow with $\alpha=1$, $\lambda=0.2$, $Re=0.1$, $Ca=10$,
$H=16R_{0}$, $L=16R_{0}$. (a) Droplet shape. (b) Streamlines: The color
represents the magnitude of velocity and red outlines in the background
represent the droplet interface. Figure 13:  Breakup type B: evolution of an
initially circular 2D droplet breakup in shear flow with $\alpha=\lambda=1$,
$Re=0.1$, $Ca=0.9$, $H=16R_{0}$, $L=16R_{0}$. (a) Droplet shape. (b)
Streamlines: The color represents the magnitude of velocity and red outlines
in the background represent the droplet interface.

When the shear is strong, the droplet is over-stretched to break up. We find
two patterns of breakup process under different viscosity ratios in
simulations. As shown in Fig. 12, when $\alpha=1$, $\lambda=0.2$, $Re=0.1$,
$Ca=10$, and $L=H=16R_{0}$, a droplet is rotated and then stripped of its main
body near the surface and gradually breaks apart. We call this breakup type A.
This type is also found in the experiment study of Grace and they call it ”tip
streaming breakup” [20]. The conditions for type A breakup happening is
exhibited in the next section. Fig. 13 shows another set of typical snapshots
of the droplet shape and flow fields in shear flow breaks when
$\alpha=\lambda=1$, $Re=0.1$, $Ca=0.9$ and $L=H=16R_{0}$. In this simulation,
a droplet is stretched and its waist becomes slender and slender and finally
breaks up. We call this breakup type B.

Figure 14:  Evolution of an initially spherical 3D droplet breakup in shear
flow when $\alpha=\lambda=1$, $Re=0.1$, $Ca=0.46$, $H=2.857R_{0}$,
$L=16R_{0}$: particle distribution (left) and interface (right).

To encompass the breakup of a 3D droplet with a large elongation, we employ a
rather long computational domain with $L=32R_{0}$. Fig. 14 shows the dynamics
of the breakup with $Re=0.1$, $H=2.857R_{0}$, $Ca=0.46$ and
$\alpha=\lambda=1$. Left side are SPH particle distributions and right side
are corresponding contour interfaces processed by SPH kernel interpolation
into mesh cells. The color represents the magnitude of velocity. We adopt the
same $Ca$ and $R_{0}/H$ as the experiment in creeping flow by Sibillo et al.
[57]. The shape of the droplet in the breakup process of our simulation is
very close to their experimental observation. Only a slight difference appears
in the final stage: in the experiment, the droplet is divided into three main
parts, while in our simulation the middle part continues to split into two
smaller droplets. In contrast to the 2D case, a 3D droplet has a more slender
shape before breaking up.

###  3.3 Phase diagram

To clearly visualize the states of a droplet in different conditions, we
consider a range of Reynolds numbers, capillary numbers, and
confinements/density/viscosity ratios and summarize our simulation results
into phase diagrams. Thereafter, we may estimate the critical capillary number
$Ca_{c}$ that segments the intact and breakup states and further investigate
how it is influenced by other dimensionless parameters.

Figure 15: Phase diagram of 2D droplets states under different confinement,
$Re$ and $Ca$ when $\alpha=\lambda=1$, $L=16R_{0}$.

(a) $H=2.4R_{0}$

(b) $H=4R_{0}$

(c) $H=8R_{0}$

(d) $H=16R_{0}$

Figure 16: Phase diagram of 2D droplets states under different confinement,
$Re$ and $Ca$ when $\alpha=\lambda=1$, $L=16R_{0}$. (a) $H=2.4R_{0}$.(b)
$H=4R_{0}$. (c) $H=8R_{0}$. (d) $H=16R_{0}$

For $\lambda=\alpha=1$, we perform a group of 2D simulations with different
Reynolds number $Re=0.01,0.1,1,10$ and confinement
$H=2.4R_{0},4R_{0},8R_{0},16R_{0}$ with $Ca\in[0.1,1.1]$ and $L=16R_{0}$. For
a general overview, the states of the droplet are summarized in Fig. 15. To
get a clear view, we slice the phase diagram by two perspectives. Firstly, we
divide results into groups of the same confinement to reveal the influence of
$Re$ on $Ca_{c}$ as shown in Fig. 16. Overall it is apparent that a higher
$Re$ reduces $Ca_{c}$. Three scenarios are special: under confinement
$H=2.4R_{0}$, $4R_{0}$ and $8R_{0}$, we can not differentiate $Ca_{c}$ between
$Re=0.01$ and $0.1$.

(a) $Re=0.01$

(b) $Re=0.1$

(c) $Re=1$

(d) $Re=10$

Figure 17: Phase diagram of 2D droplets states under different confinement,
$Re$ and $Ca$ when $\alpha=\lambda=1$, $L=16R_{0}$. (a) $Re=0.01$. (b)
$Re=0.1$. (c) $Re=1$. (c) $Re=10$

From another perspective of $Ca$ versus confinement ratio for each $Re$ on
Fig. 17, we are not able to find a universal pattern. Under $Re=0.01$,
$Ca_{c}$ decreases with $R_{0}/H$ while under $Re=10$, $Ca_{c}$ increases with
$R_{0}/H$. Whereas, under $Re=0.1$ and $1$, $Ca_{c}$ has no monotonic relation
with $R_{0}/H$.

Furthermore, we investigate effects of viscosity ratio
$\lambda=\mu_{d}/\mu_{c}\in[0.1,10]$ on the droplet dynamics for $Re=0.1$ and
three confinement ratios $H=4R_{0}$, $8R_{0}$ and $16R_{0}$. The results are
shown in Fig. 18. For breakup type A, the droplet rotates and is stripped off
as described in Sec. 3.2; Breakup type B represents that a droplet is
stretched and breaks up in the middle. Under $Re=0.1$, type A is observed only
if the droplet has a much smaller viscosity compare to the matrix fluid.
Overall, $Ca_{c}$ decreases with the increase of $\lambda$. However, we notice
a flatten trend or even a reverse trend with small difference for $Ca_{c}$
from $\lambda=5$ to $\lambda=10$, as shown on the insets of Fig. 18. According
to the study of Karam et al. and Grace [16, 20], a maximum transfer of energy
takes place across an interface, which demands this trend.

(a) $H=4R_{0}$

(b) $H=8R_{0}$

(c) $H=16R_{0}$

Figure 18: Phase diagram of 2D droplets states under different confinement,
viscosity ratios $\lambda$ and $Ca$ when $\alpha=1$, $Re=0.1$, $L=16R_{0}$.
(a) $H=4R_{0}$. (b) $H=8R_{0}$. (c) $H=16R_{0}$. (d) two patterns of breakup
Figure 19:  Phase diagram of 3D droplets states under different $Re$ and $Ca$
when $\alpha=\lambda=1$, $H=4R_{0}$, $L=32R_{0}$, $W=8R_{0}$.

Due to highly computational cost in 3D, we only consider a moderate
confinement $H/R_{0}=4$ and perform a group of simulations to draw a phase
diagram in the plane of $Ca$ and $Re$, as shown in Fig. 19. The size of the
simulation box is $L=32R_{0}$, $W=8R_{0}$, $H=4R_{0}$. As in 2D case, the
critical $Ca_{c}$ decreases with increasing $Re$ in 3D, as shown in in Fig.
16(b). However, the critical capillary number $Ca_{c}$ in 3D case is
significantly smaller than that of 2D case.

###  3.4 Water droplet in air flow

As one specific application, we employ our method to predict the breakup of a
water droplet in shear flow of air. The critical capillary number or the shear
rate determined is helpful to design an effective atomization device. Actual
physical properties of water and air around 20∘C are adopted:
$\rho_{d}=998.2\;kg\cdot m^{-3}$, $\mu_{d}=1.0087\times 10^{-3}\;Pa\cdot s$
and $\rho_{c}=1.205\;kg\cdot m^{-3}$, $\mu_{c}=1.81\times 10^{-5}\;Pa\cdot s$
are set for water (dispersed) phase and air (continuous) phase, respectively;
surface tension coefficient $\sigma=72.75\times 10^{-3}\;N\cdot m^{-1}$ is set
for the water-air interface.

Figure 20: Diagram for states of water droplets in air shear flow under
different $Re$ and $Ca$ when $H=16R_{0}$, $L=16R_{0}$.

We perform a relative large range of Reynolds numbers and depict a phase
diagram on the plane of $Re$ and $Ca$ in logarithmic-logarithmic scales on
Fig. 20. This allows us to connect the results with the same droplet size and
observe its behavior while changing $Re$ and $Ca$. Points on each dotted line
represent the droplet of the same radius, as marked in the figure. For
example, we have a line of dynamics for the droplet with $R_{0}=10\mu m$ under
shear rates of $1\times 10^{6}s^{-1}$, $2\times 10^{6}s^{-1}$, $5\times
10^{6}s^{-1}$, $1\times 10^{7}s^{-1}$, $2\times 10^{7}s^{-1}$; another line of
dynamics for the droplet with $R_{0}=100\mu m$ under shear rates of $5\times
10^{4}s^{-1}$, $1\times 10^{5}s^{-1}$, $2\times 10^{5}s^{-1}$, $5\times
10^{5}s^{-1}$, $1\times 10^{6}s^{-1}$. Furthermore, we observe that if the
$Re$ is on the order of $100$, the critical $Ca$ for breakup is very sensitive
to $Re$. We also perform a group of 3D simulations for a droplet with
$R_{0}=50\mu m$ under shear rates of $1\times 10^{5}s^{-1}$, $2\times
10^{5}s^{-1}$, $5\times 10^{5}s^{-1}$, $1\times 10^{6}s^{-1}$, $2\times
10^{6}s^{-1}$. The 3D results for the critical point of breakup is close to
that of the 2D results.

##  4 Conclusions and disucssions

In this study, we employed a multi-phase SPH method to simulate droplet
deformation and breakup subjected to a simple shear flow in an extensive range
of physical parameters. We performed both 2D and 3D simulations and validated
them by benchmarks: transient deformations and steady shapes of droplets are
compared with previous simulations, analytical derivations and experimental
data. These results indicate that the method is reliable to simulate droplet
dynamics in general. We wish to emphasize the convenience of SPH method in
simulating multi-phase problems, as we can leverage on its Lagrangian nature
and differentiate different phases by particle species. In addition, the
algorithm and data structure for 2D and 3D simulations have tiny difference
and therefore, it is a simple task to extend the code from 2D to 3D.
Economical 2D simulations allow us to investigate a wide range of physical
parameters in five dimensions, which serve as a guide to 3D realistic
situations. From the results, we come to the following conclusions.

(1) A larger Reynolds number $Re$ or capillary number $Ca$ leads to a more
considerable deformation of the droplet. The transient and steady-state
deformations of the droplet in our study are in good agreement with the
previous studies but beyond their time limits [54, 55].

(2) Under low Reynolds number ($Re=0.1$), a stronger confinement due to the
walls enhances the steady-state deformation in both 2D and 3D simulations.
When the walls are separated further apart, the Taylor deformation parameter
is almost linear with respect to $Ca$. The influence of confinement on the
deformation of a droplet has been studied by Shapira and Haber by a first-
order analytical solution based on Lorentz’s reflection method. They proved
that the walls do not influence the shape of deformed droplet but increases
the deformation magnitude with a term of order $\left(R_{0}/H\right)^{3}$
[56]. The experiment data of Sibillo et al. illustrate satisfactory agreement
with the predictions of Shapira and Haber except for the droplet being within
a small gap, where the reflection analysis is expected to fail [57]. Our 3D
simulation results resemble the whole set of experiment data even when the
droplet is within the small gap, which suggests the method as an applicative
tool for more realistic situations in microfluidics.

(3) The effects of wall confinement on the critical capillary number $Ca_{c}$
are not universal under different $Re$. When $Re=0.1$, a closer gap of walls
reduces $Ca_{c}$. This is because a closer gap of walls increases the
deformation as described above. But when $Re$ is larger, the relation between
$Ca_{c}$ and the confinement ratio is unclear. From our observation, this non-
monotonic relation results from an interplay of influences by the shear
strength and the stability of the whole flow field. On the one hand, the shear
stress transferred to the droplet from the wall is more pronounced in stronger
confinement [56], thus closer walls reduce the $Ca_{c}$. On the other hand,
the narrower channel reduces the instability of the flow and restricts droplet
movements, thus increases the $Ca_{c}$.

(4) Under $Re=0.1$ and the range of viscosity ratio
$\lambda\in\left[0.1,1\right]$, a higher $\lambda$ causes a larger
deformation. The effect of $\lambda$ on $Ca_{c}$ is not monotonic when
$\lambda>1$ and there is a minimum value of $Ca_{c}$ between $\lambda=1$ and
$\lambda=10$. The existence of a minimal $Ca_{c}$ among different $\lambda$
has also been found by previous experiment studies [16, 20], when $\lambda$ is
about $1$. The discrepancy between our results and the previous ones are
attributed to the difference between 2D and 3D cases. At the same $Re$, the
influence of density ratio on droplet deformation is much smaller compared
with that of the viscosity ratio.

(5) As an application, a phase diagram obtained by actual physical parameters
of water and air is depicted to predict the magnitude of shear rate for
breaking a droplet of certain size, which is helpful in the designing
atomization nozzles.

## Acknowledgements

K. Wang and X. Bian acknowledge the national natural science foundation of
China under grant number: 12172330. This work is partially supported by
Hangzhou Shiguangji Intelligent Electronics Technology Co., Ltd, Hangzhou,
China.

## Acknowledgements

K. Wang and X. Bian acknowledge the national natural science foundation of
China under grant number: 12172330. This work is partially supported by
Hangzhou Shiguangji Intelligient Electronics Technology Co., Ltd, Hangzhou,
China.

## References

  * [1] S. L. Anna, Droplets and bubbles in microfluidic devices, Annual Review of Fluid Mechanics 48 (2016) 285–309.
  * [2] O. J. Dressler, X. Casadevall i Solvas, A. J. DeMello, Chemical and biological dynamics using droplet-based microfluidics, Annual Review of Analytical Chemistry 10 (2017) 1–24.
  * [3] J. Liu, M. Hao, S. Chen, Y. Yang, J. Li, Q. Mei, X. Bian, K. Liu, Numerical evaluation of face masks for prevention of covid-19 airborne transmission, Environmental Science and Pollution Research (2022) 1–15.
  * [4] D. Lohse, Fundamental fluid dynamics challenges in inkjet printing, Annual review of fluid mechanics 54 (2022) 349–382.
  * [5] O. Aydin, R. Unal, Experimental and numerical modeling of the gas atomization nozzle for gas flow behavior, Computers & Fluids 42 (1) (2011) 37–43.
  * [6] C.-r. Si, X.-j. Zhang, J.-b. Wang, Y.-j. Li, Design and evaluation of a laval-type supersonic atomizer for low-pressure gas atomization of molten metals, International Journal of Minerals, Metallurgy, and Materials 21 (6) (2014) 627–635.
  * [7] Z. Xu, T. Wang, Z. Che, Droplet deformation and breakup in shear flow of air, Physics of Fluids 32 (5) (2020) 052109.
  * [8] G. I. Taylor, The viscosity of a fluid containing small drops of another fluid, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 138 (834) (1932) 41–48.
  * [9] G. I. Taylor, The formation of emulsions in definable fields of flow, Proceedings of the Royal Society of London. Series A, containing papers of a mathematical and physical character 146 (858) (1934) 501–523.
  * [10] W. Bartok, S. Mason, Particle motions in sheared suspensions: Viii. singlets and doublets of fluid spheres, Journal of Colloid Science 14 (1) (1959) 13–26.
  * [11] F.-D. Rumscheidt, S. Mason, Particle motions in sheared suspensions xii. deformation and burst of fluid drops in shear and hyperbolic flow, Journal of Colloid Science 16 (3) (1961) 238–261.
  * [12] S. Torza, C. Henry, R. Cox, S. Mason, Particle motions in sheared suspensions. xxvi. streamlines in and around liquid drops, Journal of Colloid and Interface Science 35 (4) (1971) 529–543.
  * [13] C. E. Chaffey, H. Brenner, A second-order theory for shear deformation of drops, Journal of Colloid and Interface Science 24 (2) (1967) 258–269.
  * [14] D. Barthes-Biesel, A. Acrivos, Deformation and burst of a liquid droplet freely suspended in a linear shear field, Journal of Fluid Mechanics 61 (1) (1973) 1–22.
  * [15] E. Hinch, A. Acrivos, Long slender drops in a simple shear flow, Journal of fluid mechanics 98 (2) (1980) 305–328.
  * [16] H. Karam, J. Bellinger, Deformation and breakup of liquid droplets in a simple shear field, Industrial & Engineering Chemistry Fundamentals 7 (4) (1968) 576–581.
  * [17] R. W. Flumerfelt, Drop breakup in simple shear fields of viscoleastic fluids, Industrial & Engineering Chemistry Fundamentals 11 (3) (1972) 312–318.
  * [18] H. A. Stone, B. Bentley, L. Leal, An experimental study of transient effects in the breakup of viscous drops, Journal of Fluid Mechanics 173 (1986) 131–158.
  * [19] S. Guido, M. Villone, Three-dimensional shape of a drop under simple shear flow, Journal of rheology 42 (2) (1998) 395–415.
  * [20] H. P. Grace, Dispersion phenomena in high viscosity immiscible fluid systems and application of static mixers as dispersion devices in such systems, Chemical Engineering Communications 14 (3-6) (1982) 225–277.
  * [21] H. A. Stone, L. G. Leal, The influence of initial deformation on drop breakup in subcritical time-dependent flows at low reynolds numbers, Journal of fluid mechanics 206 (1989) 223–263.
  * [22] A. Vananroye, P. Van Puyvelde, P. Moldenaers, Effect of confinement on droplet breakup in sheared emulsions, Langmuir 22 (9) (2006) 3972–3974.
  * [23] A. Vananroye, P. Van Puyvelde, P. Moldenaers, Effect of confinement on the steady-state behavior of single droplets during shear flow, Journal of rheology 51 (1) (2007) 139–153.
  * [24] M. R. Kennedy, C. Pozrikidis, R. Skalak, Motion and deformation of liquid drops, and the rheology of dilute emulsions in simple shear flow, Computers & fluids 23 (2) (1994) 251–278.
  * [25] E. Toose, B. Geurts, J. Kuerten, A boundary integral method for two-dimensional (non)-newtonian drops in slow viscous flow, Journal of non-newtonian fluid mechanics 60 (2-3) (1995) 129–154.
  * [26] W. Uijttewaal, E. Nijhof, The motion of a droplet subjected to linear shear flow including the presence of a plane wall, Journal of Fluid Mechanics 302 (1995) 45–63.
  * [27] J. Li, Y. Y. Renardy, M. Renardy, Numerical simulation of breakup of a viscous drop in simple shear flow through a volume-of-fluid method, Physics of fluids 12 (2) (2000) 269–282.
  * [28] A. Amani, N. Balcázar, J. Castro, A. Oliva, Numerical study of droplet deformation in shear flow using a conservative level-set method, Chemical Engineering Science 207 (2019) 153–171.
  * [29] H. Xi, C. Duncan, Lattice boltzmann simulations of three-dimensional single droplet deformation and breakup under simple shear flow, Physical Review E 59 (3) (1999) 3022.
  * [30] R. G. van der Sman, S. Van der Graaf, Emulsion droplet deformation and breakup with lattice boltzmann model, Computer Physics Communications 178 (7) (2008) 492–504.
  * [31] S. Farokhirad, T. Lee, J. F. Morris, Effects of inertia and viscosity on single droplet deformation in confined shear flow, Communications in Computational Physics 13 (3) (2013) 706–724.
  * [32] A. Komrakova, O. Shardt, D. Eskin, J. Derksen, Lattice boltzmann simulations of drop deformation and breakup in shear flow, International Journal of Multiphase Flow 59 (2014) 24–43.
  * [33] B. Huang, H. Liang, J. Xu, Lattice boltzmann simulation of binary three-dimensional droplet coalescence in a confined shear flow, Physics of Fluids 34 (3) (2022) 032101.
  * [34] D. Wang, D. S. Tan, B. C. Khoo, Z. Ouyang, N. Phan-Thien, A lattice boltzmann modeling of viscoelastic drops’ deformation and breakup in simple shear flows, Physics of Fluids 32 (12) (2020) 123101.
  * [35] Y. Zong, C. Zhang, H. Liang, L. Wang, J. Xu, Modeling surfactant-laden droplet dynamics by lattice boltzmann method, Physics of Fluids 32 (12) (2020) 122105\.
  * [36] J. Monaghan, Smoothed Particle Hydrodynamics and Its Diverse Applications, Annual Review of Fluid Mechanics 44 (1) (2012) 323–346.  arXiv:1012.1885.
  * [37] T. Ye, D. Pan, C. Huang, M. Liu, Smoothed particle hydrodynamics (sph) for complex fluid flows: Recent developments in methodology and applications, Physics of Fluids 31 (1) (2019) 011301.
  * [38] J. P. Morris, Simulating surface tension with smoothed particle hydrodynamics, International journal for numerical methods in fluids 33 (3) (2000) 333–353.
  * [39] X. Y. Hu, N. A. Adams, A multi-phase sph method for macroscopic and mesoscopic flows, Journal of Computational Physics 213 (2) (2006) 844–861.
  * [40] Z.-B. Wang, R. Chen, H. Wang, Q. Liao, X. Zhu, S.-Z. Li, An overview of smoothed particle hydrodynamics for simulating multiphase flow, Applied Mathematical Modelling 40 (23-24) (2016) 9625–9655.
  * [41] M. Zhang, Simulation of surface tension in 2d and 3d with smoothed particle hydrodynamics method, Journal of Computational Physics 229 (19) (2010) 7238–7259.
  * [42] A. M. Tartakovsky, A. Panchenko, Pairwise force smoothed particle hydrodynamics model for multiphase flow: surface tension and contact line dynamics, Journal of Computational Physics 305 (2016) 1119–1146.
  * [43] Q. Yang, J. Yao, Z. Huang, M. Asif, A comprehensive sph model for three-dimensional multiphase interface simulation, Computers & Fluids 187 (2019) 98–106.
  * [44] Z. Moinfar, S. Vahabi, M. Vahabi, Numerical simulation of drop deformation under simple shear flow of giesekus fluids by sph, International Journal of Numerical Methods for Heat & Fluid Flow (ahead-of-print) (2022).
  * [45] M. Vahabi, The effect of thixotropy on deformation of a single droplet under simple shear flow, Computers & Mathematics with Applications 117 (2022) 206–215.
  * [46] R. Saghatchi, M. Ozbulut, M. Yildiz, Dynamics of double emulsion interfaces under the combined effects of electric field and shear flow, Computational Mechanics 68 (4) (2021) 775–793.
  * [47] M. Hirschler, G. Oger, U. Nieken, D. Le Touzé, Modeling of droplet collisions using smoothed particle hydrodynamics, International Journal of Multiphase Flow 95 (2017) 175–187.
  * [48] X. Xu, T. Tang, P. Yu, A modified sph method to model the coalescence of colliding non-newtonian liquid droplets, International Journal for Numerical Methods in Fluids 92 (5) (2020) 372–390.
  * [49] A. Zhang, P. Sun, F. Ming, An sph modeling of bubble rising and coalescing in three dimensions, Computer Methods in Applied Mechanics and Engineering 294 (2015) 189–209.
  * [50] J. U. Brackbill, D. B. Kothe, C. Zemach, A continuum method for modeling surface tension, Journal of computational physics 100 (2) (1992) 335–354.
  * [51] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, G. Zanetti, Modelling merging and fragmentation in multiphase flows with surfer, Journal of computational physics 113 (1) (1994) 134–147.
  * [52] D. J. Price, Smoothed particle hydrodynamics and magnetohydrodynamics, Journal of Computational Physics 231 (3) (2012) 759–794.
  * [53] J. P. Morris, P. J. Fox, Y. Zhu, Modeling low reynolds number incompressible flows using sph, J. Comput. Phys. 136 (1) (1997) 214–226.
  * [54] K. S. Sheth, C. Pozrikidis, Effects of inertia on the deformation of liquid drops in simple shear flow, Comput. Fluids 24 (2) (1995) 101–119.
  * [55] H. Zhou, C. Pozrikidis, The flow of suspensions in channels: single files of drops, Phys. Fluids A 5 (2) (1993) 311–324.
  * [56] M. Shapira, S. Haber, Low reynolds number motion of a droplet in shear flow including wall effects, International journal of multiphase flow 16 (2) (1990) 305–321.
  * [57] V. Sibillo, G. Pasquariello, M. Simeone, V. Cristini, S. Guido, Drop deformation in microconfined shear flow, Physical review letters 97 (5) (2006) 054502.


# Production of $d_{N\Omega}$ dibaryon in kaon induced reactions

Jing Liu1<EMAIL_ADDRESS>Qi-Fang Lü2,3,4<EMAIL_ADDRESS>Chun-Hua
Liu1<EMAIL_ADDRESS>Dian-Yong Chen1,5111Corresponding author
<EMAIL_ADDRESS>Yu-Bing Dong6,7<EMAIL_ADDRESS>1School of Physics,
Southeast University, Nanjing 210094, China
2Department of Physics, Hunan Normal University, and Key Laboratory of Low-
Dimensional Quantum Structures and Quantum Control of Ministry of Education,
Changsha 410081, China
3Synergetic Innovation Center for Quantum Effects and Applications (SICQEA),
Hunan Normal University, Changsha 410081, China
4Research Center for Nuclear Physics (RCNP), Ibaraki, Osaka 567-0047, Japan
5Lanzhou Center for Theoretical Physics, Lanzhou University, Lanzhou 730000,
China
6Institute of High Energy Physics, Chinese Academy of Sciences, Beijing
100049, China
7School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 101408, China

###### Abstract

In this work, we propose to investigate the $d_{N\Omega}$ dibaryon production
in the process $K^{-}p\rightarrow d_{N\Omega}\bar{\Xi}^{0}$ by utilizing the
kaon beam with the typical momentum to be around 10 GeV, which may be
available at COMPASS, OKA@U-70 and SPS@CERN. The cross sections for
$K^{-}p\rightarrow d_{N\Omega}\bar{\Xi}^{0}$ are estimated and in particular,
the magnitude of the cross sections are evaluated to be several hundreds
nanobarn at $P_{K}=20$ GeV. Considering that $d_{N\Omega}$ dibaryon dominantly
decay into $\Xi\Lambda$ and $\Xi\Sigma$, we also estimate the cross sections
for $K^{-}p\to\Xi^{0}\Lambda\bar{\Xi}^{0}$ and
$K^{-}p\to\Xi^{-}\Sigma^{+}\bar{\Xi}^{0}$, where the $d_{N\Omega}$ dibaryon
can be observed in the invariant mass distributions of $\Xi^{0}\Lambda$ and
$\Xi^{-}\Sigma^{+}$, respectively.

###### pacs:

13.25.GV, 13.75.Lb, 14.40.Pq

##  I Introduction

In the strict sense, dibaryon is a sort of particle with baryon number $B=2$,
which is composed of six quarks , typically two baryons. As a typical
dibaryon, a deuteron consists of one proton and one neutron with a binding
energy to be $2.22$ MeV. Besides the deuteron, searching for the dibaryon
states composed of other baryons becomes one of intriguing topics of hadron
physics in recent decades with the developments of experimental techniques and
the accumulations of experimental data.

The theoretical investigations of the dibaryon could date back to the year of
1964, when the dibaryon states were studied by Dyson and Xuong Dyson:1964xwa
with SU(6) theory, and then in 1977 Jaffe predicted the existence of the $H$
and $H^{*}$ particles with strangeness $S=-2$, which could be the bound state
of $\Lambda\Lambda$ and $\Sigma\Sigma$, respectively. The investigations in
two different quark models indicated that there should exist dibaryon states
with strangeness to be $-3$ Goldman:1987ma , which were stable for strong
decay Goldman:1987ma . And later, the mass spectrum of the low-lying dibaryons
with strangeness to be -1 were evaluated in the quark model Konno:1988nb .
Furthermore, some other models had been extended to study the dibaryons, for
examples, the quark-cluster model Oka:1988yq , Skyrme Model Kopeliovich:1992sa
, quark potential model Wang:1995bg , the chiral SU(3) quark model
Zhang:2000sv , quark delocalization and color screening model Wu:1996fm ;
Wang:1992wi ; Chen:2007qn ; Ping:2000cb ; Pang:2003ty  and realistic
phenomenological nucleon-nucleon interaction models Julia-Diaz:2004ict ;
Froemel:2004ea . In these model, a series of dibaryon states have been
investigated, such as the nonstrange dibaryon $d^{\ast}$ with
$I(J^{P})=0(3^{+})$, dibaryons composed of $N\Xi^{\prime}$, $N\Xi_{c}$,
$N\Xi_{cc}$, $\Xi_{cc}\Xi_{cc}$ Julia-Diaz:2004ict ; Froemel:2004ea ,
$N\Omega$ and $\Delta\Omega$ Dai:2006dgq .

On the experimental side, the first breakthrough was the observations of
$d^{\ast}(2380)$, which was firstly observed in the cross sections for $np\to
d\pi\pi$ by CELSIUS/WASA Collaboration in 2009 Bashkanov:2008ih . The
discovery of $d^{\ast}(2380)$ by CELSIUS/WASA Collaboration has encouraged the
experimentists to search $d^{\ast}(2380)$ in more processes. With much better
statistics and precision data sample, the WASA-at-COSY Collaboration found
that the observed angular distributions for the deuterons and pions of the
process $np\to d\pi^{0}\pi^{0}$ in the center-of-mass system clearly preferred
$J=3$ WASA-at-COSY:2011bjg ; WASA-at-COSY:2012seb  for $d^{\ast}(2380)$, and
then the $I(J^{P})$ quantum numbers were determined to be $0(3^{+})$. Later
on, the properties of $d^{\ast}(2380)$ had been investigated in various
unpolarized $np$ collision processes with more precise data samples, for
example, $np\rightarrow np\pi^{0}\pi^{0}$, $np\rightarrow pp\pi^{0}\pi^{-}$,
$np\rightarrow d\pi^{0}\pi^{0}$ WASA-at-COSY:2014qkg ; WASA-at-COSY:2014squ ;
WASA-at-COSY:2013fzt , and polarized $np$ scattering process WASA-at-
COSY:2014dmv .

The experimental observations of $d^{\ast}(2380)$ has attracted a lot of
theorists’ attentions to explore its properties Shi:2019dpo ; Huang:2014kja ;
Dong:2015cxa ; Dong:2016rva ; Dong:2017geu ; Huang:2015nja ; Gongyo:2020pyy ;
Huang:2022qgx . There are also dissenting voices, for example, in Refs.
Ikeno:2021frl ; Molina:2021bwp , the authors claimed that the peak structure
should be tied to a triangle singularity in the last step of the reaction.
Anyhow, the experimental breakthrough has stimulate theorists’ interests in
dibaryons states. Some dibaryon systems had been investigated in various
models. For examples, in a three-body hadronic model, the authors in Ref.
Gal:2014zia  calculated the $N\Delta$ and $\Delta\Delta$ dibaryon states. In
the one pion exchange potential model, the $H$-like $\Lambda_{c}\Lambda_{c}$
Meguro:2011nr , $\Lambda_{c}N$ Liu:2011xc  dibaryons were researched. The
possible $\Lambda_{c}\Lambda_{c}/\Lambda_{b}\Lambda_{b}$ Huang:2013rla , and
$N\Sigma_{c,b}$ Huang:2013zva  dibaryon states were studied in the quark
delocalization color screening model, and the possible $\Delta^{0}\Delta^{0}$,
$\Omega\Omega$, $\Xi\Xi$ dibaryons were investigated in Refs. Lu:2020qme ;
Liu:2021gva ; Mutuk:2022zgn ; Chen:2019vdh ; Brambilla:2022ura .

|  |
---|---|---
(a) |  | (b)
Figure 1: Diagrams contributing to the process of $K^{-}p\to
d_{N\Omega}\bar{\Xi}^{0}$, where the $d_{N\Omega}$ is considered as a S-wave
$N\Omega$ dibaryon with $J^{P}=2^{+}$.

As for the dibaryon states with strangeness $S=-3$, the possibilities of the
existence of such kind of dibaryons have been evaluated by various model. For
example, the estimations in Ref. Li:1999bc  indicated that there might exist
two bound states $N\Omega$ and $\Delta\Omega$ dibaryons, respectively, and in
Ref. Dai:2006dgq , the authors found that the $N\Omega$ and $\Delta\Omega$
were weakly bound systems in the chiral quark model. The QCD sum rule
estimations in Ref. Chen:2021hxs  indicate that the $N\Omega$ dibaryon with
$J^{P}=2^{+}$ was stable for strong decay. In Ref. HALQCD:2014okw , the HAL
QCD Collaboration calculated the $N\Omega$ potential in $2+1$ flavor Lattice
QCD and one bound state with the binding energy to be about 20 MeV was found.
In the year of 2019, the HAL QCD Collaboration updated their estimations of
Ref. HALQCD:2014okw  near the physical point and they found the binding energy
of the $p\Omega(^{5}S_{2})$ became $2.46(0.34)(^{+0.04}_{-0.11})$ MeV.
Simulated by the most recent Lattice QCD estimations, the authors in Ref.
Zhang:2020dma  estimated the productions of the $\Omega$-dibaryons by
utilizing a dynamical coalescence mechanism for the relativistic heavy-ion
collisions, and the strong decays of $d_{N\Omega}$ into conventional hadrons
were evaluted Xiao:2020alj .

Actually, experimentally producing the $d_{N\Omega}$ dibaryon is the first
step of investigating its properties. Thus, searching $d_{N\Omega}$
experimentally becomes a pressing task. The key point of producing
$d_{N\Omega}$ dibaryon is the production of $\Omega$ baryon. In the high
energy heavy-ion collision or $pp$ collision processes, a large quantity of
quarks with different flavors can be produced, where three strange quarks have
chance to form a $\Omega$ baryon and produce a $d_{N\Omega}$ dibaryon by
interacting with a nucleon. In the year of 2018, the STAR Collaboration at
RIHC investigated the proton-$\Omega$ interaction by measuring the
corresponding correlation function in heavy-ion collisions at
$\sqrt{s_{NN}}=200$ GeV STAR:2018uho . By comparing the measured correlation
ratio with the theoretical calculations, they concluded that the measurements
slightly favored a proton-$\Omega$ bound system with a binding energy of
$~{}27$ MeV STAR:2018uho . Similarly, the ALICE Collaboration also proposed to
investigate the strong interaction among hadrons, including proton-$\Omega$,
by using the ultrarelativistic proton-proton collisions at LHC ALICE:2020mfd .

Besides the high energy heavy-ion collision and $pp$ collision processes, the
$\Omega$ baryon can also be produced by kaon induced reactions since there is
already a strange quark in the kaon. In Ref. Aoki:2021cqa , a project of the
extension of the Hadron Experimental Facility at J-PARC was proposed. By
utilizing the secondary beam of kaon with the typical momentum to be around 3
GeV, the $\Omega$ baryon can be produced via the
$K^{-}p\to\Omega^{-}K^{+}\bar{K}^{(\ast)0}$ reaction and then the
$d_{N\Omega}$ dibaryon can be investigated in the invariant mass spectrum of
$\Xi^{-}\Lambda$ of the process $\Omega^{-}d\to\Xi^{-}\Lambda p$.

In addition to the secondary beam of kaon at J-PARC Nagae:2008zz , the high
energy kaon beam with high quality are also available at COMPASS
Nerling:2012er , OKA@U-70 Obraztsov:2016lhp  and SPS@CERN Velghe:2016jjw .
These high energy beam, especially ones with momentum to be around 10 GeV, may
provide us another approach to directly produce $d_{N\Omega}$ via the reaction
$K^{-}p\to d_{N\Omega}\bar{\Xi}^{0}$ and then the $d_{N\Omega}$ dibaryon can
decay into $\Xi\Lambda$ and $\Xi\Sigma$. Thus, one can detect $d_{N\Omega}$ in
the $\Xi\Lambda$ and $\Xi\Sigma$ invariant mass spectra of the processes
$K^{-}p\to\bar{\Xi}^{0}\Xi\Lambda$ and $K^{-}p\to\bar{\Xi}^{0}\Xi\Sigma$,
respectively. In the present work, we evaluate the possibility of observing
$d_{N\Omega}$ dibaryon in these processes by estimating their cross sections.

This work is organized as follows. After the introduction, the mechanism of
$d_{N\Omega}$ production in kaon induced reactions is presented. In Section
III, the cross sections for the processes $K^{-}p\rightarrow
d_{N\Omega}\bar{\Xi}^{0}$, $K^{-}p\rightarrow\Xi^{0}\Lambda\bar{\Xi}^{0}$ and
$K^{-}p\rightarrow\Xi^{-}\Sigma^{+}\bar{\Xi}^{0}$ are presented, and the last
section is devoted to a short summary.

Figure 2: Diagrams contributing to $K^{-}p\to\Xi^{0}\Lambda\bar{\Xi}^{0}$
(diagram (a)) and $K^{-}p\to\Xi^{-}\Sigma^{+}\bar{\Xi}^{0}$ (diagram (b)).

.

##  II Production Processes

###  II.1 The $K^{-}p\to d_{N\Omega}\bar{\Xi}^{0}$ process

In the present work, we consider the $d_{N\Omega}$ is a dibaryon composed of
$N$ and $\Omega$ with the $J(J^{P})=\frac{1}{2}(2^{+})$, i.e.,
$|d_{N\Omega}^{0}\rangle=|p\Omega\rangle$ and
$|d_{N\Omega}^{-}\langle=|n\Omega\rangle$. The dibayon $d_{N\Omega}$ can be
produced in the high energy $K^{-}p$ interaction process. As indicated in Ref.
Xiao:2020alj , the dibaryon $d_{N\Omega}$ could decay into $\Lambda\Xi^{0}$,
which indicate the strong coupling between $d_{N\Omega}$ and $\Lambda\Xi^{0}$,
thus, there exists the $s$ channel contributions to
$K^{-}p\to\bar{\Xi}^{0}d_{N\Omega}$ as show in Fig. 1-(a), where the initial
$K^{-}p$ and final $d_{N\Omega}\bar{\Xi}^{0}$ are connected by $\Lambda$.
Moreover, there is the contributions from the $\Omega$ baryon exchange
process, which should be the $u$ channel contribution in the strict sense. In
the process $K^{-}p\to\bar{\Xi}^{0}d_{N\Omega}$, there is no $t$ channel
digram in the tree level. Empirically, in high energy $K^{-}p/\pi p$
scattering process, the contributions from the $s$ channels are strongly
suppressed Liu:2021ojf ; Liu:2020ruo ; Huang:2015xud ; Liu:2008qx , thus the
$d_{N\Omega}$ production in high energy $K^{-}p$ scattering should occur
dominantly by exchanging a $\Omega$ baryon as presented in Fig. 1-(b), while
the $s$ channel contribution is ignored.

In the present work, we estimated the cross sections for $K^{-}p\to
d_{N\Omega}\bar{\Xi}$ in an effective Lagrangian approach. The interaction of
the $d_{N\Omega}$ dibayon and its components can be described as Xiao:2020alj
,

$\displaystyle\mathcal{L}_{d_{N\Omega}N\Omega}=g_{d_{N\Omega}N\Omega}d^{\mu\nu^{\dagger}}_{N\Omega}\bar{\Omega}_{\mu}\gamma_{\nu}N^{c}+H.c.,$
(1)

where $N^{c}=C\bar{N}^{T}$, $\bar{N}^{c}=N^{T}C$, and
$C=i\gamma^{2}\gamma^{0}$ is the charge-conjugation matrix, $T$ is the
transpose transformation operator. From the above effective Lagrangian, the
tensor field of $d_{N\Omega}$ can be constructed by the appropriate
combination of a Dirac field for the nucleon and the a Rarita Schwinger field
for $\Omega$. The polarization tensor $\epsilon^{\mu\nu}(\vec{p},\lambda)$
could be constructed by the combination of the Dirac field for spin-$1/2$ and
a Rarita Schwinger field for spin-$3/2$ i.e,

$\displaystyle\epsilon^{\mu\nu}(\vec{p},\lambda)=\sum_{\alpha,\beta}\langle\frac{3}{2}\alpha\frac{1}{2}\beta|2\lambda\rangle\psi^{\mu}_{\alpha}(\vec{p})\gamma^{\nu}\psi_{\beta}^{c}(\vec{p})$
(2)

with $\lambda=(\pm 2,\pm 1,0)$, $\alpha=(\pm 3/2,\pm 1/2)$ and $\beta=\pm
1/2$, respectively. The polarization tensor satisfy,

$\displaystyle
p_{\mu}\epsilon^{\mu\nu}(\vec{p},\lambda)=p_{\nu}\epsilon^{\mu\nu}(\vec{p},\lambda)=0,$
$\displaystyle\epsilon^{\mu\nu}(\vec{p},\lambda)=\epsilon^{\nu\mu}(\vec{p},\lambda),\
\ \ \epsilon^{\mu}_{\mu}(\vec{p},\lambda)=0$
$\displaystyle\epsilon^{\mu\nu\ast}(\vec{p},\lambda)\epsilon_{\mu\nu}(\vec{p},\lambda^{\prime})=\delta_{\lambda\lambda^{\prime}}$
(3)

The effective Lagrangian for $\Omega\Xi K$ interaction reads Schutz:1994ue ;
Ronchen:2012eg ; Matsuyama:2006rp ; He:2017aps ; Machleidt:1987hj ; Liu:2001ce
,

$\displaystyle\mathcal{L}_{\Omega\Xi K}=\frac{g_{\Omega\Xi
K}}{m_{\pi}}\partial_{\beta}K\bar{\Omega}^{\beta}\Xi+H.c..$ (4)

With the above effective Lagrangians, we can obtain the amplitude
corresponding to Fig. 1-(b), which is,

$\displaystyle\mathcal{M}$ $\displaystyle=$ $\displaystyle\frac{g_{\Omega\Xi
K}}{m_{\pi}}\left(-ip_{1\beta}\right)\ g_{d_{N\Omega}N\Omega}\
d_{N\Omega}^{\mu\nu}F(k^{2},m^{2}_{\Omega})$ (5) $\displaystyle\times$
$\displaystyle\left[\bar{u}^{c}(p_{2},m_{2})\gamma_{\nu}S(k,m_{\Omega})_{\mu\beta}\nu(p_{3},m_{3})\right],$

with $S(k,m_{\Omega})_{\mu\beta}$ to be the propagator of the $\Omega$ baryon,
which is,

$\displaystyle
S(k,m_{\Omega})_{\mu\beta}=i\frac{k\\!\\!\\!/+m_{\Omega}}{k^{2}-m_{\Omega}^{2}}\left[-g_{\mu\beta}+\frac{1}{3}\gamma_{\mu\beta}+\frac{2k_{\mu}k_{\beta}}{3m_{\Omega}^{2}}\right.$
$\displaystyle\left.+\frac{\gamma_{\mu}k_{\beta}-\gamma_{\beta}k_{\mu}}{3m_{\Omega}}\right].$
(6)

To depict the internal structure and the off shell effect of the exchanged
$\Omega$ baryon, we introduce a form factor $F(k^{2},m^{2}_{\Omega})$ in the
amplitude and its concrete form will be discussed later. With the amplitude in
Eq. (5), one can obtain the cross section for $K^{-}p\to
d_{N\Omega}\bar{\Xi}^{0}$ by,

$\displaystyle\frac{d\sigma}{d\cos\theta}=\frac{1}{32\pi
s}\frac{|\vec{p}_{f}|}{|\vec{p}_{i}|}\left(\frac{1}{2}|\overline{\mathcal{M}}|^{2}\right),$
(7)

where $s=(p_{1}+p_{2})^{2}$ is the center of mass energy and $\theta$ is the
scattering angle, which refers to the angle of outgoing $d_{N\Omega}$ and the
kaon beam direction in the center of mass frame. The $\vec{p}_{i}$ and
$\vec{p}_{f}$ are three momentum of the initial kaon beam and the final
$d_{N\Omega}$ dibaryon in the center of mass frame, respectively.

###  II.2 The $K^{-}p\to\Lambda\Xi^{0}\bar{\Xi}^{0}$ and
$K^{-}p\to\Sigma^{+}\Xi^{-}\bar{\Xi}^{0}$ processes

As indicated in Ref. Xiao:2020alj , the dibayon $d_{N\Omega}$ dominantly decay
into $\Lambda\Xi$ and $\Xi\Sigma$. Thus one can detect $d_{N\Omega}$ in the
invariant mass spectrum of $\Lambda\Xi^{0}$ and $\Xi^{-}\Sigma^{+}$ of the
processes $K^{-}p\rightarrow\Xi^{0}\Lambda\bar{\Xi}^{0}$ and
$K^{-}p\rightarrow\Xi^{-}\Sigma^{+}\bar{\Xi}^{0}$, respectively. To depict
these processes, additional effective Lagrangians related to
$d_{N\Omega}\Xi\Lambda$ and $d_{N\Omega}\Xi\Sigma$ are involved. Since $\Xi$,
$\Sigma$ and $\Lambda$ have the same $J^{P}$ quantum numbers, thus these two
effective Lagrangians have the same form, which is,

Figure 3: The cross sections for the process $K^{-}p\to
d_{N\Omega}\bar{\Xi}^{0}$ depending on the beam energy (diagram (a)), and the
differential cross sections depending on $\cos(\theta)$ (diagram (b)).

$\displaystyle\mathcal{L}_{d_{N\Omega}Y_{1}Y_{2}}$ $\displaystyle=$
$\displaystyle
i\frac{G_{d_{N\Omega}Y_{1}Y_{2}}}{2M_{Y_{1}}}\bar{Y_{1}}^{c}\left(\gamma_{\mu}\partial_{\nu}+\gamma_{\nu}\partial_{\mu}\right)Y_{2}d_{N\Omega}^{\mu\nu}$
(8)
$\displaystyle+\frac{F_{d_{N\Omega}Y_{1}Y_{2}}}{(2M_{Y_{1}})^{2}}\partial_{\mu}\bar{Y}_{1}^{c}\partial_{\nu}Y_{2}d_{N\Omega}^{\mu\nu}+H.c,$

where $Y_{1}$, $Y_{2}$ could be $\Xi$, $\Sigma$ and $\Lambda$. Similar to the
case of tensor meson, we can choose $F_{d_{N\Omega}Y_{1}Y_{2}}=0$ with tensor
dominance hypothesis Renner:1970sbf , and the values of the couplings
$G_{d_{N\Omega}Y_{1}Y_{2}}$ will be discussed in next section. With this
additional effective Lagrangian, we can obtain the amplitudes corresponding to
Fig. 2-(a), which are,

$\displaystyle\mathcal{M}_{a}$ $\displaystyle=$
$\displaystyle\Bigg{[}g_{d_{N\Omega}N\Omega}\bar{u}^{c}\left(p_{2},m_{2}\right)\gamma^{\nu}S(k,m_{\Omega})_{\mu\beta}v(p_{3},m_{3})\Bigg{]}$
(9) $\displaystyle\times$ $\displaystyle\left[\frac{g_{\Omega\Xi
K}}{m_{\pi}}(-ip_{1\beta})\right]\mathcal{P}_{d_{N\Omega}}^{\mu\nu\lambda\omega}\Big{(}q,m_{d_{N\Omega}},\Gamma_{d_{N\Omega}}\Big{)}\Bigg{[}i\frac{G_{d_{N\Omega}\Xi\Lambda}}{2m_{\Xi}}$
$\displaystyle\times$
$\displaystyle\bar{u}^{c}(p_{5},m_{5})\left(\gamma_{\lambda}(-ip_{4\omega})+\gamma_{\omega}(-ip_{4\lambda})\right)u(p_{4},m_{4})\Bigg{]}$
$\displaystyle\times$ $\displaystyle
F\left(k^{2},m^{2}_{\Omega}\right)F\left(q^{2},m^{2}_{d_{N\Omega}}\right),$

where
$\mathcal{P}_{d_{N\Omega}}^{\mu\nu\lambda\omega}(q,m_{d_{N\Omega}},\Gamma_{d_{N\Omega}})$
is the propagator of the dibaryon $d_{N\Omega}$, and its concrete form is,

$\displaystyle\mathcal{P}_{d_{N\Omega}}^{\mu\nu\lambda\omega}(q,m_{d_{N\Omega}},\Gamma_{d_{N\Omega}})=\frac{i}{q^{2}-m_{d_{N\Omega}}^{2}+im_{d_{N\Omega}}\Gamma_{d_{N\Omega}}}$
$\displaystyle\quad\quad\quad\quad\quad\times\left[\frac{1}{2}\left(\tilde{g}_{\mu\lambda}\tilde{g}_{\nu\omega}+\tilde{g}_{\mu\omega}\tilde{g}_{\nu\lambda}\right)-\frac{1}{3}\tilde{g}_{\mu\nu}\tilde{g}_{\lambda\omega}\right],\qquad$
(10)

with $\tilde{g}^{\mu\nu}=-g^{\mu\nu}+q^{\mu}q^{\nu}/m^{2}$. In the above
amplitudes, an addition form factor $F(q^{2},m^{2}_{d_{N\Omega}})$ is
introduced to depict the internal structure and off shell effects of the
$d_{N\Omega}$ and its concrete form will be discussed in the next section. In
the same way, one can obtain the amplitude of
$K^{-}p\to\bar{\Xi}^{0}\Xi^{-}\Sigma^{+}$ corresponding to Fig. 2-(b). With
the amplitudes in Eq. (9), we can obtain the cross sections for the $2\to 3$
processes by,

$\displaystyle
d\sigma=\frac{1}{8(2\pi)^{4}}\frac{1}{\Phi}|\mathcal{M}|^{2}dp_{5}^{0}dp_{3}^{0}d\cos\theta
d\eta,$ (11)

where $\Phi=2\sqrt{\lambda(s,m_{1}^{2},m_{2}^{2})}=4|\vec{p}_{1}|\sqrt{s}$
with $\vec{p}_{1}$ to be the three momentums of the incident particle $K^{-}$.
$p_{3}^{0}$ is the energy of the outgoing $\bar{\Xi}^{0}$, while $p_{5}^{0}$
is the energy of $\Lambda$ in the $K^{-}p\to\Lambda\Xi^{0}\bar{\Xi}^{0}$
process, and the energy of $\Sigma^{+}$ in the
$K^{-}p\to\bar{\Xi}^{0}\Xi^{-}\Sigma^{+}$ process.

##  III NUMERICAL RESULTS and discussions

Figure 4: Cross sections for $K^{-}p\to\Xi^{0}\Lambda\bar{\Xi}^{0}$ (diagram
(a)) and $K^{-}p\to\Xi^{-}\Sigma^{+}\bar{\Xi}^{0}$ (diagram (b)) depending on
the beam energy.

###  III.1 Form factor and Coupling Constants

In the present work, we introduce two form factors to depict the internal
structures and off-shell effects of the exchanging $\Omega$ baryon and
intermediate $d_{N\Omega}$ dibaryon. The specific expression of the factor is
Chen:2013cpa ; Chen:2014ccr ; :2019rva; Kim:2015ita ,

$\displaystyle
F(q^{2},m^{2})=\frac{\Lambda^{4}}{(m^{2}-q^{2})^{2}+\Lambda^{4}},$ (12)

where $\Lambda$ is model parameter. In principle, the value of the model
parameter $\Lambda$ should be determined by comparing the theoretical
estimations with the corresponding experimental measurements. However, the
experimental data for the cross sections for the discussed processes are not
available at present. In Ref. Kim:2015ita , the authors investigated the cross
sections for $\pi^{-}p\to K^{\ast 0}\Lambda$ and the parameter $\Lambda$ is
determined to be 0.55 GeV for the $t$ channel and 0.60 GeV for the $u/s$
channels by comparing the estimated cross sections with the experimental data.
With this parameter, they also extended to estimate the cross sections for
$\pi^{-}p\to D^{\ast-}\Lambda_{c}^{+}$. In the present work, we take a very
similar parameter range, which is $0.55\ \mathrm{GeV}<\Lambda<0.65\
\mathrm{GeV}$, to calculate the cross sections for $K^{-}p\to
d_{N\Omega}\bar{\Xi}^{0}$ process.

Before the estimations of the cross sections for the discussed processes, the
relevant coupling constants should be clarified. As for the coupling constant
$g_{d_{N\Omega}N\Omega}$, it can be estimated by the compositeness condition
of the composite state. In Ref. Xiao:2020alj , the coupling constant
$g_{d_{N\Omega}N\Omega}$ is estimated to be about $1.88\sim 2.38$ with the
variation of the model parameter, where the binding energy is set to be $2.46$
MeV. In the present estimation, we take $g_{d_{N\Omega}N\Omega}=1.97$. With
this coupling constant, the partial decay widths of
$d_{N\Omega}\to\Lambda\Xi^{0}$ and $d_{N\Omega}\to\Sigma^{+}\Xi^{-}$ are
estimated to be $582$ and $22.8$ keV Xiao:2020alj , respectively. Together
with the effective Lagrangian in Eq. (8), one can obtain amplitudes of
$d_{N\Omega}\to\Lambda\Xi^{0}$ and $d_{N\Omega}\to\Sigma^{+}\Xi^{-}$, and then
the corresponding partial width can be estimated by,

$\displaystyle\Gamma_{d_{N\Omega}\to...}=\frac{1}{(2J+1)8\pi}\frac{|\vec{k_{f}}|}{M_{d_{N\Omega}}^{2}}|\overline{\mathcal{M}_{d_{N\Omega}\to...}}|^{2},$
(13)

where $J=2$ is the angular momentum of $d_{N\Omega}$, $|\vec{k}_{f}|$ is the
three momentum of the daughter particles in $d_{N\Omega}$ rest frame. From Eq.
(13) and the partial widths obtained in Ref. Xiao:2020alj , one can estimate
the corresponding effective coupling constants, which are
$g_{d_{N\Omega}\Lambda\Xi}=7.9\times 10^{-2}$ and
$g_{d_{N\Omega}\Sigma\Xi}=2.5\times 10^{-2}$, respectively. As for the
coupling constant $g_{K\Xi\Omega}$, we take $g_{K\Xi\Omega}=2.12$ as in Ref.
Schutz:1994ue .

###  III.2 Cross sections for $K^{-}p\to d_{N\Omega}\bar{\Xi}^{0}$

With the above preparation, we can evaluate the cross sections for $K^{-}p\to
d_{N\Omega}\bar{\Xi}^{0}$ depending on the beam energy and the model parameter
$\Lambda$. In Fig 3-(a), the solid curve refers to the cross sections for
$K^{-}p\to d_{N\Omega}\bar{\Xi}^{0}$ with $\Lambda=0.60$ GeV, while the
uncertainties of the cross sections are obtained by the variation of the
parameter $\Lambda$ from 0.55 GeV to 0.65 GeV. From the figure one can find
the cross sections for $K^{-}p\to d_{N\Omega}\bar{\Xi}^{0}$ increase sharply
near the threshold of $d_{N\Omega}\bar{\Xi}^{0}$, however, when the beam
energy is greater than 9 GeV, the cross sections increase very slowly with the
increase of the beam energy. In the considered parameter range, the cross
sections are estimated to be $404^{+358}_{-202}\ \mathrm{nb}$ at $P_{K}=20$
GeV, where the center value is estimated with $\Lambda=0.60$, while the
uncertainties are resulted from the variation of $\Lambda$ from $0.55$ GeV to
0.65 GeV. In Ref. Kim:2015ita , the cross sections for $\pi p\to
D^{\ast-}\Lambda_{c}^{+}$ is estimated to be about 13 nb, while the present
estimation indicate the cross sections for $K^{-}p\to
d_{N\Omega}\bar{\Xi}^{0}$ can reach up to 400 nb with the same model
parameter, which is about 30 times larger than the one for $\pi p\to
D^{\ast-}\Lambda_{c}^{+}$ Kim:2015ita .

In addition, we also present the differential cross sections depending on
$\cos(\theta)$ in Fig. 3-(b). We select three typical beam energies as
examples, which are $P_{K}=10,\ 12,\ 14$ GeV with $\Lambda=0.6$ GeV. From the
figure one can find that more $d_{N\Omega}$ dibaryon are concentrated in the
forward angle area even in the case of $P_{K}=10$ GeV, which are resulted from
the $\Omega$ exchange.

###  III.3 $K^{-}p\to\Xi^{0}{\Lambda}\bar{\Xi}^{0}$ and
$K^{-}p\to\Xi^{-}\Sigma^{+}\bar{\Xi}^{0}$

Besides the $K^{-}p\to d_{N\Omega}\bar{\Xi}^{0}$ process, we also estimate the
beam energy dependences of the cross sections for
$K^{-}p\to\Xi^{0}{\Lambda}\bar{\Xi}^{0}$ and
$K^{-}p\to\Xi^{-}\Sigma^{+}\bar{\Xi}^{0}$, where $\Xi^{0}\Lambda$ and
$\Xi^{-}\Sigma^{+}$ are the daughter particles of $d_{N\Omega}$. As indicated
in Ref. Xiao:2020alj , the dibaryon $d_{N\Omega}$ dominantly decays into
$\Lambda\Xi^{0}$, and the branching ratio is estimated to be about $95\%$,
thus one can experimentally detect $d_{N\Omega}$ in the $\Lambda\Xi^{0}$
invariant mass distributions of the process
$K^{-}p\to\Xi^{0}{\Lambda}\bar{\Xi}^{0}$ as shown in Fig. 2-(a), where
$\Lambda$ can be reconstructed by $p\pi^{-}$ and $n\pi^{0}$, while $\Xi^{0}$
can be reconstructed by the cascade decay processes
$\Xi^{0}\to\Lambda\pi^{0}\to p\pi^{-}\pi^{0}$ or $\Xi^{0}\to\Lambda\pi^{0}\to
n\pi^{0}\pi^{0}$. Our estimations indicate that the cross sections for
$K^{-}p\to\Xi^{0}\Lambda\bar{\Xi}^{0}$ increase sharply near the threshold and
become very weakly dependent on the beam energy. In particular, the the cross
section is estimated to be $13^{+20}_{-7}\ \mathrm{nb}$ at $P_{K}=20$ GeV,
where the center value is estimated with $\Lambda=0.6$, while the
uncertainties are resulted from the variation of $\Lambda$ from 0.55 GeV to
0.65 GeV.

Furthermore, the branching ratio of $d_{N\Omega}\to\Xi^{-}\Sigma^{+}$ is also
sizable, and the final states are charged, which may be easier to be detected.
Thus, in the present work, we also estimate the cross sections for
$K^{-}p\to\Xi^{-}\Sigma^{+}\bar{\Xi}^{0}$. As shown in Fig. 2-(b), the beam
energy dependences of the cross sections for
$K^{-}p\to\Xi^{-}\Sigma^{+}\bar{\Xi}^{0}$ are very similar to the ones for
$K^{-}p\to\Xi^{0}\Lambda\bar{\Xi}^{0}$ and the magnitude of the cross section
is estimated to be $0.5^{+0.5}_{-0.2}\ \mathrm{nb}$ at $P_{K}=20$ GeV.

##  IV Summary

The production of the dibaryon $d_{N\Omega}$ is the crucial step of
investigating its properties experimentally. The STAR Collaboration at RHIC
have detected the dibaryon $d_{N\Omega}$ by measuring the proton-$\Omega$
correlation function in high energy heavy-ion collision. The ALICE
Collaboration at LHC also proposed to detect $d_{N\Omega}$ in a very similar
processes. By utilizing the secondary kaon beam with the typical momentum to
be around 3 GeV, the Hadron experimental Facility at J-PARC proposed to
produce $d_{N\Omega}$ by two step reactions, i.e.,
$K^{-}p\to\Omega^{-}K^{+}K^{(\ast)0}$, $\Omega d\to\Xi^{-}\Lambda p$, where
the dibaryon $d_{N\Omega}$ is expected to be observed in the $\Xi^{-}\Lambda$
invariant mass spectrum.

Besides the above two kinds of production processes, we propose to directly
produce $d_{N\Omega}$ in the $K^{-}p\to d_{N\Omega}\bar{\Xi}^{0}$ process
using a secondary kaon beam with the typical momentum to be around 10 GeV in
the present work. The cross sections for $K^{-}p\to d_{N\Omega}\bar{\Xi}^{0}$
are estimated and we find the cross section is about several hundreds nanobarn
at $P_{K}=20$ GeV. Moreover, the estimated differential cross sections
indicate that the produced $d_{N\Omega}$ are concentrated in the forward angle
area.

Considering the fact that $d_{N\Omega}$ dominantly decays into $\Xi\Lambda$,
we also estimate the cross sections for $K^{-}p\to\Xi^{0}\Lambda\bar{\Xi}^{0}$
and $K^{-}p\to\Xi^{+}\Sigma^{-}\bar{\Xi}^{0}$, where the $d_{N\Omega}$
dibaryon can be detected in $\Xi^{0}\Lambda$ and $\Xi^{-}\Sigma^{+}$ invariant
mass spectrum, respectively. In particular, the cross sections for
$K^{-}p\to\Xi^{0}\Lambda\bar{\Xi}^{0}$ are estimated to be about ten nanobarn
at $P_{K}=20$ GeV, while the cross section for
$K^{-}p\to\Xi^{-}\Sigma^{+}\bar{\Xi}^{0}$ is about 20 times smaller than the
one for $K^{-}p\to\Xi^{0}\Lambda\bar{\Xi}^{0}$.

ACKNOWLEDGEMENTS

This work is supported by the National Natural Science Foundation of China
under Grants No.11705056, No. 12175037, No.11947224, No.11475192, No.11975245,
and No.U1832173. This work is also supported by the Key Project of Hunan
Provincial Education Department under Grant No. 21A0039, the State Scholarship
Fund of China Scholarship Council under Grant No. 202006725011, the Sino-
German CRC 110 ”Symmetries and the Emergence of Structure in QCD” project by
NSFC under the Grant No. 12070131001, the Key Research Program of Frontier
Sciences, CAS, under the Grant No. Y7292610K1, and the National Key Research
and Development Program of China under Contracts No. 2020YFA0406300.

## References

  * (1) F. Dyson and N. H. Xuong, Phys. Rev. Lett. 13, no.26, 815-817 (1964) doi:10.1103/PhysRevLett.13.815
  * (2) J. T. Goldman, K. Maltman, G. J. Stephenson, Jr., K. E. Schmidt and F. Wang, Phys. Rev. Lett. 59 (1987), 627 doi:10.1103/PhysRevLett.59.627
  * (3) N. Konno, H. Nakamura and H. Noya, Phys. Rev. D 37, 154-158 (1988) doi:10.1103/PhysRevD.37.154
  * (4) M. Oka, Phys. Rev. D 38, 298 (1988) doi:10.1103/PhysRevD.38.298
  * (5) V. B. Kopeliovich, B. Schwesinger and B. E. Stern, Nucl. Phys. A 549, 485-497 (1992) doi:10.1016/0375-9474(92)90661-3
  * (6) F. Wang, J. l. Ping, G. h. Wu, L. j. Teng and J. T. Goldman, Phys. Rev. C 51, 3411 (1995) doi:10.1103/PhysRevC.51.3411 [arXiv:nucl-th/9512014 [nucl-th]].
  * (7) Z. Y. Zhang, Y. W. Yu, C. R. Ching, T. H. Ho and Z. D. Lu, Phys. Rev. C 61, 065204 (2000) doi:10.1103/PhysRevC.61.065204
  * (8) G. H. Wu, L. J. Teng, J. L. Ping, F. Wang and J. T. Goldman, Phys. Rev. C 53, 1161-1166 (1996) doi:10.1103/PhysRevC.53.1161
  * (9) F. Wang, G. h. Wu, L. j. Teng and J. T. Goldman, Phys. Rev. Lett. 69, 2901-2904 (1992) doi:10.1103/PhysRevLett.69.2901 [arXiv:nucl-th/9210002 [nucl-th]].
  * (10) L. Chen, H. Pang, H. Huang, J. Ping and F. Wang, Phys. Rev. C 76, 014001 (2007) doi:10.1103/PhysRevC.76.014001 [arXiv:nucl-th/0703103 [nucl-th]].
  * (11) J. l. Ping, F. Wang and J. T. Goldman, Nucl. Phys. A 688, 871-881 (2001) doi:10.1016/S0375-9474(00)00593-5 [arXiv:nucl-th/0007040 [nucl-th]].
  * (12) H. r. Pang, J. l. Ping, F. Wang, J. T. Goldman and E. g. Zhao, Phys. Rev. C 69, 065207 (2004) doi:10.1103/PhysRevC.69.065207 [arXiv:nucl-th/0306043 [nucl-th]].
  * (13) F. Froemel, B. Julia-Diaz and D. O. Riska, Nucl. Phys. A 750, 337-356 (2005) doi:10.1016/j.nuclphysa.2005.01.022 [arXiv:nucl-th/0410034 [nucl-th]].
  * (14) B. Julia-Diaz and D. O. Riska, Nucl. Phys. A 755, 431-434 (2005) doi:10.1016/j.nuclphysa.2005.03.050 [arXiv:nucl-th/0405061 [nucl-th]].
  * (15) L. R. Dai, D. Zhang, C. R. Li and L. Tong, [arXiv:nucl-th/0606056 [nucl-th]].
  * (16) M. Bashkanov, C. Bargholtz, M. Berlowski, D. Bogoslawsky, H. Calen, H. Clement, L. Demiroers, E. Doroshkevich, D. Duniec and C. Ekstrom, et al. Phys. Rev. Lett. 102, 052301 (2009) doi:10.1103/PhysRevLett.102.052301 [arXiv:0806.4942 [nucl-ex]].
  * (17) P. Adlarson et al. [WASA-at-COSY], Phys. Rev. Lett. 106, 242302 (2011) doi:10.1103/PhysRevLett.106.242302 [arXiv:1104.0123 [nucl-ex]].
  * (18) P. Adlarson et al. [WASA-at-COSY], Phys. Lett. B 721 (2013), 229-236 doi:10.1016/j.physletb.2013.03.019 [arXiv:1212.2881 [nucl-ex]].
  * (19) P. Adlarson et al. [WASA-at-COSY], Phys. Rev. C 88, no.5, 055208 (2013) doi:10.1103/PhysRevC.88.055208 [arXiv:1306.5130 [nucl-ex]].
  * (20) P. Adlarson et al. [WASA-at-COSY], Phys. Lett. B 743, 325-332 (2015) doi:10.1016/j.physletb.2015.02.067 [arXiv:1409.2659 [nucl-ex]].
  * (21) P. Adlarson et al. [WASA-at-COSY], Phys. Rev. C 91, no.1, 015201 (2015) doi:10.1103/PhysRevC.91.015201 [arXiv:1408.5744 [nucl-ex]].
  * (22) P. Adlarson et al. [WASA-at-COSY], Phys. Rev. Lett. 112 (2014) no.20, 202301 doi:10.1103/PhysRevLett.112.202301 [arXiv:1402.6844 [nucl-ex]].
  * (23) P. P. Shi, F. Huang and W. L. Wang, Eur. Phys. J. C 79 (2019) no.4, 314 doi:10.1140/epjc/s10052-019-6809-1 [arXiv:1904.06018 [nucl-th]].
  * (24) F. Huang, Z. Y. Zhang, P. N. Shen and W. L. Wang, Chin. Phys. C 39 (2015) no.7, 071001 doi:10.1088/1674-1137/39/7/071001 [arXiv:1408.0458 [nucl-th]].
  * (25) Y. Dong, P. Shen, F. Huang and Z. Zhang, Phys. Rev. C 91 (2015) no.6, 064002 doi:10.1103/PhysRevC.91.064002 [arXiv:1503.02456 [nucl-th]].
  * (26) Y. Dong, F. Huang, P. Shen and Z. Zhang, Phys. Rev. C 94 (2016) no.1, 014003 doi:10.1103/PhysRevC.94.014003 [arXiv:1603.08748 [hep-ph]].
  * (27) Y. Dong, F. Huang, P. Shen and Z. Zhang, Phys. Lett. B 769 (2017), 223-226 doi:10.1016/j.physletb.2017.03.058 [arXiv:1702.03658 [nucl-th]].
  * (28) F. Huang, P. N. Shen, Y. B. Dong and Z. Y. Zhang, Sci. China Phys. Mech. Astron. 59, no.2, 622002 (2016) doi:10.1007/s11433-015-5767-3 [arXiv:1505.05395 [nucl-th]].
  * (29) S. Gongyo et al. [HAL QCD], Phys. Lett. B 811, 135935 (2020) doi:10.1016/j.physletb.2020.135935 [arXiv:2006.00856 [hep-lat]].
  * (30) F. Huang, Rev. Mex. Fis. Suppl. 3, no.3, 0308031 (2022) doi:10.31349/SuplRevMexFis.3.0308031
  * (31) N. Ikeno, R. Molina and E. Oset, Phys. Rev. C 104, no.1, 014614 (2021) doi:10.1103/PhysRevC.104.014614 [arXiv:2103.01712 [nucl-th]].
  * (32) R. Molina, N. Ikeno and E. Oset, [arXiv:2102.05575 [nucl-th]].
  * (33) A. Gal and H. Garcilazo, Nucl. Phys. A 928 (2014), 73-88 doi:10.1016/j.nuclphysa.2014.02.019 [arXiv:1402.3171 [nucl-th]].
  * (34) W. Meguro, Y. R. Liu and M. Oka, Phys. Lett. B 704, 547-550 (2011) doi:10.1016/j.physletb.2011.09.088 [arXiv:1105.3693 [hep-ph]].
  * (35) Y. R. Liu and M. Oka, Phys. Rev. D 85, 014015 (2012) doi:10.1103/PhysRevD.85.014015 [arXiv:1103.4624 [hep-ph]].
  * (36) H. Huang, J. Ping and F. Wang, Phys. Rev. C 89, no.3, 035201 (2014) doi:10.1103/PhysRevC.89.035201 [arXiv:1311.4732 [hep-ph]].
  * (37) H. Huang, J. Ping and F. Wang, Phys. Rev. C 87, no.3, 034002 (2013) doi:10.1103/PhysRevC.87.034002
  * (38) Z. Liu, H. T. An, Z. W. Liu and X. Liu, Phys. Rev. D 105, no.3, 034006 (2022) doi:10.1103/PhysRevD.105.034006 [arXiv:2112.02510 [hep-ph]].
  * (39) Z. T. Lu, H. Y. Jiang and J. He, Phys. Rev. C 102, no.4, 045202 (2020) doi:10.1103/PhysRevC.102.045202 [arXiv:2007.15878 [nucl-th]].
  * (40) X. H. Chen, Q. N. Wang, W. Chen and H. X. Chen, Chin. Phys. C 45, no.4, 041002 (2021) doi:10.1088/1674-1137/abdfbe [arXiv:1906.11089 [hep-ph]].
  * (41) H. Mutuk and K. Azizi, Phys. Rev. D 105, no.9, 094021 (2022) doi:10.1103/PhysRevD.105.094021 [arXiv:2204.03050 [hep-ph]].
  * (42) N. Brambilla, H. X. Chen, A. Esposito, J. Ferretti, A. Francis, F. K. Guo, C. Hanhart, A. Hosaka, R. L. Jaffe and M. Karliner, et al. [arXiv:2203.16583 [hep-ph]].
  * (43) Q. B. Li and P. N. Shen, Eur. Phys. J. A 8, 417-421 (2000) doi:10.1007/s100530050052 [arXiv:nucl-th/9910060 [nucl-th]].
  * (44) X. H. Chen, Q. N. Wang, W. Chen and H. X. Chen, Phys. Rev. D 103 (2021) no.9, 094011 doi:10.1103/PhysRevD.103.094011 [arXiv:2103.09739 [hep-ph]].
  * (45) F. Etminan et al. [HAL QCD], Nucl. Phys. A 928 (2014), 89-98 doi:10.1016/j.nuclphysa.2014.05.014 [arXiv:1403.7284 [hep-lat]].
  * (46) S. Zhang and Y. G. Ma, Phys. Lett. B 811 (2020), 135867 doi:10.1016/j.physletb.2020.135867 [arXiv:2007.11170 [hep-ph]].
  * (47) C. J. Xiao, Y. B. Dong, T. Gutsche, V. E. Lyubovitskij and D. Y. Chen, Phys. Rev. D 101, 114032 (2020) doi:10.1103/PhysRevD.101.114032 [arXiv:2004.12415 [hep-ph]].
  * (48) J. Adam et al. [STAR], Phys. Lett. B 790, 490-497 (2019) doi:10.1016/j.physletb.2019.01.055 [arXiv:1808.02511 [hep-ex]].
  * (49) S. Acharya et al. [ALICE], Nature 588 (2020), 232-238 [erratum: Nature 590 (2021), E13] doi:10.1038/s41586-020-3001-6 [arXiv:2005.11495 [nucl-ex]].
  * (50) K. Aoki, H. Fujioka, T. Gogami, Y. Hidaka, E. Hiyama, R. Honda, A. Hosaka, Y. Ichikawa, M. Ieiri and M. Isaka, et al. [arXiv:2110.04462 [nucl-ex]].
  * (51) T. Nagae, The J-PARC project, Nucl. Phys. A 805, 486 (2008).
  * (52) F. Nerling [COMPASS Collaboration], Hadron Spectroscopy with COMPASS: Newest Results, EPJ Web Conf. 37 (2012) 01016.
  * (53) V. Obraztsov [OKA Collaboration], High statistics measurement of the $K^{+}\to\pi^{0}e^{+}\nu$(Ke3) decay formfactors, Nucl. Part. Phys. Proc. 273-275, 1330 (2016).
  * (54) B. Velghe [NA62-RK and NA48/2 Collaborations], $K^{\pm}\to\pi^{\pm}\gamma\gamma$ Studies at NA48/2 and NA62-RK Experiments at CERN, Nucl. Part. Phys. Proc. 273-275, 2720 (2016).
  * (55) J. Liu, Q. Wu, J. He, D. Y. Chen and T. Matsuki, Phys. Rev. D 101 (2020) no.1, 014003 doi:10.1103/PhysRevD.101.014003 [arXiv:2001.00212 [hep-ph]].
  * (56) J. Liu, D. Y. Chen and J. He, Eur. Phys. J. C 81 (2021) no.11, 965 doi:10.1140/epjc/s10052-021-09766-6 [arXiv:2108.00148 [hep-ph]].
  * (57) X. H. Liu, Q. Zhao and F. E. Close, Phys. Rev. D 77 (2008), 094005 doi:10.1103/PhysRevD.77.094005 [arXiv:0802.2648 [hep-ph]].
  * (58) Y. Huang, J. He, X. Liu, H. F. Zhang, J. J. Xie and X. R. Chen, Phys. Rev. D 93 (2016) no.3, 034022 doi:10.1103/PhysRevD.93.034022 [arXiv:1512.00981 [hep-ph]].
  * (59) D. Ronchen, M. Doring, F. Huang, H. Haberzettl, J. Haidenbauer, C. Hanhart, S. Krewald, U. G. Meissner and K. Nakayama, Eur. Phys. J. A 49, 44 (2013) doi:10.1140/epja/i2013-13044-5 [arXiv:1211.6998 [nucl-th]].
  * (60) J. He, Phys. Rev. D 95, no.7, 074031 (2017) doi:10.1103/PhysRevD.95.074031 [arXiv:1701.03738 [hep-ph]].
  * (61) R. Machleidt, K. Holinde and C. Elster, Phys. Rept. 149, 1-89 (1987) doi:10.1016/S0370-1573(87)80002-9
  * (62) A. Matsuyama, T. Sato and T. S. H. Lee, Phys. Rept. 439, 193-253 (2007) doi:10.1016/j.physrep.2006.12.003 [arXiv:nucl-th/0608051 [nucl-th]].
  * (63) W. Liu, C. M. Ko and Z. W. Lin, Phys. Rev. C 65, 015203 (2002) doi:10.1103/PhysRevC.65.015203
  * (64) C. Schutz, J. W. Durso, K. Holinde and J. Speth, Phys. Rev. C 49 (1994), 2671-2687 doi:10.1103/PhysRevC.49.2671
  * (65) B. Renner, Phys. Lett. B 33 (1970), 599-600 doi:10.1016/0370-2693(70)90359-X Copy to ClipboardDownload
  * (66) J. He and D. Y. Chen, Eur. Phys. J. C 79 (2019) no.11, 887 doi:10.1140/epjc/s10052-019-7419-7 [arXiv:1909.05681 [hep-ph]].
  * (67) D. Y. Chen, X. Liu and T. Matsuki, Phys. Rev. D 87 (2013) no.9, 094010 doi:10.1103/PhysRevD.87.094010 [arXiv:1304.0372 [hep-ph]].
  * (68) D. Y. Chen, X. Liu and T. Matsuki, Phys. Rev. D 90 (2014) no.3, 034019 doi:10.1103/PhysRevD.90.034019 [arXiv:1406.6763 [hep-ph]].
  * (69) S. H. Kim, A. Hosaka, H. C. Kim and H. Noumi, Phys. Rev. D 92, no.9, 094021 (2015) doi:10.1103/PhysRevD.92.094021 [arXiv:1509.03567 [hep-ph]].