# Drag coefficient of a rising bubble in a shear-thinning fluid using the
power-law scheme coupled with a Cahn-Hilliard equation with a variable
mobility: A lattice Boltzmann study and comparison with experiment

Amirabbas Ghorbanpour Arani Reza Haghani-Hassan-Abadi Mohammad Majidi
Mohammad-Hassan Rahimian<EMAIL_ADDRESS>School of Mechanical Engineering,
College of Engineering, University of Tehran, Tehran, Iran

###### Abstract

This study aims to investigate the behavior of multicomponent fluid flows
consisting of Newtonian and non-Newtonian components, especially terminal
velocity of a rising bubble in a power-law fluid. A recent lattice Boltzmann
(LB) model is extended using power-law scheme to be able to simulate both
Newtonian and non-Newtonian fluid flows at high density and viscosity ratios.
Also, a variable mobility is introduced in this study to minimize the
unphysical error around small bubbles in the domain. A three-component fluid
flow system is examined using a constant and variable mobility. It is shown
that each component has more stability using variable mobility while constant
mobility causes interface dissipation, leading to mass loss gradually. In
addition, two test cases including power-law fluid flows driven between two
parallel plates are conducted to show the accuracy and capability of the
model. To find a grid-independent computational domain, a grid independency
test is carried out to show that a $200\times 400$ domain size is suitable for
our computations. Then, terminal velocity of a rising bubble is compared to an
existing correlation in the literature, indicating that the results are in
good agreement with existing study so that average relative error in six
different cases is 5.66 $\%$. Also, the simulated examples show good
conformity to experimental results over a range of the Reynolds and Eötvös
numbers.

###### keywords:

non-Newtonian , phase-field , lattice Boltzmann method , Variable mobility ,
Terminal velocity

##  1 Introduction

Multiphase flows, especially rising bubble in Newtonian and non-Newtonian
liquids, are widely applicable in the industry. Multiphase flow systems exist
in many contexts such as water treatment processes, biomedical engineering,
pharmaceutical and food industries, and so on [1, 2]. Also, multiphase flow
systems are found in the droplet-based microfluidic devices which is
applicable in drug delivery and cell analysis [3, 4, 5]. Therefore, many
research efforts are targeted on predicting interaction between liquid and gas
not only in Newtonian liquids but also in non-Newtonian ones. Since gas phase
is in the form of bubble in many systems, it is of importance to find the
behaviors and characteristics of the bubble in Newtonian and non-Newtonian
liquids. As a matter of fact, understanding two-phase flow behavior can help
improve and maintain such systems. Terminal velocity of a bubble is an
important parameter in two-phase flow systems which plays an important role in
optimization and designing of gas-liquid systems [6], showing the necessity of
investigating bubble dynamics in Newtonian and non-Newtonian liquids. In
addition, terminal velocity can exhibit equilibrium in a system, as it occurs
when drag and buoyant force are balanced.

Numerical analysis has become very popular because of its prediction and
comprehension of phenomena in nature. Many numerical approaches have been
applied to analyze multiphase flows in macroscopic scale like front-tracking
[7], volume of fluid [8], level set method [9], etc. Battistella et al. [6]
used a front-tracking model to simulate rising bubble in a power-law fluid.
They used a drag closure proposed by dijkhuizen et al. [10] to obtain terminal
velocity of bubble, and compare it with numerical simulation results. In most
cases, the accuracy of the model was 20$\%$. Li et al. [11] experimentally
investigated rising bubble in a non-Newtonian liquid. They proposed a drag
coefficient correlation by considering the Reynolds number, aspect ratio, and
rheological properties of the liquid. Also, they plotted a map of shape
regimes for bubbles in a non-Newtonian fluid using Reynolds, Eötvös, and
Morton numbers. Islam et al. [12] used volume-of-fluid (VOF) approach to
investigate rising bubble behavior in a Newtonian and non-Newtonian liquid.
Behavior of a non-Newtonian liquid was described by power-law equation. They
studied terminal velocity and deformation shape of bubble using numerical
simulations. They showed that their results are in good agreement with
empirical correlations. Ohta et al. [13] used coupled level-set and volume-of-
fluid (CLSVOF) to simulate bubble rising in a non-Newtonian liquid. The
viscosity of liquid was determined by Carreau model. They investigated bubble
rise velocity and bubble shape in non-Newtonian fluids by changing Carreau
model parameters. They concluded that in case of implementing effective
parameters, Bhaga-Weber map can predict bubble shape in the non-Newtonian
fluid. Premlata et al. [14] carried out VOF method using the Carreau-Yasuda
model to simulate rising bubble in a non-Newtonian fluid. They validated their
model by comparing with results in the literature. They showed that rising
velocity increases by increasing shear-thinning tendency and deformation of
bubble increases at high Gallilei number and low Eötvös number. Pang and Lu
[15] used VOF method to simulate bubble rising in Carreau fluid. They
investigated the influence of surface tension, rheological index, and
different gravity field on bubble dynamics. They found that bubble deformation
increases by increasing Gallilei and Eötvös numbers and decreases with the
increase of rheological index. Also, they showed that terminal velocity
increases by increasing Gallilei number and decreasing rheological index.

In the recent decades, lattice Boltzmann method (LBM) has been used as an
alternative way to simulate fluids dynamics in mesoscopic scale. The LBM has
several advantages over traditional computational fluid dynamics (CFD) methods
such as ease of implementation, capability of scaling on parallel computers,
and handling complex geometries [16, 17, 18]. Several LB models have been
proposed for two-phase fluids such as color gradient [19], pseudopotential
[20], free energy [21], and mean-field [22]. All of these early models
confront with some drawbacks like limited density ratio and spurious currents
(non-physical velocity fields). Inamuro et al. [23] proposed a model based on
free energy that can capture large density ratio, but the Poisson equation in
the model reduced simplicity and efficiency of the LBM. Lee [24] proposed a
two-phase model using Cahn-Hilliard (CH) interface equation that can handle
large density and viscosity ratios. Also, spurious currents is in order of
$10^{-14}$ which are negligible. Haghani-Hassan-Abadi et al. [25] developed
Lee’s two-phase model to three-component model, implementing bulk free energy
and CH equation proposed in [26]. Their model maintained characteristics of
Lee’s binary model and can be used for one, two, and three component fluids.
This model was validated by several cases and capable of analyzing phenomena
in nature and industry [27]. Beside above models, many efforts have been done
to capture non-Newtonian behavior of fluids in LBM. Gabbanelli et al. [28]
extended LBM to non-Newtonian formulation using a truncated power-law model to
capture non-Newtonian characteristics of single-component fluids. Boyd et al.
[29] proposed a second order LB model to simulate non-Newtonian fluids. Their
model is efficient computationally and more accurate compared to previous
models. Moreover, further works were done for multiphase non-Newtonian fluids
[30, 31, 32].

To the authors’ best knowledge, the model introduced in [25] is strong due to
many numerical tests conducted by the authors. Therefore, in this study, this
model is used for the computation of multiphase flows. In this paper, we
enhance the model using variable mobility to minimize little numerical
dissipation around small bubbles reported in [27]. Also, Variable and constant
mobilities are compared to show how mass conservation is satisfied in the
present study. Variable mobility can effectively reduce dissolutions of bubble
and droplets to surrounding fluid [26, 33, 34]. Also, we couple the mentioned
model to the power-law scheme, using the truncated power-law model for
viscosity distribution, to simulate non-Newtonian fluids at high density and
viscosity ratios.

The paper is organized as follows. Section 2 represents lattice Boltzmann
equations (LBEs) together with CH equations and explanations about the
truncated power-law model. Benefits of variable mobility and validation of
method, which includes two non-Newtonian test cases studied by the present
method and compared with analytical solution to show accuracy and robustness
of the extension, are carried out in section 3. Then in section 4, grid
convergence test is carried out and rising behavior of a bubble in a shear-
thinning fluid is examined. Also, terminal velocity and deformation shape of
the bubble are compared with empirical results. Finally, this work is
concluded in section 5.

##  2 Methodology

###  2.1 Interface capturing model

Consider a domain $\varOmega$ of three incompressible, immiscible fluids. Each
fluid is distinguished by its phase-field variable ($C_{i}\quad i=1-3$) along
with the following constraint:

$\displaystyle C_{1}+C_{2}+C_{3}=1$ (1)

Free energy for three-component flows is as follows [26]:

$\displaystyle
F(C_{1},C_{2},C_{3})=\int_{\varOmega}\left[\frac{12}{\xi}E_{0}(C_{1},C_{2},C_{3})+\frac{3}{8}\xi\Sigma_{1}\left|\bm{\nabla}C_{1}\right|^{2}+\frac{3}{8}\xi\Sigma_{2}\left|\bm{\nabla}C_{2}\right|^{2}+\frac{3}{8}\xi\Sigma_{3}\left|\bm{\nabla}C_{3}\right|^{2}\right]d\varOmega$
(2)

where $\xi$ is the interface thickness which is usually of 3 to 5 lattices,
$\Sigma_{i}$ are capillary coefficients, and $E_{0}$ is the bulk free energy.
For consistency with diphasic systems, the capillary coefficients are defined
as $\Sigma_{i}=\sigma_{ij}+\sigma_{ik}-\sigma_{jk}\quad(i,j,k=1-3)$ where
$\sigma_{ij}$ is the surface tension coefficient between fluid $i$ and fluid
$j$. The bulk free energy is defined as

$\displaystyle
E_{0}=\sigma_{12}C_{1}^{2}C_{2}^{2}+\sigma_{13}C_{1}^{2}C_{3}^{2}+\sigma_{23}C_{2}^{2}C_{3}^{2}+C_{1}C_{2}C_{3}(\Sigma_{1}C_{1}+\Sigma_{2}C_{2}+\Sigma_{3}C_{3})+\varLambda
C_{1}^{2}C_{2}^{2}C_{3}^{2}$ (3)

where $\varLambda$ is a positive parameter for total spreading scenarios.

The three-component CH equations can be written as following:

$\displaystyle\dfrac{\partial C_{i}}{\partial
t}+\bm{\nabla}\cdot(\bm{u}C_{i})=\bm{\nabla}\cdot(D_{i}\bm{\nabla}\mu_{i})\quad
i=1-3$ (4)

$\displaystyle\mu_{i}=\frac{4\Sigma_{T}}{{\xi}}\sum_{j\neq
i}\dfrac{1}{\Sigma_{j}}[\partial_{i}E_{0}-\partial_{j}E_{0}]-\dfrac{3}{4}\xi\Sigma_{i}\nabla^{2}C_{i}\quad
i=1-3$ (5)

where $t$ is the time, $\bm{u}$ is the velocity, and $\mu_{i}$ is the chemical
potential of the fluid $i$. the mobility of the fluid $i$, $D_{i}$, and
$\Sigma_{T}$ are defined as

$D_{i}=\left\\{\begin{array}[]{rcl}\frac{D_{0}}{\Sigma_{i}}&Constant\
Mobility&\\\ \frac{D_{0}}{\Sigma_{i}}C_{i}(1-C_{i})&Variable\
Mobility\end{array}\right.\quad i=1-3$ (6)
$\displaystyle\frac{3}{\Sigma_{T}}=\frac{1}{\Sigma_{1}}+\frac{1}{\Sigma_{2}}+\frac{1}{\Sigma_{3}}$
(7)

where $D_{0}=\xi/300$ is an auxiliary parameter [25].

###  2.2 Lattice Boltzmann equations

The LBE with the single relaxation time can be written as [24]

$f_{\alpha}(\bm{x}+\bm{e}_{\alpha}\delta t,t+\delta
t)-f_{\alpha}(\bm{x},t)=-\dfrac{1}{\lambda}(f_{\alpha}-f_{\alpha}^{\rm
eq})+\dfrac{(\bm{e}_{\alpha}-\bm{u})\cdot\bm{F}f_{\alpha}^{\rm eq}}{\rho
c_{s}^{2}}$ (8)

where $f_{\alpha}$ is the particle distribution function in the
$\alpha$-direction, $\lambda$ is the relaxation time, $\rho$ is the density of
mixture, and $\bm{F}$ is the intermolecular force modeling non-ideal gas
effects. $f_{\alpha}^{\rm eq}$ is the equilibrium distribution function which
are defined as follows:

$f_{\alpha}^{\rm
eq}=w_{\alpha}\rho\left[1+\frac{(\bm{e}_{\alpha}\cdot\bm{u})}{c_{s}^{2}}+\frac{(\bm{e}_{\alpha}\cdot\bm{u})^{2}}{2c_{s}^{4}}-\frac{(\bm{u}\cdot\bm{u})}{2c_{s}^{2}}\right]$
(9)

In this study, the D2Q9 lattice structure is used with the following weight
coefficients, $w_{\alpha}$, and microscopic velocity, $\bm{e}_{\alpha}$, set:

$w_{\alpha}=\begin{cases}4/9,&\alpha=0\\\ 1/9,&\alpha=1-4\\\
1/36,&\alpha=5-8\end{cases}$ (10)
$\bm{e}_{\alpha}=c\begin{cases}(0,0),&\alpha=0\\\
(\cos[(\alpha-1)\pi/2],\sin[(\alpha-1)\pi/2]),&\alpha=1-4\\\
(\cos[(2\alpha-9)\pi/4],\sin[(2\alpha-9)\pi/4])\sqrt{2},&\alpha=5-8\end{cases}$
(11)

In the above equations, $c_{s}=\dfrac{c}{\sqrt{3}}$ is the speed of sound in
the D2Q9 lattice in which $c=\delta x/\delta t$ is the lattice speed with
$\delta x$ and $\delta t$ being the length scale and time scale, respectively.
The force term $\bm{F}$ is expressed as

$\bm{F}=\bm{\nabla}\rho
c_{s}^{2}-\bm{\nabla}p+\sum_{i=1}^{3}\mu_{i}\bm{\nabla}C_{i}+\bm{F}_{b}$ (12)

where $\bm{F}_{b}=-\rho g\overrightarrow{j}$ is the body force with $g$ being
the gravitational acceleration.

To obtain a LBE for recovering the hydrodynamics properties of the mixture,
the following transformation is used:

$g_{\alpha}=f_{\alpha}c_{s}^{2}+(p-\rho c_{s}^{2})\varGamma_{\alpha}(0)$ (13)

where $\varGamma_{\alpha}(\bm{u})=f_{\alpha}^{\rm eq}/\rho$. Taking the total
derivative of $g_{\alpha}$ and using the trapezoidal rule results in the
following LBE:

$\displaystyle\bar{g}_{\alpha}(\bm{x}+\bm{e}_{\alpha}\delta t,t+\delta
t)-\bar{g}_{\alpha}(\bm{x},t)$ (14)
$\displaystyle=-\frac{1}{\tau+0.5}(\bar{g}_{\alpha}-\bar{g}_{\alpha}^{\rm
eq})\rvert_{(\bm{x},t)}$ $\displaystyle+\delta
t(\bm{e}_{\alpha}-\bm{u})\cdot\left[\bm{\nabla}\rho
c_{s}^{2}(\varGamma_{\alpha}-\varGamma_{\alpha}(0))+(\sum_{i=1}^{3}\mu_{i}\bm{\nabla}C_{i}+\bm{F}_{b})\varGamma_{\alpha}\right]_{(\bm{x},t)}$

where the dimensionless relaxation time $\tau=\lambda/\delta t$ is related to
kinematic viscosity $\nu$ by $\tau=\nu/c_{s}^{2}$. In the above equations,
where $\bar{g}_{\alpha}$ and $\bar{g}_{\alpha}^{\rm eq}$ are modified particle
and equilibrium distribution functions, respectively, which are defined as:

$\begin{split}&\bar{g}_{\alpha}=g_{\alpha}+\frac{1}{2\tau}(g_{\alpha}-g_{\alpha}^{\rm
eq})-\frac{\delta t}{2}(\bm{e}_{\alpha}-\bm{u})\cdot\left[\bm{\nabla}\rho
c_{s}^{2}(\varGamma_{\alpha}-\varGamma_{\alpha}(0))+(\sum_{i=1}^{3}\mu_{i}\bm{\nabla}C_{i}+\bm{F}_{b})\varGamma_{\alpha}\right]\end{split}$
(15) $\displaystyle\bar{g}_{\alpha}^{\rm eq}=g_{\alpha}^{\rm eq}-\frac{\delta
t}{2}(\bm{e}_{\alpha}-\bm{u})\cdot\left[\bm{\nabla}\rho
c_{s}^{2}(\varGamma_{\alpha}-\varGamma_{\alpha}(0))+(\sum_{i=1}^{3}\mu_{i}\bm{\nabla}C_{i}+\bm{F}_{b})\varGamma_{\alpha}\right]$
(16)

In order to track the interfaces among three fluids, we use the following
distribution functions

$h_{\alpha}^{i}=\frac{C_{i}}{\rho}f_{\alpha}\quad i=1,2$ (17)

Note that $C_{3}$ is obtained via the constraint in Eq. (1).

Like previously, by taking the total derivative of $h_{\alpha}^{i}$ and using
the trapezoidal rule, we have the following LBE for interfaces capturing:

$\begin{split}&\bar{h}_{\alpha}^{i}(\bm{x}+\bm{e}_{\alpha}\delta t,t+\delta
t)-\bar{h}_{\alpha}^{i}(\bm{x},t)=-\frac{1}{\tau_{h}+0.5}(\bar{h}_{\alpha}^{i}-\bar{h}_{\alpha}^{i,\rm
eq})\rvert_{(\bm{x},t)}\\\ &+\delta
t(\bm{e}_{\alpha}-\bm{u})\cdot\left[\bm{\nabla}C_{i}-\frac{C_{i}}{\rho
c_{s}^{2}}(\bm{\nabla}p-\sum_{i=1}^{3}\mu_{i}\bm{\nabla}C_{i}-\bm{F}_{b})\right]\varGamma_{\alpha}\rvert_{(\bm{x},t)}\\\
&+\frac{\delta
t}{2}D_{i}\nabla^{2}\mu_{i}\varGamma_{\alpha}\rvert_{(\bm{x},t)}+\frac{\delta
t}{2}D_{i}\nabla^{2}\mu_{i}\varGamma_{\alpha}\rvert_{(\bm{x}+\bm{e}_{\alpha}\delta
t,t)}\\\ &+\frac{\delta
t}{2}\frac{D_{i}}{C_{i}(1-C_{i})}\bm{\nabla}\mu_{i}\varGamma_{\alpha}[\bm{\nabla}C_{i}-\bm{\nabla}C_{i}^{2}]\rvert_{(\bm{x},t)}+\frac{\delta
t}{2}\frac{D_{i}}{C_{i}(1-C_{i})}\bm{\nabla}{\mu}_{i}\varGamma_{\alpha}[\bm{\nabla}C_{i}-\bm{\nabla}C_{i}^{2}]\rvert_{(\bm{x}+\bm{e}_{\alpha}\delta
t,t)}\\\ &\quad i=1,2\end{split}$ (18)

where $\tau_{h}=0.5$ is the dimensionless relaxation time.
$\bar{h}_{\alpha}^{i}$ and $\bar{h}_{\alpha}^{i,\rm eq}$ are the modified
particle and equilibrium distribution functions, respectively, defined as

$\displaystyle\bar{h}_{\alpha}^{i}=h_{\alpha}^{i}+\frac{1}{2\tau}(h_{\alpha}^{i}-h_{\alpha}^{i,\rm
eq})-\frac{\delta
t}{2}(\bm{e}_{\alpha}-\bm{u})\cdot\left[\bm{\nabla}C_{i}-\frac{C_{i}}{\rho
c_{s}^{2}}(\bm{\nabla}p-\sum_{i=1}^{3}\mu_{i}\bm{\nabla}C_{i}-\bm{F}_{b})\right]\varGamma_{\alpha}$
(19)

$\begin{split}&\bar{h}_{\alpha}^{i,\rm eq}=h_{\alpha}^{i,\rm eq}-\frac{\delta
t}{2}(\bm{e}_{\alpha}-\bm{u})\cdot\left[\bm{\nabla}C_{i}-\frac{C_{i}}{\rho
c_{s}^{2}}(\bm{\nabla}p-\sum_{i=1}^{3}\mu_{i}\bm{\nabla}C_{i}-\bm{F}_{b})\right]\varGamma_{\alpha}\end{split}$
(20)

where $h_{\alpha}^{i,\rm eq}=C_{i}\varGamma_{\alpha}$.

As mentioned earlier, in this study we employ a variable mobility as it shows
better mass conservation properties. Equation (6) defines both constant and
variable mobilities. As can be seen, the variable mobility is not constant and
dependent of $C_{i}$. As such, when taking the total derivative of Eq. (17)
and considering the variable mobilities, the diffusion term on the right-hand
side of Eq. (4) simplifies to

$\begin{split}&\bm{\nabla}\cdot(D_{i}\bm{\nabla}\mu_{i})=\bm{\nabla}\cdot[\frac{D_{0}}{\Sigma_{i}}C_{i}(1-C_{i})\bm{\nabla}\mu_{i}]\\\
&=\frac{D_{0}}{\Sigma_{i}}C_{i}(1-C_{i})\nabla^{2}\mu_{i}+\frac{D_{0}}{\Sigma_{i}}(\bm{\nabla}C_{i}-\bm{\nabla}C_{i}^{2})\bm{\nabla}\mu_{i}\\\
&=D_{i}\nabla^{2}\mu_{i}+\frac{D_{i}}{C_{i}(1-C_{i})}(\bm{\nabla}C_{i}-\bm{\nabla}C_{i}^{2})\bm{\nabla}\mu_{i}\quad
i=1,2\end{split}$ (21)

which the last term is related to the variable mobilities and also appears in
the interface capturing LBE, i.e., Eq. (18).

In the above equations, second-order central and mixed finite differences and
directional derivatives of a scaler variable are implemented to enhance the
numerical stability [27].

The phase-field variables, momentum and pressure can be calculated by taking
the moments of the modified distribution functions

$C_{i}=\sum_{\alpha}\bar{h}_{\alpha}^{i}\quad i=1,2$ (22)
$\rho\bm{u}=\frac{1}{c_{s}^{2}}\sum_{\alpha}\bm{e}_{\alpha}\bar{g}_{\alpha}+\frac{\delta
t}{2}(\sum_{i=1}^{3}\mu_{i}\bm{\nabla}C_{i}+\bm{F}_{b})$ (23)
$p=\sum_{\alpha}\bar{g}_{\alpha}+\bm{u}\cdot\bm{\nabla}\rho c_{s}^{2}$ (24)

###  2.3 Strain rate tensor and relaxation time calculation

In the continuum fluid flows, the Navier-stokes (NS) equations are recovered
from LBE through the Chapman-Enskog (CE) analysis. The strain rate tensor is
calculated by the second moment of non-equilibrium distribution function as
follows [35, 36]:

$\displaystyle
D=\frac{1}{2}(\bm{\nabla}\bm{u}+\nabla{\bm{u}}^{t})=-\frac{\Pi^{(1)}}{2(\tau+0.5)\rho{c_{s}}^{2}}$
(25)

where $\bm{u}^{t}$ denotes the transposed matrix to $\bm{u}$ and $\Pi^{(1)}$
as first-order momentum flux for the D2Q9 lattice is given as

$\Pi^{(1)}_{\alpha\beta}(\bm{x},t)=\sum_{i=0}^{8}f^{(1)}_{i}(\bm{x},t)e_{\alpha
i}e_{\beta i}$ (26)

where $f_{i}^{(1)}$ is a first-order perturbation distribution function [35].
Using the transformation in Eq. (13), the relationship between non-equilibrium
particle distribution functions can be obtained as follows:

$g_{\alpha}-g_{\alpha}^{\rm eq}=(f_{\alpha}-f_{\alpha}^{\rm eq})c^{2}_{s}$
(27)

By using Eqs. (15) and (16) we have

$g_{\alpha}-g_{\alpha}^{\rm
eq}=\frac{2\tau}{2\tau+1}(\bar{g}_{\alpha}-\bar{g}_{\alpha}^{\rm eq})$ (28)

As such, the non-equilibrium particle distribution function of $f$ can be
replaced as follows:

$f_{\alpha}-f_{\alpha}^{\rm
eq}=\frac{1}{c^{2}_{s}}\frac{2\tau}{2\tau+1}(\bar{g}_{\alpha}-\bar{g}_{\alpha}^{\rm
eq})$ (29)

As a result, knowing the magnitude of non-equilibrium distribution function of
$f$, $f_{\alpha}-f_{\alpha}^{\rm eq}$, Eq. (26) and Eq. (25) can be solved at
each time step to obtain strain rate tensor, i.e., $D$.

Constitutive equation for a generalized Newtonian fluid is defined as [37]

$T=2\eta(\dot{\gamma})D$ (30)

where $\eta(\dot{\gamma})$ and $T$ are the apparent viscosity function of
shear rate and shear stress tensor, respectively. The non-Newtonian fluid in
this study is modeled as a power-law fluid with consistency $k$ and index $n$.
Therefore, the apparent viscosity in Eq. (30) is defined as follows:

$\eta(\dot{\gamma})=k\dot{\gamma}^{n-1}$ (31)

The shear rate, as an invariant of strain rate tensor, is expressed as
$\dot{\gamma}=\sqrt{2D:D}$, and can be calculated locally by the non-
equilibrium particle distribution function using Eq. (25).

In the non-Newtonian fluids, viscosity is not constant, so relaxation
parameter is calculated at each lattice site by recalling the relationship
between fluid viscosity and relaxation parameter:

$\tau=\frac{\eta({\dot{\gamma}}(\tau))}{\rho c^{2}_{s}}$ (32)

Therefore, the relaxation parameter obtained by Eq. (32) is used in
Eq.(14).The above equation is an implicit equation which can be solved by
iterative method numerically. However, in this study, $\dot{\gamma}(\tau)$ is
updated by the relaxation parameter at the previous time step [35, 38]. Such
approach uses much less computational cost.

Apparent viscosity corresponds to shear-thinning fluid, shear-thickening
fluid, and Newtonian fluid when $n$ is less than 1, more than 1, and equal to
1, respectively. Therefore, when shear rate increases, the viscosity of shear-
thinning fluid decreases. As such, the shear rate in this study is computed by
particle distribution function rather than using derivatives of velocities by
finite difference scheme. This approach, deriving parameters locally, is more
consistent with the philosophy of the LBM [35].

Direct implementation of a power-law fluid flow in the LBM confronts an
important problem. In a shear-thinning fluid ($n<1$) and at zero shear rates
($\dot{\gamma}=0$), the relaxation parameter diverges. On the other hand, in a
shear-thickening fluid ($n>1$) and at zero shear rates ($\dot{\gamma}=0$), the
relaxation parameter becomes zero, leading to instability in solution as well.
Hence, in this study the truncated power-law viscosity model is implemented to
limit the range of viscosity as follows [28]:

$\eta(\dot{\gamma})=\left\\{\begin{array}[]{rcl}k\dot{\gamma}^{n-1}_{0}&\mbox{for}&\dot{\gamma}<\dot{\gamma}_{0}\\\
k\dot{\gamma}^{n-1}&\mbox{for}&\dot{\gamma}_{0}<\dot{\gamma}<\dot{\gamma}_{\infty}\\\
k\dot{\gamma}^{n-1}_{\infty}&\mbox{for}&\dot{\gamma}_{\infty}<\dot{\gamma}\end{array}\right.$
(33)

Therefore, in the truncated power-law model, range of relaxation parameter is
limited and prevents solution from instability.

##  3 Validation

###  3.1 Effect of variable mobility

In this section, we examine the effect of mobility in our simulations. We set
our variable and constant mobility as in Eq. (6). As can be seen, the
$C_{i}(1-C_{i})$ term, which was applied in [39], is multiplied by the
constant mobility introduced in [25]. In fact, instead of a constant mobility
applied to the whole domain, a variable mobility is implemented on the
interfaces only. Figure 1 shows initial condition of three components so that
red, green, and blue colors correspond to component 1, 2, and 3, respectively.
Green bubble rises because of buoyant force to stick to the red droplet. This
situation, bubble-droplet rising, is a critical situation in our simulation
which can cause mass loss. Table 1 shows properties of the three-component
fluid flow system in a $200\times 400$ computational domain we investigate.

Figure 1: Initial condition of three-component fluid flows. $\rho_{1}$ | $\rho_{2}$ | $\rho_{3}$ | g | $\frac{\tau_{1}}{\tau_{3}}=\frac{\tau_{2}}{\tau_{3}}$ | $\sigma_{12}=\sigma_{13}=\sigma_{23}$
---|---|---|---|---|---
$1$ | $0.01$ | $1$ | $10^{-6}$ | $2$ | $10^{-4}$
Table 1: Physical properties of three-component fluid flow system

Seven line contours at different times were plotted to compare our simulation
using variable mobility and constant mobility. As shown in Fig. 2, red lines
and black lines represent level contour of 1 of components obtained by
constant mobility and variable mobility, respectively. The droplet, which is
under no buoyant force, preserves its mass with variable mobility while with
constant mobility its mass is dissipated. Therefore, it can be seen in the
Fig. 2 that red line droplet shrink inside of the black line droplet over
time. Also, the rising bubble loses its mass using constant mobility, so the
black line bubble get more buoyant force, indicating that it deforms sooner
and rises faster than red line bubble.

To compare our results precisely, concentration of each component is recorded
when it is above 0.5. As Fig. 3 shows, summation of each $C_{i}$ in all of the
lattices demonstrates that green bubble and red droplet gradually lose their
mass, and the losing mass adds to the blue component mass. When the number of
lattices are increased, transformations of mass between components decrease
with constant mobility while variable mobility can omit this transformation
even in a small domain. Since variable mobility conserves mass exceedingly, we
perform our simulations using variable mobility in the present study.

Figure 2: Line contour of bubble and droplet with constant and variable
mobility Figure 3: Variation of total concentration of each component over
time.

###  3.2 Single phase power-law fluid flows between two parallel plates

In this subsection and the following one, two benchmarks are conducted to
validate the implementation of our code. First benchmark is a pressure-driven
flow between two fixed parallel plates at a distance of $2L$ in the
$y$-direction while fluid flow moves in the $x$-direction as shown in Fig. 4.

Figure 4: Schematic of a single phase power-law fluid flow between two
parallel plates.

The analytical solution for steady state velocity profile with power-law index
is as follows: [6]

$u_{x}(y)=u_{in}\frac{2n+1}{n+1}(1-{\lvert\frac{y}{L}\rvert}^{\frac{n+1}{n}})$
(34)

where $u_{in}$ is the inlet velocity, $y$ is the distance from the center of
two parallel plates, and $n$ is the power-law index. The simulations are
performed for various power-law indices including shear-thickening, shear-
thinning, and Newtonian fluids. The No-slip boundary condition is imposed on
the top and bottom walls. Velocity inlet [40] and convective boundary
conditions [41] are applied on the left and right, respectively. The
simulations are performed in a two dimensional rectangular domain where width
of the channel is 60 and the length is 240 lattice unit. Initially,
compositions of $C_{1}$ and $C_{2}$ are equal to zero so that the channel is
filled completely by the third component. A constant velocity in the
$x$-direction ($u_{in}=0.1$) is imposed at the left boundary. Several
simulations by various power-law indices are carried out and theoretical and
numerical velocity profiles are compared at steady state as shown in Fig. 5.

Figure 5: Comparison of analytical and LBM velocity profiles in a 2D channel
for power-law fluids with indices of $n=0.5,0.8,1,1.5,1.8$.

The relative $L_{2}$-norm of the velocity has been calculated at steady state
condition between analytical and numerical results as presented in the
following equation:

$E_{rel}=\frac{\lvert\lvert u_{x}-u^{analytical}_{x}\rvert\rvert}{\lvert\lvert
u^{analytical}_{x}\rvert\rvert}$ (35)

where $u_{x}$ is the $x$ component of velocity, $u^{analytical}_{x}$ is the
$x$ component of analytical velocity, and $E_{rel}$ is the relative
$L_{2}$-norm of the velocity.

Table 2 shows the relative errors of the power-law fluid velocity between two
parallel plates at different $n$. It can be concluded that the simulation
results are in good agreement with the analytical solution.

n | 0.5 | 0.8 | 1 | 1.5 | 1.8
---|---|---|---|---|---
$E_{rel}$ | $0.03$ | $0.02$ | $0.01$ | $0.01$ | $0.01$
Table 2: Relative errors of the velocity at different power-law indices.

###  3.3 Two-phase pressure driven power-law fluid flows between two parallel
plates

In this subsection, two-phase layered flow is considered between two parallel
plates which is driven by pressure gradient. This benchmark is used to assess
the ability of the present LB model in phases interaction. So, one of the
component is set to zero initially to obtain a two-phase flow. As illustrated
schematically in Fig. 6, each component fills half of the 2D channel and the
composition 2 is set to zero. The height of the channel is $2H$ and the
component 3 is placed at $\lvert y\rvert<Y_{i}$. Periodic boundary condition
is applied on the left and right boundaries while bounce-back scheme is
imposed on the top and bottom. Each fluid is driven by a constant body force
($-\frac{\partial P}{\partial x}=F_{b}$) whose value is $1.5\times 10^{-8}$.
The simulation is performed in a domain with $5\times 200$ lattice units.
Density ratio is equal to unity but different viscosity ratios are
investigated. The viscosity ratio is defined as $\lambda={\mu_{3}}/{\mu_{1}}$
which is used for Newtonian systems. For non-Newtonian systems new parameter
is defined as $\lambda^{p}=\mu^{p}_{3}/\mu^{p}_{1}$ which is used instead of
$\lambda$ [31]. Truncated power-law model is applied in the simulation and the
relaxation parameter limit, $\tau_{min}$ and $\tau_{max}$, are set. The
surface tension is $\sigma_{13}=5\times 10^{-4}$ and the interface thickness
is $\xi=4$. The simulation is run until it reaches the steady state which is
determined by $\frac{{\lvert\lvert u(x,t)-u(x,t-\delta
t)\rvert\rvert}_{2}}{{\lvert\lvert u(x,t)\rvert\rvert}_{2}}<10^{-8}$. As
discussed in [31], the analytical velocity in the $x$-direction is obtained as
follows:

$u^{*}_{x}(y)=\left\\{\begin{array}[]{rcl}\frac{n_{R}}{n_{R}+1}(-\frac{1}{\mu^{R}_{P}}\frac{\partial
P}{\partial x})^{\frac{1}{n_{R}}}(Y^{\frac{1+n_{R}}{n_{R}}}_{i}-{\lvert
y\rvert}^{\frac{1+n_{R}}{n_{R}}})+\frac{n_{B}}{n_{B}+1}(-\frac{1}{\mu^{B}_{P}}\frac{\partial
P}{\partial x})^{\frac{1}{n_{B}}}(Y^{\frac{1+n_{B}}{n_{B}}}_{i}-{\lvert
y\rvert}^{\frac{1+n_{B}}{n_{B}}})&\mbox{if}&\lvert y\rvert<Y_{i}\\\
\frac{n_{B}}{n_{B}+1}(-\frac{1}{\mu^{B}_{P}}\frac{\partial P}{\partial
x})^{\frac{1}{n_{B}}}(Y^{\frac{1+n_{B}}{n_{B}}}_{i}-{\lvert
y\rvert}^{\frac{1+n_{B}}{n_{B}}})&\mbox{otherwise}&\end{array}\right.$ (36)

It should be noted that the analytical maximum velocity is $u^{*}_{x}(y=0)$.
The maximum velocity should not be greater than $10^{-2}$ because of
compressibility error [31].

Figure 6: Schematic of the two-phase layered fluid flow between two parallel
plates. Component 3 is placed at the center, component 2 is equal to zero, and
component 1 is placed at top and bottom.

Several simulations with different $\lambda^{P}$ and power-law indices are
performed. Analytical and numerical velocity profiles are compared as shown in
Fig. 7. The relative $L_{2}$-norm which was presented in Eq. (35) is obtained
in each simulation to investigate the accuracy of this model. Table 3 shows
the relative errors of the simulated and analytical velocity. As it can be
seen in the table 3, $\lambda^{P}$ and $\lambda$ are used for non-Newtonian
and Newtonian systems, respectively.

Figure 7: Comparison of analytical and LBM velocity profiles in a 2D channel for two-phase power-law layered flow at different $\lambda^{P}$ and power-law indices. (a) $n_{3}=0.5$, $n_{1}=1$, $\lambda^{P}=0.002$, (b) $n_{3}=0.5$, $n_{1}=1$, $\lambda^{P}=0.02$, (c) $n_{3}=0.5$, $n_{1}=1$, $\lambda^{P}=0.004$, (d) $n_{3}=1$, $n_{1}=1$, $\lambda=0.3$, (e) $n_{3}=1$, $n_{1}=1$, $\lambda=1$, (f) $n_{3}=1.5$ ,$n_{1}=1$, $\lambda^{P}=2000$, (g) $n_{3}=1.5$, $n_{1}=1$, $\lambda^{P}=5000$, and (h) $n_{3}=1.5$, $n_{1}=1$, $\lambda^{P}=10000$.  $\lambda^{P}$ | $\lambda$ | $n_{3}$ | $n_{1}$ | $E_{rel}$
---|---|---|---|---
$0.002$ |  | $0.5$ | $1$ | $0.04$
$0.02$ |  | $0.5$ | $1$ | $0.01$
$0.004$ |  | $0.5$ | $1$ | $0.01$
| $0.3$ | $0.5$ | $1$ | $0.01$
| $1$ | $1$ | $1$ | $0.01$
$2000$ |  | $1$ | $1$ | $0.01$
$5000$ |  | $1.5$ | $1$ | $0.01$
$10000$ |  | $1.5$ | $1$ | $0.01$
Table 3: Relative errors of the velocity for different power-law indices and
different $\lambda^{P}$ and $\lambda$ in the two-phase power-law layered flow
driven by pressure gradient.

The relative errors are in order of $10^{-2}$ in each simulation which shows
that the numerical results are in good agreement with the analytical solution.
It should be noted that as long as the maximum velocity decreases, the
compressibility error will be also smaller, therefore, the relative errors are
lower.

##  4 Results

###  4.1 Drag coefficient and dimensionless numbers of a bubble rising in a
liquid

In this section, the terminal velocity of a bubble rising in a power-law fluid
is investigated. When the bubble rising (subscript $g$) in a liquid (subscript
$l$) reaches its terminal velocity, buoyancy and drag forces exerting on the
bubble become equal. Therefore, the drag coefficient can be written as [10,
11]:

$C_{D}=\frac{4dg(\rho_{l}-\rho_{g})}{3\rho_{l}U^{2}_{t}}$ (37)

The Reynolds number, Morton number, and Eötvös number in power-law fluids are
defined as [11, 42]

$Re=\frac{\rho_{l}U^{2-n}d^{n}_{e}}{k}$ (38)

$Mo=\frac{g\Delta\rho}{\rho^{2}_{l}\sigma^{3}}[k{(\frac{U}{d_{e}})}^{n-1}]^{4}$
(39) $Eo=\frac{(\rho_{l}-\rho_{g})gd^{2}_{e}}{\sigma}$ (40)

Where $d_{e}$ is the equivalent diameter which is defined as
$d_{e}=\sqrt[3]{\frac{6V_{b}}{\pi}}$ using bubble total volume, $V_{b}$. In
the Newtonian fluids, the drag coefficient is only a function of Reynolds
number, but in the power-law fluids the index $n$ has also effect on the drag
coefficient ($C_{D}=f(Re,n)$)[42].

Li et al. [11] proposed a drag correlation for a bubble rising in a shear-
thinning liquid, considering deformation of bubble and rheological properties
of fluids besides Reynolds number. They showed that the proposed correlation
is capable of handling spherical and deformed bubble.

$C_{D}=\frac{16}{Re}(1+0.43Re^{0.44})(1+3.868n^{0.7528}(1-E)^{0.6810}).$ (41)

where $E$ is the aspect ratio of the bubble, expressing the ratio of vertical
and horizontal diameters of the bubble. They showed that this correlation have
great agreement with experimental results under the finite range of
$0.05<Re<300$. It should be noted that in our study, a bubble rising by
buoyant force, the bubble experiences $E\geq 1$ because the only force exerted
on the bubble is from bubble rear. But if a bubble rises in other non-
Newtonian liquids applying forces from both left and right like viscoelastic
fluids, bubble aspect ratio can be less than unity.

###  4.2 Grid independency

Grid independency test in four domains is performed to find the grid
independent computational domain. A domain of $L\times 2L$ is considered which
$L$ is the length of it. The no-slip boundary condition is imposed on the top
and bottom and periodic boundary condition is applied on the left and right
boundaries. Time evolution of gravity center of the single bubble, which is
under buoyant force, is plotted for each grid in Fig. 8. Table 4 shows
rheological properties of liquid and physical properties of gas-liquid system,
where $\tau_{g}$, $d_{g}$, and $g$ represent relaxation time of the disperesed
component, the rising bubble, diameter of the bubble, and gravitational
acceleration, respectively.

Figure 8: Gravity center of the bubble as function of time for different grid resolution.  Grid number | Domain Size | $\rho_{l}/\rho_{g}$ | $k$ | $n$ | $\tau_{g}$ | g | $d_{g}$ | $\sigma$
---|---|---|---|---|---|---|---|---
$1$ | $50\times 100$ | $800$ | $2\times 10^{-4}$ | $0.5$ | $0.5$ | $8\times 10^{-6}$ | $\frac{2L}{5}$ | $1.044\times 10^{-4}$
$2$ | $100\times 200$ | $800$ | $10^{-4}$ | $0.5$ | $0.5$ | $10^{-6}$ | $\frac{2L}{5}$ | $5.22\times 10^{-5}$
$3$ | $200\times 400$ | $800$ | $5\times 10^{-5}$ | $0.5$ | $0.5$ | $1.25\times 10^{-7}$ | $\frac{2L}{5}$ | $2.61\times 10^{-5}$
$4$ | $400\times 800$ | $800$ | $2.5\times 10^{-5}$ | $0.5$ | $0.5$ | $1.56\times 10^{-8}$ | $\frac{2L}{5}$ | $1.3\times 10^{-5}$
Table 4: Rheological and physical properties of different computational
domain.

The parameters used in the Table 4 are selected so that dimensionless numbers,
Reynolds, Eötvös, and Morton numbers, remain constant in each grid. The gas-
liquid system is assumed to be quiescent initially and initial condition for
concentration is as follows:

$\phi_{d}(\bm{x},0)=0.5-0.5\tanh(\frac{\sqrt{(x-0.5L)^{2}+(y-0.6L)^{2}}-0.5{d_{b}}}{0.5{\xi}})$
(42)

In this study, the velocity of the bubble and gravity center of the bubble are
recorded in average-form as follows:

$\bar{V}=\frac{\sum C_{g}V}{\sum C_{g}}$ (43) $\bar{y}=\frac{\sum C_{g}y}{\sum
C_{g}}$ (44)

Also, we define dimensionless form for velocity and location for better
comparison by dividing them to $\sqrt{gd_{g}}$ and $2L$, respectively. To
choose a suitable domain, dimensionless velocity and center of gravity of the
bubble are examined. Magnitude of dimensionless terminal velocity of the
bubbles are $0.198$, $0.243$, $0.268$, and $0.273$ for grid 1, grid 2, grid 3,
and grid 4, respectively. Since the error between grid 3 and grid 4 is less
than $1\%$, grid 3 is suitable for the rest of our simulations. Also, curves
correspond to the grid 3 and 4 are so close as shown in Fig. 8

###  4.3  Bubble rising in a shear-thinning liquid

In this section, the terminal velocity of a bubble in a shear-thinning liquid,
whose viscosity decreases by applying shear stress, is investigated. To
investigate prediction of correlation proposed by Li et al. [11], some
numerical simulation with different rheological properties are conducted using
present LB model. All the simulations represent shear-thinning behavior of the
liquid ($n<1$). It should be noted that the liquid contains solely a bubble
and is completely pure. Properties of gas-liquid systems were provided in
table 5

Case | $\rho_{c}/\rho_{d}$ | $n$ | $k$ | $g$ | $\tau_{d}$ | $d_{b}$ | $\sigma$
---|---|---|---|---|---|---|---
$1$ | $800$ | $0.5$ | $1\times 10^{-4}$ | $5\times 10^{-7}$ | $0.5$ | $\frac{L}{4}$ | $2.61\times 10^{-6}$
$2$ | $800$ | $0.5$ | $1\times 10^{-4}$ | $5\times 10^{-7}$ | $0.5$ | $\frac{L}{4}$ | $2.61\times 10^{-5}$
$3$ | $800$ | $0.5$ | $2\times 10^{-4}$ | $5\times 10^{-7}$ | $0.5$ | $\frac{L}{4}$ | $5\times 10^{-4}$
$4$ | $800$ | $0.6$ | $5\times 10^{-4}$ | $2\times 10^{-7}$ | $0.5$ | $\frac{L}{4}$ | $2.61\times 10^{-5}$
$5$ | $800$ | $0.63$ | $5\times 10^{-4}$ | $2\times 10^{-7}$ | $0.5$ | $\frac{L}{4}$ | $2.61\times 10^{-5}$
$6$ | $800$ | $0.8$ | $1\times 10^{-2}$ | $5\times 10^{-7}$ | $0.5$ | $\frac{L}{4}$ | $2.61\times 10^{-5}$
Table 5: Physical properties of the gas-liquid systems.

Since the behavior of the bubble is determined using dimensionless numbers,
the Reynolds and Eötvös numbers for each case are shown in table 6:

Case | $Re$ | $Eo$
---|---|---
$1$ | $4.05$ | $478.9$
$2$ | $4.29$ | $47.89$
$3$ | $1.36$ | $2.5$
$4$ | $0.44$ | $19.15$
$5$ | $0.99$ | $19.15$
$6$ | $0.47$ | $47.89$
Table 6: Dimensionless numbers at different cases.

Velocity of the bubble increases gradually and reaches a constant value
(steady state). To show this process, for instance, in case 4, lattice unit
and dimensionless velocity of the bubble are shown in Fig. 9 (a) and (b),
respectively.

Figure 9: Dimensionless (a) and lattice unit (b) velocity of the bubble in the
case 4

As drag coefficient is dependent on the Reynolds number, the Reynolds number
plays an important role in determining the location of the bubble. The effect
of Reynolds number is shown in Fig. 10 by recording the gravity center of the
bubble, which was made dimensionless by the height of the domain, for all the
cases. Initially, the bubble starts to move gradually, demonstrating that the
slope of lines in all the cases is equal to zero. Then, the slope of line
increases to reach a constant value. At higher Reynolds number, the slope of
curves and gravity center are higher, indicating that the bubble associated
with the higher Reynolds number goes upward faster.

Figure 10: Dimensionless gravity center of the bubble in six cases.

To compare our results with the mentioned correlation, a relative error,
$(|numerical-correlation|/correlation)$, for each simulation is provided in
Table 7 with the maximum relative error of 19$\%$, minimum relative error of
0, and average relative error of 5.6$\%$. Table 7 shows that terminal velocity
of the bubble at different conditions are in good agreement with empirical
correlation proposed by Li et al. [11]. One of the reasons behind small
deviation between the results is that in the present model, the shear-thinning
liquid is assumed to be pure, unlike experimental results which were done
using carboxymethyl cellulose (CMC) solutions.

Case | $n$ | $E$ |  $U_{t}$ (Li et. al [11]) |  $U_{t}$ (present) | Relative Error
---|---|---|---|---|---
$1$ | $0.5$ | $0.77$ | $0.00159$ | $0.00156$ | $0.01$
$2$ | $0.5$ | $0.73$ | $0.00159$ | $0.00157$ | $0.01$
$3$ | $0.5$ | $0.98$ | $0.0012$ | $0.0011$ | $0.08$
$4$ | $0.6$ | $0.96$ | $0.00046$ | $0.00046$ | $0$
$5$ | $0.63$ | $0.9$ | $0.0006$ | $0.00063$ | $0.05$
$6$ | $0.8$ | $0.94$ | $0.00071$ | $0.00085$ | $0.19$
Table 7: Comparison of terminal velocity obtained by present LB model and
correlation.

Different shapes of bubble are obtained at different Reynolds and Eötvös
numbers. Since surface tension tries to maintain minimum surface area of the
bubble, higher Eötvös number represents more deformation tendency of the
bubble. Also, higher Reynolds number represents higher momentum force in the
liquid, resulting in higher exerted force on the bubble. To analyze better,
concentration, viscosity, and velocity contours are plotted at steady state
condition for each case.

In the case 1, the surface tension is very low; therefore, the bubble tends to
deform more as shown in Fig. 11 (a). Kinematic viscosity distribution of the
liquid shows that when the bubble rises, viscosity decreases around the bubble
because of the shear-thinning behavior. As bubble rises, shear rate around the
bubble increases, leading to decrease in viscosity. The farther away from the
bubble, the higher the viscosity. When the bubble gets far from the bottom
boundary, the viscosity in the vicinity of the boundary increases or when the
bubble gets close to the top boundary, the viscosity nearby the top boundary
becomes lower gradually which indicate shear-thinning behavior of the liquid.
Also, left and right boundaries are periodic, which indicate less shear rate
in the left and right of the bubble, causing a thin line of high viscosity in
the left and right of the bubble.

Figure 11: Final (a) concentration, (b) kinematic viscosity, and (b) velocity
contours of case 1.

In the case 2, the physical properties are the same as the case 1 but the
surface tension. Therefore, the final shape of the bubble in the case 2, which
has more surface tension, is less deformed with respect to the final shape of
the bubble in the case 1 as shown in Fig. 12 (a). Also, Fig. 12 (b)
demonstrates that by increasing surface tension, the kinematic viscosity
around the bubble increases because a low deformation indicates a low shear
rate.

Figure 12: Final (a) concentration, (b) kinematic viscosity, and (b) velocity
contours of case 2.

In the case 3, Reynolds and Eötvös numbers decrease by increasing consistency
index and surface tension coefficient, respectively. As a result, the bubble
is capable of maintaining spherical shape over time as shown in Fig. 13 (a).
Also, Fig. 13 (b) shows that kinematic viscosity distribution is higher with
respect to the cases 1 and 2. Also, the two regions on the left and right of
the bubble, having large viscosity, is bigger with respect to the cases 1 and
2.

Figure 13: Final (a) concentration, (b) kinematic viscosity, and (b) velocity
contours of case 3.

In the cases 4 and 5, higher power-law indices, $n=0.6$ and $n=0.63$, are
implemented while all other parameters are the same. Since in our simulation
the shear rate magnitude is less than unity, shear thinning behavior increases
by increasing power-law index [43]. Therefore, as Fig. 14 (b) and Fig. 15 (b)
shows, $n=0.63$ exhibits more shear-thinning characteristics than $n=0.6$. In
addition, due to the kinematic viscosity distribution in the domain, the
Reynolds number in the case 5 is approximately two times larger than that in
the case 4.

Figure 14: Final (a) concentration, (b) kinematic viscosity, and (b) velocity
contours of case 4.  Figure 15: Final (a) concentration, (b) kinematic
viscosity, and (b) velocity contours of case 5.

In some cases, for instance case 6, as shown in Fig. 16 (b), two regions at
the top and bottom have constant, high kinematic viscosity because these
regions experience minimum and equal shear rate at specific times. When the
bubble rises, it exerts large shear stress around itself, resulting in a
decrease in the viscosity of surrounding shear-thinning fluid. Therefore, the
bottom region retrieve its high, constant kinematic viscosity when the bubble
take away from it, and this constant-viscosity region get larger as long as
the bubble goes upward. As a result, when the bubble gets close to the top
boundary, shear rate becomes variable and larger at the top region,
transforming the top region to a small, high-viscosity region.

Figure 16: Final (a) concentration, (b) kinematic viscosity, and (b) velocity
contours of case 6.

The rising bubble pushes the surrounding fluid upward, forcing it to
circulate. As a result, two wakes appear, experiencing low velocity as shown
in Fig. 11 (c)- 16 (c), which make the bubble rear flattened or skirted,
especially in high Eötvös numbers.

##  5 Conclusions

In this paper, a power-law scheme was implemented in a LB model to capture
both shear-thinning and Newtonian behavior of a fluid. A variable mobility was
introduced in this study to reduce interface dissipation between phases. A
comparison between constant and variable mobility in a three-component system
was presented, indicating that the variable mobility can preserve mass of
components better.

Two cases were investigated to test the accuracy of the model. A single phase
power-law fluid flow driven by a constant inlet velocity between two parallel
plates was simulated to show the small relative error between the numerical
and analytical fully-developed velocity profile. Also, a two-phase power-law
fluid flow driven by pressure between two parallel plates was investigated at
different power-law indices and viscosity ratios to show the small relative
errors between numerical simulation and analytical solution of velocity
profiles. The results showed the strength of the present model in simulating
power-law and Newtonian fluids with wide range of power-law indices.

The dynamics of a rising bubble in the shear-thinning liquid was examined
numerically using a LB model. Results show that the rising bubble deforms at
high Eötvös number and reaches a steady state shape and velocity. Also, the
bubble moves upward faster and has higher location at high Reynolds number.
The terminal velocity obtained by an experimental drag coefficient was
compared to the simulated terminal velocity, demonstrating that the results of
the present study are in good agreement with the experiment. The average,
maximum, and minimum relative errors for six cases were $5.6\%$, $19\%$, and
$0$, respectively. When the bubble rises in a shear-thinning liquid, the
viscosity around the bubble decreases because of the increase of shear rate.
Therefore, the viscosity nearby the bubble, which experiences the highest
shear rate, is the lowest and it increases gradually when gets farther away
from the bubble. Also, the regions which are too far from the bubble have
constant and maximum viscosity. The velocity contour and streamlines show that
the rising bubble generates two wakes around itself, containing low or
negative velocity.

In future study, we will use multiple relaxation time (MRT) to capture higher
viscosity ratios. Also, higher Reynolds numbers should be captured to
investigate wider range of dimensionless numbers.

## Conflict of interest

The authors declare no conflict of interest.

## Data Availability

The data used to generate the results in the manuscript are available from the
corresponding author upon reasonable request.

## References

  * [1] R. P. Chhabra, Bubbles, drops, and particles in non-Newtonian fluids, CRC press, 2006.
  * [2] R. P. Chhabra, J. F. Richardson, Non-Newtonian flow and applied rheology: engineering applications, Butterworth-Heinemann, 2011.
  * [3] G. Pontrelli, E. J. Carr, A. Tiribocchi, S. Succi, Modeling drug delivery from multiple emulsions, Phys. Rev. E 102 (2020) 023114.  doi:10.1103/PhysRevE.102.023114.
  * [4] J. Farhadi, A. Sattari, P. Hanafizadeh, Passage of a rising bubble through a liquid-liquid interface: A flow map for different regimes, The Canadian Journal of Chemical Engineering 100 (2) (2022) 375–390.  doi:10.1002/cjce.24078.
  * [5] Y. Wang, Z. Chen, F. Bian, L. Shang, K. Zhu, Y. Zhao, Advances of droplet-based microfluidics in drug discovery, Expert Opinion on Drug Discovery 15 (8) (2020) 969–979.  doi:10.1080/17460441.2020.1758663.
  * [6] On the terminal velocity of single bubbles rising in non-newtonian power-law liquids, Journal of Non-Newtonian Fluid Mechanics 278 (2020) 104249.  doi:j.jnnfm.2020.104249.
  * [7] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, Y.-J. Jan, A front-tracking method for the computations of multiphase flow, Journal of Computational Physics 169 (2) (2001) 708–759.  doi:10.1006/jcph.2001.6726.
  * [8] C. Hirt, B. Nichols, Volume of fluid (vof) method for the dynamics of free boundaries, Journal of Computational Physics 39 (1) (1981) 201–225.  doi:10.1016/0021-9991(81)90145-5.
  * [9] S. Pillapakkam, P. Singh, A level-set method for computing solutions to viscoelastic two-phase flow, Journal of Computational Physics 174 (2) (2001) 552–578.  doi:10.1006/jcph.2001.6927.
  * [10] W. Dijkhuizen, M. van Sint Annaland, J. Kuipers, I. Roghair, Dns of gas bubbles behaviour using an improved 3d front tracking model—drag force on isolated bubbles and comparison with experiments (2010).  doi:10.1016/j.ces.2009.10.021.
  * [11] S. Li, Y. Ma, S. Jiang, T. Fu, C. Zhu, H. Z. Li, The Drag Coefficient and the Shape for a Single Bubble Rising in Non-Newtonian Fluids, Journal of Fluids Engineering 134 (8) (07 2012).  doi:10.1115/1.4007073.
  * [12] M. T. Islam, P. B. Ganesan, J. Cheng, M. S. Uddin, Single bubble rising behaviors in newtonian and non-newtonian fluids with validation of empirical correlations: A computational fluid dynamics study, Engineering Reports 2 (1) (2020) e12100.  doi:10.1002/eng2.12100.
  * [13] M. Ohta, Y. Yoshida, M. Sussman, A computational study of the dynamic motion of a bubble rising in carreau model fluids, Fluid Dynamics Research 42 (2) (2009) 025501.  doi:10.1088/0169-5983/42/2/025501.
  * [14] A. Premlata, M. K. Tripathi, B. Karri, K. C. Sahu, Dynamics of an air bubble rising in a non-newtonian liquid in the axisymmetric regime, Journal of Non-Newtonian Fluid Mechanics 239 (2017) 53–61.  doi:10.1016/j.jnnfm.2016.12.003.
  * [15] M. Pang, M. Lu, Numerical study on dynamics of single bubble rising in shear-thinning power-law fluid in different gravity environment, Vacuum 153 (2018) 101–111.  doi:10.1016/j.vacuum.2018.04.011.
  * [16] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. M. Viggen, The lattice boltzmann method.
  * [17] M. Majidi, M. A. Bijarchi, A. G. Arani, M. H. Rahimian, M. B. Shafii, Magnetic field-induced control of a compound ferrofluid droplet deformation and breakup in shear flow using a hybrid lattice boltzmann-finite difference method, International Journal of Multiphase Flow 146 (2022) 103846.  doi:10.1016/j.ijmultiphaseflow.2021.103846.
  * [18] J. Farhadi, V. Bazargan, Marangoni flow and surfactant transport in evaporating sessile droplets: A lattice boltzmann study, Physics of Fluids 34 (3) (2022) 032115\.
  * [19] A. K. Gunstensen, D. H. Rothman, S. Zaleski, G. Zanetti, Lattice boltzmann model of immiscible fluids, Phys. Rev. A 43 (1991) 4320–4327.  doi:10.1103/PhysRevA.43.4320.
  * [20] X. Shan, H. Chen, Lattice boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815–1819.  doi:10.1103/PhysRevE.47.1815.
  * [21] M. R. Swift, E. Orlandini, W. R. Osborn, J. M. Yeomans, Lattice boltzmann simulations of liquid-gas and binary fluid systems, Phys. Rev. E 54 (1996) 5041–5052.  doi:10.1103/PhysRevE.54.5041.
  * [22] A lattice boltzmann scheme for incompressible multiphase flow and its application in simulation of rayleigh–taylor instability, Journal of Computational Physics 152 (2) (1999) 642–663.  doi:https://doi.org/10.1006/jcph.1999.6257.
  * [23] T. Inamuro, T. Ogata, S. Tajima, N. Konishi, A lattice boltzmann method for incompressible two-phase flows with large density differences, Journal of Computational Physics 198 (2) (2004) 628–644.  doi:10.1016/j.jcp.2004.01.019.
  * [24] T. Lee, Effects of incompressibility on the elimination of parasitic currents in the lattice boltzmann equation method for binary fluids, Computers & Mathematics with Applications 58 (5) (2009) 987–994, mesoscopic Methods in Engineering and Science.  doi:10.1016/j.camwa.2009.02.017.
  * [25] R. Haghani Hassan Abadi, A. Fakhari, M. H. Rahimian, Numerical simulation of three-component multiphase flows at high density and viscosity ratios using lattice boltzmann methods, Phys. Rev. E 97 (2018) 033312.  doi:10.1103/PhysRevE.97.033312.
  * [26] Boyer, Franck, Lapuerta, Céline, Study of a three component cahn-hilliard flow model, ESAIM: M2AN 40 (4) (2006) 653–687.  doi:10.1051/m2an:2006028.
URL https://doi.org/10.1051/m2an:2006028

  * [27] A. Ghorbanpour-Arani, M.-H. Rahimian, R. Haghani-Hassan-Abadi, Numerical simulation of dissolved air flotation using a lattice boltzmann method, Phys. Rev. E 101 (2020) 023105.  doi:10.1103/PhysRevE.101.023105.
  * [28] S. Gabbanelli, G. Drazer, J. Koplik, Lattice boltzmann method for non-newtonian (power-law) fluids, Phys. Rev. E 72 (2005) 046312.  doi:10.1103/PhysRevE.72.046312.
  * [29] J. Boyd, J. Buick, S. Green, A second-order accurate lattice boltzmann non-newtonian flow model, Journal of Physics A: Mathematical and General 39 (46) (2006) 14241–14247.  doi:10.1088/0305-4470/39/46/001.
URL https://doi.org/10.1088/0305-4470/39/46/001

  * [30] Y. Ba, N. Wang, H. Liu, Q. Li, G. He, Regularized lattice boltzmann model for immiscible two-phase flows with power-law rheology, Phys. Rev. E 97 (2018) 033307.  doi:10.1103/PhysRevE.97.033307.
URL https://link.aps.org/doi/10.1103/PhysRevE.97.033307

  * [31] N. Wang, H. Liu, C. Zhang, Deformation and breakup of a confined droplet in shear flows with power-law rheology, Journal of Rheology 61 (4) (2017) 741–758.  doi:10.1122/1.4984757.
  * [32] C. Xie, J. Zhang, V. Bertola, M. Wang, Lattice boltzmann modeling for multiphase viscoplastic fluid flow, Journal of Non-Newtonian Fluid Mechanics 234 (2016) 118–128.  doi:10.1016/j.jnnfm.2016.05.003.
  * [33] H. Liu, A. J. Valocchi, Y. Zhang, Q. Kang, Phase-field-based lattice boltzmann finite-difference model for simulating thermocapillary flows, Phys. Rev. E 87 (2013) 013010.  doi:10.1103/PhysRevE.87.013010.
  * [34] Cahn-hilliard vs singular cahn-hilliard equations in phase field modeling, Communications in Computational Physics 7 (2) (2010) 362–382.  doi:0.4208/cicp.2009.09.016.
  * [35] T. N. Phillips, G. W. Roberts, Lattice Boltzmann models for non-Newtonian flows, IMA Journal of Applied Mathematics 76 (5) (2011) 790–816.  doi:10.1093/imamat/hxr003.
  * [36] S. A. Hosseini, F. Huang, D. Thévenin, Lattice boltzmann model for simulation of flow in intra-cranial aneurysms considering non-newtonian effects, Physics of Fluids (2022).
  * [37] F. Irgens, Rheology and non-newtonian fluids, Vol. 1, Springer, 2014.
  * [38] O. MALASPINAS, G. COURBEBAISSE, M. DEVILLE, Simulation of generalized newtonian fluids with the lattice boltzmann method, International Journal of Modern Physics C 18 (12) (2007) 1939–1949.  doi:10.1142/S0129183107011832.
  * [39] J. Kim, A numerical method for the cahn–hilliard equation with a variable mobility, Communications in Nonlinear Science and Numerical Simulation 12 (8) (2007) 1560–1571.
  * [40] Y. Q. Zu, S. He, Phase-field-based lattice boltzmann model for incompressible binary fluid systems with density and viscosity contrasts, Phys. Rev. E 87 (2013) 043301.  doi:10.1103/PhysRevE.87.043301.
URL https://link.aps.org/doi/10.1103/PhysRevE.87.043301

  * [41] Q. Lou, Z. Guo, B. Shi, Evaluation of outflow boundary conditions for two-phase lattice boltzmann equation, Phys. Rev. E 87 (2013) 063301.  doi:10.1103/PhysRevE.87.063301.
URL https://link.aps.org/doi/10.1103/PhysRevE.87.063301

  * [42] R. Chhabra, Motion of spheres in power law (viscoinelastic) fluids at intermediate reynolds numbers: a unified approach, Chemical Engineering and Processing: Process Intensification 28 (2) (1990) 89–94.  doi:10.1016/0255-2701(90)80004-O.
  * [43] N. Ioannou, H. Liu, M. S. Oliveira, Y. Zhang, Droplet dynamics of newtonian and inelastic non-newtonian fluids in confinement, Micromachines 8 (2) (2017) 57.


XXXX-XXXX

# Direction-sensitive dark matter search with a low-background gaseous
detector NEWAGE-0.3b”

Tomonori Ikeda<EMAIL_ADDRESS>Department of Physics,
Graduate School of Science, Kobe University, Rokkodai-cho, Nada-ku, Kobe-shi,
Hyogo, 657-8501, Japan  Kiseki Nakamura  Kamioka Observatory, Institute for
Cosmic Ray Research, the University of Tokyo, Higashi-Mozumi, Kamioka-cho,
Hida-shi, Gifu, 506-1205, Japan  Takuya Shimada  Ryota Yakabe  Takashi
Hashimoto  Hirohisa Ishiura  Takuma Nakamura  Hiroshi Ito  Koichi Ichimura
Research Center for Neutrino Science, Tohoku University, Sendai 980-8578,
Japan  Ko Abe  Kavli Institute for the Physics and Mathematics of the Universe
(WPI), the University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa-shi, Chiba,
277-8582, Japan  Kazuyoshi Kobayashi  Waseda Research Institute for Science
and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo 169-8555,
Japan  Toru Tanimori  Division of Physics and Astronomy, Graduate School of
Science, Kyoto University, Kitashirakawaoiwake-cho, Sakyo-ku, Kyoto-shi,
Kyoto, 606-8502, Japan  Hidetoshi Kubo  Atsushi Takada  Hiroyuki Sekiya
Atsushi Takeda  Kentaro Miuchi

###### Abstract

NEWAGE is a direction-sensitive dark matter search using a low-pressure
gaseous time projection chamber. A low alpha-ray emission rate micro pixel
chamber had been developed in order to reduce background for dark matter
search. We conducted the dark matter search at the Kamioka Observatory in
2018. The total live time was 107.6 days corresponding to an exposure of 1.1
kg$\cdot$days. Two events remained in the energy region of 50-60 keV which was
consistent with 2.5 events of the expected background. A directional analysis
was carried out and no significant forward-backward asymmetry derived from the
WIMP-nucleus elastic scatterings was found. Thus a 90% confidence level upper
limit on Spin-Dependent WIMP-proton cross section of 50 pb for a WIMP mass of
100 GeV/$c^{2}$ was derived. This limit is the most stringent yet obtained
from direction-sensitive dark matter search experiments.

Dark matter, MPGD, $\mu$TPC

##  1 Introduction

Dark matter is one of the biggest puzzles of the modern cosmology and the
particle physics. A number of experimental efforts aiming to find the Weakly
Interacting Massive Particle (WIMP) dark matter through direct searches which
observe the scatterings of the WIMP and nuclei have been carried out XENON ;
LUX ; CDMS ; DarkSide ; DEAP3600 ; EDELWEISS ; PANDAX ; XMASS . However, the
dark matter has not been discovered yet. In the direct search, the annual
modulation and the directional signature would be two possible signals among
the characteristic signals of the dark matter. The annual modulation is caused
by the orbital motion of the Earth around the Sun. The modulation amplitude is
expected to be a few % BAUM2019262 . On the other hand, the directional
signature is due to the circular motion of the solar system around the galaxy
center. The forward-backward ratio in nuclear recoil angular distribution
derived from the WIMP-nucleus elastic scatterings could be an order of
magnitude Spergel . In addition, the directional method could discover the
WIMP dark matter beyond the neutrino floor, which represents the ultimate
background by the coherent neutrino-nucleus scatterings PhysRevD.92.063518 ,
and reveal the astrophysical and particle properties of the dark matter
PhysRevD.83.075002 ; Lee_2012 ; NAGAO2020100426 .

NEWAGE (NEw generation WIMP search with an Advanced Gaseous tracker
Experiment) is a direction-sensitive WIMP dark matter search experiment using
a low pressure gas micro time projection chamber ($\mu$TPC) Miuchi_2003 . A
direction-sensitive experiment needs to detect the direction of the recoil
nuclei. Hence a TPC at low pressure gas and a readout device, $\mu$-PIC Takada
, which is one variation of micro pattern gaseous detectors, are used. In
2015, an underground measurement was performed and the best directional
constraint testing the forward-backward asymmetry in the nuclear recoil
angular distribution was achieved PTEPNEWAGE . Then we increased the
statistics by a factor of more than ten and the first 3d-vector directional
dark matter search was performed Yakabe2020 . However, a certain amount of
radioactivity, which potentially contributed to the background, was later
found inside the $\mu$-PIC. Therefore the surface material of the $\mu$-PIC,
which was the dominant background source, was replaced with less radioactive
material. This newly developed low-background $\mu$-PIC was called a low
alpha-ray emission rate micro pixel chamber (LA$\mu$-PIC) LAu-PIC . In this
paper, the first results of the direction-sensitive dark matter search using
the LA$\mu$-PIC are reported.

##  2 NEWAGE-0.3b” detector

The NEWAGE-0.3b’ detector was upgraded to NEWAGE-0.3b” by replacing the
readout device from the standard $\mu$-PIC to LA$\mu$-PIC. Most part of the
detector system is unchanged and we briefly summarize the structure of the
NEWAGE-0.3b” detector and its performance in this section.

###  2.1 System

The NEWAGE-0.3b” detector is a low-pressure gaseous $\mu$TPC that is comprised
of a LA$\mu$-PIC LAu-PIC , a gas electron multiplier (GEM Sauli ) and a TPC
cage. The schematic drawing is shown in Fig. 1. The LA$\mu$-PIC was
manufactured by Dai Nippon Printing Co. Lt.. It has 768 $\times$ 768 pixels
with a pitch of 400 $\mu$m forming the detection area of 30.7 $\times$ 30.7
cm2. These electrodes are connected by 768 anode strips and 768 cathode
strips. The anode and cathode strips of LA$\mu$-PIC are orthogonally formed
and thus a two-dimensional position of a hit pixel can be known. The structure
of the pixel electrode is same as that of a standard $\mu$-PIC Takada , while
the material facing the detection volume was changed. Measurements of the U/Th
contamination in the $\mu$-PIC components using a high purity germanium
detector indicated that the glass cloth-sheet has a large amount of
radioactive contamination. Therefore, a compound of epoxy and polyimide
without the glass cloth-sheet, which is a factor of one hundred less
contaminated by isotopes of 238U and 232Th, was chosen as a new surface
material. Performance test showed that the newly developed $\mu$–PIC works as
the standard $\mu$-PIC does. Details can be found in Ref. HASHIMOTO2020 . A
GEM with an effective area of 32 $\times$ 31 cm2 is used as a first stage
amplifier in order to obtain a sufficient gas gain while keeping a stable
operation. The substrate of the GEM is a 100 $\mu$m thick liquid crystal
polymer. Cylindrical holes are formed using the laser etching technique. A
hole size and pitch are 70 $\mu$m and 140 $\mu$m, respectively TAMAGAWA . The
TPC field cage, which is made of four plates of polyether ether ketone plastic
plates, was installed in order to make a uniform electric field. Copper wires
with a spacing of 1 cm are placed on the side walls and chained by resistors.
The TPC cage has a length of 41 cm. The vessel was filled with CF4 at 76 Torr.
CF4 is chosen because of its small diffusion and a large cross section for the
Spin-Dependent (SD) interaction of fluorine. A gas circulation system with a
cooled charcoal (TSURUMICOAL 2GS) of 100 g was installed in order to remove
radons.

Figure 1: Schematic drawings of the NEWAGE-0.3b” detector. The left-top
microscope picture shows the pixel electrode structure of the LA$\mu$-PIC.
Left-bottom drawing is the cross section of the LA$\mu$-PIC. Surface material
is combination of epoxy and polyimide without the glass cloth-sheet (GC).
NEWAGE-0.3b” chamber is shown in the right panel. LA$\mu$-PIC has a detection
area of 31 $\times$ 31 cm2 and the drift length is 41 cm. The gap region is
defined as a volume between the LA$\mu$-PIC and the GEM. The 10B glass plate
is set inside the TPC field cage for the energy calibration.

The charge signals are read out from the anode and cathode electrodes with
Amplifier-Shaper-Discriminator chips (SONY CXA3653Q ASD_orito ). These signals
are divided into digitized and analog signals. The digitized signals are sent
to the FPGA-based encoding system and the hit-patterns are recorded with a
clock of 100 MHz. In addition, the time-over-thresholds (TOTs) of each strip
are recorded. The track-length and direction are reconstructed using these
information (see Ref. Yakabe2020  for details). Analog signals of the 768
cathode strips are grouped into four channels and each channel is then divided
into two. One of the divided signals is directly connected to a waveform
digitizer (REPIC RPV160, 100 MHz), while the other is attenuated by a factor
of three and read by the waveform digitizer. The waveforms are mainly used to
determine the energy deposit of each event.

The energy calibration is performed with alpha-rays generated in a
${}^{10}{\rm B}(n,\alpha)^{7}{\rm Li}$ reaction. A glass plate with a thin 10B
layer is installed inside the $\mu$TPC. The detector is irradiated with
neutrons from a 252Cf fission source placed outside of the vessel and the
neutrons are thermalized by polyethylene blocks. A continuous spectrum with a
maximum edge at 1.5 MeV because of the thickness of the 10B layer is obtained.
The measured spectrum is fitted with simulated ones and the calibration factor
converting the charge to energy is determined.

###  2.2 Performance

The following four analysis cuts were applied in order to select nuclear
recoil events. Here the 252Cf neutron calibration data was taken without
polyethylene blocks in order to induce the fast neutron-nuclear scatterings.

  * •

Fiducial cut: The fiducial area of XY plane was defined as 28 $\times$ 24 cm2.
The whole part of the track is required to be in the fiducial volume. This cut
removes the charged particles from the wall and the 10B glass plate.

  * •

Length-Energy cut: The stopping powers of nuclei are larger than those of
electrons. Hence the track-length vs. energy distribution can be used to
identify the nuclear recoil tracks. Figure 2 (a) shows the track-length vs.
energy distributions for the 252Cf calibration (black points) and the 137Cs
calibration (blue points). From the 252Cf calibration data, we determined the
nuclear band by fitting with Gaussian function for every energy bin with a
width of 10 keV. This cut removes electrons and alpha-rays.

  * •

TOT-Energy cut: Energy deposition on a single strip is stored as TOT. Since
nuclear recoil events have larger energy losses than electron events, the sum
of TOT (TOT-sum) for given energy tends to be large. In addition, since TOT-
sum was expected to be linear with respect to the total energy deposit, a
parameter defined by dividing the TOT-sum by energy was suitable for this
purpose. Figure 3 (a) shows TOT-sum/energy vs. energy distributions for the
252Cf calibration (black points) and the 137Cs calibration (blue points). From
the 252Cf calibration data, we determined the nuclear band by fitting with
gaussian function every energy bin with a width of 10 keV. This cut rejects
electron events.

  * •

Roundness cut: The roundness is defined as a reduced chi-square value of a
linear fit to the track. While the self-triggering mode TPC is not able to
measure the absolute z position, this parameter is known to show some
correlation with the absolute z position. For instance, small z events have
short drift lengths and thus, the effect of electron diffusion is negligible.
The track information keeps its original shape and fitted with a straight line
well. On the other hand, large z events are largely affected by electron
diffusions, and the fitness to a straight line would be worse. In other words,
these large z events have large roundness. Hence the roundness correlates with
the absolute z position. Especially, this parameter is useful to remove “gap
events”. The gap events are defined as events depositing their all energy in
the gap region (Fig. 1). These gap events are not amplified by the GEM and
thus the measured charge is smaller by a factor of the GEM gain than those in
the detection volume. In order to reproduce the gap events, we irradiated the
detector with neutrons from a 252Cf fission source without a drift electric
field. Figure 4 (a) shows the roundness vs. energy distribution for the gap
events. The roundness of gap events (red points) was found to be small. The
events with roundness $>$ 0.05 are selected in order to remove gap events.

Figure 2: Track-length vs. energy distributions. (a) The gradation and red
points represent the 252Cf neutron calibration data and the 137Cs electron
calibration data, respectively. The red solid and dotted lines indicate the
median and $\pm$3$\sigma$ quantiles of the neutron calibration. (b) Scientific
RUN data (DM RUN) after the fiducial cut. Figure 3: TOT-sum/energy vs. energy
distributions. (a) The gradation and red points represent the 252Cf neutron
calibration data and the 137Cs electron calibration data, respectively. The
red solid and dotted lines indicate the median and $\pm$3$\sigma$ quantiles of
the neutron calibration. (b) Scientific RUN data (DM RUN) after the fiducial
cut. Figure 4: Roundness vs. energy distribution. (a) The gradation and red
points represent the 252Cf neutron calibration with a drift electric field and
without a drift electric field, respectively. The cyan dotted line
(roundness$=0.05$) indicate the nuclear selection line. (b) Scientific RUN
data (DM RUN) after the fiducial cut.

The nuclear detection efficiency was evaluated by irradiating the $\mu$TPC
with neutrons from a 252Cf fission source. In order to cancel the x, y and z
position dependence of $\mu$TPC and evaluate overall efficiency, we needed
homogeneous irradiation throughout the detection volume. We confirmed by Monte
Carlo simulations that homogeneous irradiation in the TPC is made when 252Cf
fission sources were set six positions. In the actual measurement, we set a
252Cf source at one of the six positions and took the data, then move the
source to another position. We performed the measurements six times in series
and combined the data in the analysis to realize homogeneous irradiation. The
detection efficiency was calculated by comparing the measured energy and the
simulated one. The measured detection efficiencies of the nuclear recoil
events are shown in Fig. 5. The black and red lines are fitted lines of
experimental data and indicate the detection efficiency when we applied only
the fiducial cut and all event selections. The efficiency with all event
selections was 14% at 50 keV.

The $\mu$TPC has a non-isotropic response with regard to the nuclear recoil
track direction because of the track-reconstruction algorithm. Thus, we need
to measure the relative direction-dependent efficiency of nuclear recoil in
the energy range of 50-100 keV. Figure 6 shows the measured distribution of
the elevation angle $\theta_{\rm ele}$ and the azimuth angle $\phi_{\rm azi}$
in the detector coordinate. Although an ideal directional detector is expected
to have a uniform distribution, our detector is not optimized yet in terms of
the uniformity for the directional response. The efficiency is low around XY
plane, XZ plane and YZ plane in any energy region because of the poor track
reconstruction along anode/cathode strips. In addition, there are higher
efficiency areas along the diagonal lines from the direction of the anode and
the cathode strips. This is because the current tracking algorithm tends to
recognize the diffused tracks to such directions. We evaluated the gamma
rejection power, or the detection efficiency of electrons, by irradiating the
detector with gamma-rays from a 137Cs source. The gamma rejection power or the
electron detection efficiency for the energy range of 50-60 keV was
$1.3^{+3.0}_{-1.1}\times 10^{-6}$.

Figure 5: Typical detection efficiency for nuclear recoil events. The gray and
black points with errors are experimental data after the fiducial cut and the
roundness cut, respectively. The fitted black solid line is the detection
efficiency in the fiducial volume. The fitted red solid line is the detection
efficiency after all event selections. Figure 6: Relative direction-dependent
efficiencies in each energy region. The axis with white and black label are
the azimuth angle $\phi_{\rm{azi}}$ and the elevation angle
$\theta_{\rm{ele}}$ in the detector coordinate, respectively. The x-axis and
y-axis correspond to $(\phi_{\rm azi},\theta_{\rm ele})=(0,0)$ and $(\phi_{\rm
azi},\theta_{\rm ele})=(90,0)$, respectively.

The angular resolution was evaluated using the neutron-nuclei elastic
scatterings of same method in previous study Nakamura_2012 . The angular
$\theta$ is defined as the angle between the direction of the scattered nuclei
and the neutron source. We evaluated the angular resolution by the comparison
of measured and simulated distributions of the recoil angle $\cos\theta$. The
obtained angular resolution was 48.0${}^{+6.8}_{-2.2}$ degree in the energy
range of 50-100 keV. The expected forward-backward ratio of
$|\cos\theta_{\rm{cygnus}}|$ for the 100 GeV/c2 WIMP is 30%. Here
$\theta_{\rm{cygnus}}$ is defined as the angle between the WIMP-wind direction
and the measured direction of the recoil nucleus.

##  3 Experiment

A directional dark matter search (RUN22) was carried in Laboratory B of the
Kamioka Observatory (36.25’N, 137.18’E) located at 2700 m water equivalent
underground. The LA$\mu$-PIC plane was placed vertically and the $z$-axis is
aligned to the direction of S30∘E. The first sub RUN was performed from Jun.
6th 2018 to 24th Aug. 2018 and the second sub RUN was carried out from 20th
Sep. 2018 to 14th Nov. 2018. The target gas is CF4 at 76 Torr (0.1 bar) and
the target mass in the fiducial volume of 28 $\times$ 24 $\times$ 41 cm3 (28
L) is 10 g. The total live time is 107.6 days corresponding to an exposure of
1.1 kg$\cdot$days.

The energy calibration and the detection efficiency measurement were carried
out every two weeks. The gas gain at the beginnings of the sub RUNs was about
1100 and time-dependent variation was observed due to the gas deterioration.
The energy scale of the data was corrected considering the time-dependence of
the gas gain. The energy resolution was $13.2\pm 2.3\%$ above 50 keV. The
measured drift velocity at the beginnings of the sub RUNs was 9.6 cm/$\mu$s
and the gas deterioration gave $+3.3\%$ and $-12\%$ of the maximum
uncertainty.

Length-Energy, the TOT-Energy and the Roundness-Energy distributions after the
fiducial cut for the dark matter search data are shown in Fig. 2 (b), 3 (b)
and 4 (b), respectively. A large fraction of the events has long track lengths
and small TOTs, which indicates that most of the measured events are
electrons. These events are effectively reduced by analysis cuts introduced in
Sec. 2.2. Energy spectra at each cut stage are shown in Fig. 7 (a). The final
event sample was reduced to 17 events in the energy region of 50-400 keV owing
to analysis cuts using track information without any shields like Pb. Since
the directionality is lost in the energy below 50 keV because of the short
tracks and the diffusion effect, the lower energy bound is set at 50 keV. The
energy spectrum for the final sample unfolded by the nuclear detection
efficiency is shown in Fig. 7 (b) together with one of the previous results
RUN14 PTEPNEWAGE  using the standard $\mu$-PIC. The main background of RUN14
in 50-100 keV are alpha-rays radiated from the surface material of the
standard $\mu$-PIC. These backgrounds were reduced in RUN22 by about factor 10
thanks to the LA$\mu$-PIC whose surface material is less contaminated with
238U and 232Th. It is demonstrated that the LA$\mu$-PIC works as expected to
reduce the alpha-ray backgrounds.

Figure 8 shows the skymap of the final sample on the detector coordinate. We
calculated the nuclear recoil distribution $|\cos\theta_{\rm{cygnus}}|$ for
the energy region of 50-100 keV in order to evaluate the forward-backward
asymmetry. Figure. 9 shows the measured and expected
$|\cos\theta_{\rm{cygnus}}|$ distribution binned into two for each energy
range. As a model-independent interpretation of these results, we first set
upper limits on the forward-backward asymmetry parameters based on the raw
counting numbers in each energy bin. The forward-backward asymmetry parameter
is defined as the ratio of the first and second bins. The 90% confidence level
upper limits of these values in 50-60 keV, 60-70 keV, 70-80 keV and 80-90 keV
are 2.3, 3.9, 5.3 and 3.9, respectively. The systematic uncertainties of the
expected rate for the WIMP were summarized in Table 1. The angular resolution
gives the dominant systematic uncertainty, which impacts the shape of the
nuclear recoil distribution, and is considered in the following statistic
test.

Figure 7: (a) Energy spectra of the dark matter search at each selection step. The black, red, blue and green lines are no-cut energy spectrum, the one after the Fiducial cut, Length-Energy cut and TOT-Energy cut, respectively. The black points with error bars are the final event sample after the Roundness cut. (b) Final energy spectrum considering the detection efficiency. The black and gray points with error bars represent RUN22 (using the LA$\mu$-PIC) and RUN14 (using the standard $\mu$-PIC), respectively. Error bars indicate statistic poisson errors. The blue dotted line shows the expected spectrum of the WIMP-nucleus scattering with the WIMP mass of 100 GeV/c2, the WIMP-proton cross section of $\sigma_{\chi}=20~{}{\rm pb}$ and the energy resolution of 13.2%. Figure 8: Skymap of the final sample on the detector coordinate. The black, red and blue points show the energy region of 50-100 keV, 100-200 keV and 200-400 keV, respectively. Figure 9: Measured and expected $|\rm{cos}\theta_{\rm{cygnus}}|$ distribution for each energy range. The black points with errors are measured data. The cyan lines show the WIMP expected signal of the WIMP-nucleus scattering with the WIMP mass of 100 GeV/c2 and the WIMP-proton cross section of $\sigma_{\chi}=20~{}{\rm pb}$. Table 1: Systematic uncertainties of the expected rate for the WIMP mass $m_{\chi}=100$ GeV/c2. Source |  $\cos\theta$ range | Relative uncertainty (%)
---|---|---
Energy resolution | [ 0, 1 ] | $<0.1$
Drift velocity | [ 0, 0.5 ] | $<0.2$
[ 0.5, 1 ] | $<0.2$
Angular resolution | [ 0, 0.5 ] | +5.4 $-2.2$
[ 0.5, 1 ] | +1.7 $-4.2$

##  4 Results

In order to obtain a possible anisotropic $|\rm{cos}\theta_{\rm{cygnus}}|$
distribution, a binned likelihood-ratio method was used PhysRevD.67.012002 .
The minimized statistic value $\chi^{2}$ was defined as,

$\chi^{2}=2\sum_{i=0}^{n}\biggl{[}(N_{i}^{\rm{exp}}-N_{i}^{\rm{data}})+N_{i}^{\rm{data}}{\rm{ln}}\biggl{(}\frac{N_{i}^{\rm{data}}}{N_{i}^{\rm{exp}}}\biggr{)}\biggr{]}+\alpha^{2},$
(1)

where the subscript $i$ is the bin number of $|\cos\theta_{\rm{cygnus}}|$
distribution, $N_{i}^{\rm{data}}$ is the measured number of events and
$N_{i}^{\rm{exp}}$ is the expected number of events. A nuisance parameter
$\alpha$ $(=\xi/\sigma_{\kappa})$ was introduced to consider the systematic
uncertainty of the angular resolution $\sigma_{\kappa}$. A possible angular-
resolution shift is $\xi$. The expected event number $N_{i}^{\rm{exp}}$ was
obtained by a signal Monte Carlo simulation. It depends on the WIMP mass
$m_{\chi}$, the WIMP-proton cross section $\sigma_{\chi-p}$ and the
astrophysical parameters. In addition, the nuclear quenching factor and the
detector responses were considered. Here the nuclear quenching factor was
simulated by SRIM SRIM , which represented the experimental alpha-ray data of
the previous experiment NISHIMURA2009185 . The energy bin was divided into 10
keV/bin considering the energy resolution. Since each energy bin had low
statistics, the measured $|\cos\theta_{\rm{cygnus}}|$ were binned into 2-bin.

The measurement data was fitted by minimizing $\chi^{2}$ using an anisotropic
WIMP model. Here, no background (BG) was added in order to prevent any
uncertainty of the BG distribution. The WIMP-proton cross section
$\sigma_{\chi-p}$ and the nuisance parameter $\alpha$ were treated as fitting
parameters. The minimum $\chi^{2}$ value for the 50-60 keV bin was 3.3 where
$\sigma_{\chi-p}$ and $\alpha$ were 18.5 pb and 0.12, respectively. Figure 10
shows the measured $|\cos\theta_{\rm{cygnus}}|$ distribution in the energy
region of 50-60 keV along with the expected one using best-fit values. In
order to calculate the p-value, we made a $\chi^{2}$ distribution of an
isotropic BG model and an anisotropic WIMP model from dummy samples. One
thousand dummy samples were produced by Monte Carlo simulations and $\chi^{2}$
value of each dummy sample was calculated. P-values for both of the WIMP and
BG model were 3.3%. Hence we cannot claim the detection of WIMP dark matter
with sufficient significance from the observed data. This is a natural result
because of the large statistic error and the small expected anisotropic ratio.
Since no significant amplitude was found, a 90% confidence level (C.L.) upper
limit was set on the SD cross section. The 90% C.L. upper limit on the SD
cross section was obtained as 50 pb for 100 GeV WIMPs. Figure 11 shows the 90%
C.L. upper limit on the SD WIMP-proton cross section as a function of the WIMP
mass. This result marked a new best sensitivity record of the SD WIMP search
with the direction-sensitive method. This result improved the constraint by
about 15 times compared to the previous result of RUN14-18. This improvement
owes to the surface background reduction of the $\mu$-PIC detector.

Figure 10: The $|\cos\theta_{\rm{cygnus}}|$ distribution in the energy range
of 50-60 keV. The black points show observed data. The solid and dotted red
lines are the simulated distribution using best-fit values and excluded values
of 90% C.L., respectively. Figure 11: 90% C.L. upper limits on the SD WIMP-
proton cross section as a function of the WIMP mass. The red thick solid line
is the result of the directional method in this work. The blue solid lines of
RUN14 and RUN14-18 are our previous results PTEPNEWAGE ; Yakabe2020 . The
light blue and black solid lines show the results from the directional
analysis of DMTPC DMTPC  and the conventional analysis without the directional
sensitivity of DRIFT BATTAT20151 , respectively. Allowed region from
DAMA/LIBRA experiment PhysRevD.70.123513  is shown by the green area.

##  5 Discussion

One of the milestones of the NEWAGE is to search an allowed region from DAMA
experiment PhysRevD.70.123513  with the directional method. More than two
orders magnitude of improvements is required in order to cover the entire DAMA
region. Since the sensitivity is limited by the remaining backgrounds, we
investigated the origin of these backgrounds in RUN22 using Geant4 Geant4
simulation.

The main internal background candidates are alpha-rays coming from the decay
of radons and the LA$\mu$-PIC surface. The alpha-rays cannot be discriminated
from nuclear recoil events around 50 keV in NEWAGE-0.3b” detector by analysis
cuts. Hence screening the material of the TPC and the LA$\mu$-PIC is needed. A
broad peak is observed in the energy spectrum around 6 MeV corresponding to
the 220Rn and 222Rn decay during the dark matter search IkedaDron . Estimated
contamination of 220Rn and 222Rn are 5.7 $\pm$ 0.3 mBq/m3 and (5.3 $\pm$ 2.1)
$\times$ 10-1 mBq/m3, respectively. In Geant4 simulation, we generated alpha-
rays from decay chains of 220Rn and 222Rn inside the TPC according to the
branching ratio and estimated the number of events or rate with the analysis
cut same to RUN22. The expected number of events (rate) due to 220Rn and 222Rn
contamination in the low energy region of 50-60 keV were $(6.1\pm 0.7)\times
10^{-1}$ (5.5 $\times$ 10-2 dru111dru$=$counts/sec/keV/kg) and $(5.3\pm
2.2)\times 10^{-2}$ (4.7 $\times$ 10-3 dru), respectively. The remaining
alpha-ray emission rate of the LA$\mu$-PIC was (2.1 $\pm$ 0.5) $\times$ 10-4
alpha/cm2/hr HASHIMOTO2020  and the dominant component was 5.3 MeV alpha-rays
from 210Po decay. The expected number of events (rate) in 50-60 keV region
from this background was less than $1.2\times 10^{-1}$ (1.1 $\times$ 10-2
dru).

Ambient gamma-rays and neutrons from rocks in the mine are two of the main
components of the external background. Contributions of the cosmic-ray muons
are negligible compared with ambient gamma-rays and neutrons. Measured ambient
gamma-rays flux in the Laboratory B Nishimura  was used for the estimation.
The expected counts (rate) in the energy range of 50-60 keV was $1.5\pm 1.5$
(1.4 $\times$ 10-1 dru). The ambient neutrons flux was measured using a 3He
proportional counter and an energy spectrum of ambient neutrons produced by
$(\alpha,n)$ reactions and spontaneous fission was predicted Mizukoshi . The
expected counts (rate) in the energy range of 50-60 keV is $(3.5\pm 0.9)\times
10^{-1}$ (3.1 $\times$ 10-2 dru). The expected number of background events in
the energy region of 50-60 keV is summarized in Table 2. The result without
roundness-cut is also shown in order to confirm the alpha-ray backgrounds from
the LA$\mu$-PIC surface, directly. The measured number of events are in good
agreement with the predicted one within errors in both with and without
roundness-cut.

The 222Rn backgrounds were reduced by the gas circulation system with cooled
charcoal and their contribution was found to be negligible. Ambient gamma-rays
and ambient neutrons contribute some part of the backgrounds. These
backgrounds can be reduced by the external shields comprising of materials
like lead (Pb) and water (H2O). On the other hand, internal backgrounds of
220Rn and the LA$\mu$-PIC surface would remain with the external shields and
they would be dominant backgrounds. A straightforward way of reducing these
backgrounds is to replace the detector components with radiopure materials. A
$\mu$-PIC with a further background reduction is being developed. Another
approach to reduce the background is to detect the absolute $z$ position of
the events. The majority of the remaining backgrounds are known to locate at
low $z$ (around the LA$\mu$-PIC and the GEM) and high $z$ position (around the
drift plate). A discovery of minority carriers in CS2 \+ O2 gas mixtures by
the DRIFT group opened the potential of an absolute $z$ measurement in self-
triggering TPCs BATTAT20151 . We have recently demonstrated a three-
dimensional tracking with a spatial resolution of 130 $\mu$m using the
$\mu$-PIC in a negative ion gas SF6. Simultaneously, the absolute $z$
coordinate was determined with a location accuracy of 16 mm Ikeda2020 . The
negative ion gas TPC enables us to reduce these backgrounds using the $z$
fiducialization, effectively. Thus the use of the negative ion gas TPC is
another promising approach to reduce the surface alpha-ray backgrounds. It
should be noted that the surface background will be one of the ultimate
background sources even after significant efforts of material selection since
radioactive isotopes can be embedded by the decays of radons in a normal
atmosphere even after the production. Thus it is important to take both
possible ways to reduce the background so as to start investigating the DAMA
region and further searches.

Table 2: Summary of the expected numbers of background events and measured numbers in the energy region of 50-60 keV. Source | w/ roundness | w/o roundness
---|---|---
Ambient gamma-rays | $1.5\pm 1.5$ | $4.6\pm 2.7$
Ambient neutrons | $(3.5\pm 0.9)\times 10^{-1}$ | $(4.8\pm 1.2)\times 10^{-1}$
222Rn | $(5.3\pm 2.2)\times 10^{-2}$ | $(8.6\pm 3.5)\times 10^{-2}$
220Rn | $(6.1\pm 0.7)\times 10^{-1}$ | $1.1\pm 0.1$
LA$\mu$-PIC surface | $<1.2\times 10^{-1}$ | $9.1\pm 2.3$
Total background | $2.5\pm 1.5$ | $15\pm 3.5$
Measurement | $2.0\pm 1.4$ | $12\pm 3.5$

##  6 Conclusion

We developed a low background $\mu$TPC detector, NEWAGE-0.3b”, for the
directional dark matter search with an LA$\mu$-PIC, which is a two-dimensional
tracking gaseous detector made of low radioactive materials. A directional
dark matter search in the Kamioka Observatory was carried out in 2018. Total
exposure was 1.1 kg$\cdot$days and the number of the observed events in the
energy region of 50-60 keV was two which was consistent with 2.5 events of the
expected background. No significant forward-backward asymmetry of a WIMP
signal was found, therefore we derived a 90% confidence level upper limits on
the SD WIMP-proton cross section of 50 pb for 100 GeV/c2 WIMPs. We improved
the constraint of the previous result PTEPNEWAGE  by a factor of 15 and marked
the best direction-sensitive limit.

## Acknowledgment

This work was partly supported by JSPS (Japan Society for the promotion of
Science) KAKENHI (Grant-in-Aids for Scientific Research) (grant nos. 16H02189,
19684005, 23684014 and 26104005, 19H05806), JSPS Bilateral Collaborations
(Joint Research Projects and Seminars) program, ICRR Joint-Usage, Program for
Advancing Strategic International Networks to Accelerate the Circulation of
Talented Researches (R2607), and JSPS Research Fellow ( grant nos. 17J03537).

## References

  * (1) E. Aprile et al., Phys. Rev. Lett., 121, 111302 (2018).
  * (2) D. S. Akerib et al., Phys. Rev. Lett., 118, 021303 (2017).
  * (3) R. Agnese et al., Phys. Rev. Lett., 121, 051301 (2018).
  * (4) P. Agnes et al., Phys. Rev. Lett., 121, 081307 (2018).
  * (5) R. Ajaj et al., Phys. Rev. D, 100, 022004 (2019).
  * (6) E. Armengaud et al., Journal of Cosmology and Astroparticle Physics, 05, 019 (2016).
  * (7) X. Cui et al., Phys. Rev. Lett., 119, 181302 (2017).
  * (8) K. Abe et al., Physics Letters B, 759, 272–276 (2016).
  * (9) S. Baum et al., Physics Letters B, 789, 262–269 (2019).
  * (10) D. Spergel, Phys. Rev. D, 37, 1353–1355 (1988).
  * (11) C. A. J. O’Hare et al., Phys. Rev. D, 92, 063518 (2015).
  * (12) J. Billard et al., Phys. Rev. D, 83, 075002 (2011).
  * (13) S.K. Lee and A.H.G Peter, Journal of Cosmology and Astroparticle Physics, 04, 029–029 (2012).
  * (14) K. I. Nagao et al., Physics of the Dark Universe, 27, 100426 (2020).
  * (15) K. Miuchi et al., IEEE Transactions on Nuclear Science, 50, 825–830 (2003).
  * (16) A Takada et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 573(1), 195–199 (2007).
  * (17) K. Nakamura et al., Progress of Theoretical and Experimental Physics, 043F01 (2015).
  * (18) R. Yakabe et al., Progress of Theoretical and Experimental Physics, ptaa147 (2020).
  * (19) T. Hashimoto et al., AIP Conference Proceedings, 1921(1), 070001 (2018).
  * (20) F. Sauli, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 805, 2–24 (2016).
  * (21) T. Hashimoto et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 977, 164285 (2020).
  * (22) T. Tamagawa et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 608(3), 390–396 (2009).
  * (23) R. Orito et al., IEEE Transactions on Nuclear Science, 51, 1337–1342 (2004).
  * (24) K. Nakamura et al., Journal of Instrumentation, 7(02), C02023 (2012).
  * (25) J. Conrad et al., Phys. Rev. D, 67, 012002 (2003).
  * (26) J. F. Ziegler et al., Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 268(11), 1818–1823 (2010).
  * (27) H. Nishimura et al., Astroparticle Physics, 31(3), 185–191 (2009).
  * (28) S. Ahlen et al., Physics Letters B, 695(1), 124–129 (2011).
  * (29) J.B.R. Battat et al., Physics of the Dark Universe, 9-10, 1–7 (2015).
  * (30) C. Savage et al., Phys. Rev. D, 70, 123513 (2004).
  * (31) S. Agostinelli et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 506(3), 250–303 (2003).
  * (32) T. Ikeda, Doctor Thesis, Kobe University (2020).
  * (33) H. Nishimura, Doctor Thesis, Kyoto University (2008).
  * (34) K. Mizukoshi et al., Progress of Theoretical and Experimental Physics, 12, 123C01 (2018).
  * (35) T. Ikeda et al., Journal of Instrumentation, 15(07), P07015 (2020).