# Critical scaling for dense granular flow between parallel plates near
jamming

Michio Otsuki<EMAIL_ADDRESS>Graduate School of Engineering
Science, Osaka University, 1-3 Machikaneyama, Toyonaka, Osaka, 560-8531, Japan
Kenta Hayashi  Graduate School of Engineering Science, Osaka University, 1-3
Machikaneyama, Toyonaka, Osaka, 560-8531, Japan  Kiwamu Yoshii  Graduate
School of Engineering Science, Osaka University, 1-3 Machikaneyama, Toyonaka,
Osaka, 560-8531, Japan  Department of Physics, Nagoya University, Furo-cho,
Chikusa-ku, Nagoya, Aichi, 464-8602, Japanķ

###### Abstract

We numerically study the flow of dense granular materials between parallel
plates driven by an external force. The granular materials form a jammed
solid-like state when the external force is below a critical force, while they
flow like fluids above the critical force. The transition is characterized by
the mass flux. The critical force depends on the average packing fraction and
the distance between the plates. The scaling laws for the critical force and
the mass flux are predicted theoretically based on a continuum model. They are
numerically verified.

##  I Introduction

The dynamics of dense amorphous particles such as granular materials,
suspensions, foams, and emulsions has been the subject of considerable
investigation [1, 2, 3]. These materials can form a solid-like jammed state
when the packing fraction $\phi$ exceeds a critical fraction $\phi_{\rm c}$,
while they behave like fluids for $\phi<\phi_{\rm c}$. Such a change in a
rheological property is known as the jamming transition. The jammed particles
above $\phi_{\rm c}$ exhibit an elastic response under small shear stress
$\sigma$, but they flow when $\sigma$ exceeds the yield stress $\sigma_{\rm
Y}$ [4]. In the vicinity of $\phi_{\rm c}$, the pressure, shear modulus, yield
stress, and viscosity exhibit continuous transitions, which are characterized
by critical power-law behavior [5, 6, 7, 8, 9, 10, 11, 12]. In particular, it
has been reported that the flow curve [13, 14, 15, 16, 17, 18], complex shear
modulus [19, 20, 21], relaxation modulus [22], mean square displacement [23,
24] obey critical scaling laws, similar to those for equilibrium critical
phenomena.

Most previous studies of the jamming transition have focused on the critical
behavior in homogeneous systems characterized by a constant packing fraction
$\phi$ and shear rate $\dot{\gamma}$. However, in many manufacturing processes
of dense amorphous particles, including granular materials, the formation of a
solid-like jammed state is observed in inhomogeneous systems [25]. This
phenomenon manifests in various forms, including clogging [26, 27], avalanches
[28, 29, 30, 31], and the separation of fluidized surface layers and frozen
bulk regions [32, 33, 34, 35, 36, 37, 38]. Critical behavior is expected for
jamming in these inhomogeneous systems. However, little is known about the
critical behavior in inhomogeneous systems.

To analyze inhomogeneous granular flows, a continuum theory based on the
$\mu(I)$-rheology is useful [25]. In the $\mu(I)$-rheology, a local
constitutive equation is assumed, where the bulk friction

$\displaystyle\mu=|\sigma|/p$ (1)

with the shear stress $\sigma$ and the normal pressure $p$ is described as a
function of the inertia number:

$\displaystyle I=\dot{\gamma}d/\sqrt{p/\rho_{s}},$ (2)

where $\dot{\gamma}$ is the shear rate, $\rho_{\rm s}$ is the solid grain
density, and $d$ is the mean particle diameter [25, 39, 40, 41, 42, 43, 44].
Recent studies have improved the $\mu(I)$-rheology to account for the non-
local effect [45, 44, 46, 47, 48, 49, 50]. The $\mu(I)$-rheology is known to
reproduce the velocity profiles of various granular flows [34, 31, 51, 38,
52]. Therefore, we use it to analyze the critical behavior of the
inhomogeneous flow in the vicinity of jamming.

In this paper, we numerically study the inhomogeneous steady flow of granular
materials driven by an external force between parallel rough plates with the
mean packing fraction $\phi_{0}$. In Sec. II, we explain the setup of our
model. Section III deals with the numerical results of the granular flow. In
Sec. IV, we analytically and numerically investigate scaling laws. The details
of the continuum model are presented in Sec. IV.1. In Sec. IV.2, the scaling
law of the critical force is derived. In Sec. IV.3, the critical behavior of
the mass flux is studied. We discuss and conclude our results in Sec. V. In
Appendix A, we present the numerical results on the pressure and estimate the
critical fraction. The details of the derivation of the scaling laws are
presented in Appendix B.

##  II Setup

Figure 1:  Schematic of granular flows between rough parallel plates. Red line
indicates schematic velocity profiles. Blue particles represent the rough
plates. The size of the fluidized regions are represented by $\ell$.

We consider two-dimensional granular materials consisting of $N$ frictionless
particles with identical density $\rho_{\rm s}$ between two parallel rough
plates with the mean packing fraction $\phi_{0}$ as shown in Fig. 1. The
plates with length $L$ along the $x$-direction are placed at $z=\pm H/2$ with
distance $H$. We impose periodic boundary conditions in the $x$-direction. The
particles are driven by an external force along the $x$-direction.

The equation of motion is given by

$\displaystyle m_{i}\dfrac{d^{2}}{dt^{2}}\bm{r}_{i}=\sum_{j\neq
i}\bm{F}_{ij}+m_{i}f\bm{e}_{x},$ (3)

where $m_{i}$ and $\bm{r}_{i}=(x_{i},z_{i})$ are the mass and position of
particle $i$, respectively. Here, $f$ is the external force per unit mass, and
$\bm{e}_{x}=(1,0)$ is the unit vector along the $x$-direction. The interaction
force is given by

$\displaystyle\bm{F}_{ij}=-\left(ku^{\rm(n)}_{ij}+\eta
v^{\rm(n)}_{ij}\right)\Theta(d_{ij}-r_{ij})\bm{n}_{ij},$ (4)

where $\Theta(x)$ is the Heaviside step function satisfying $\Theta(x)=1$ for
$x\geq 0$ and $\Theta(x)=0$ for otherwise. The elastic constant, viscosity,
and diameter of particle $i$ are denoted by k, $\eta$, and $d_{i}$,
respectively. Here, the relative displacement, the normal relative velocity,
and the normal vector are given by $u^{\rm(n)}_{ij}=r_{ij}-d_{ij}$,
$v^{\rm(n)}_{ij}=\dfrac{d}{dt}u^{\rm(n)}_{ij}$, and
$\bm{n}_{ij}=\bm{r}_{ij}/r_{ij}$, with $d_{ij}=(d_{i}+d_{j})/2$,
$\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$, and $r_{ij}=\left|\bm{r}_{ij}\right|$,
respectively.

The granular materials consist of an equal number of particles with diameters
$d_{0}$ and $sd_{0}$, where $s$ is the size ratio of two types of particles.
The mean diameter is given by $d=(1+s)d_{0}/2$. The mean packing fraction of
this system is given by $\phi_{0}=\pi N(1+s^{2})d_{0}^{2}/(8LH)$. The rough
plates consist of particles with diameter $d_{0}$ placed at $(x,z)=(nd_{0},\pm
H/2)$ with an integer $n$. We first apply an external force large enough to
generate flow, and the external force is changed to a given $f$. The steady
parallel flow is realized after a sufficiently long relaxation time, as shown
in Fig. 1. Particles in a long vertical pipe connected to a hopper [52] and
those between pressure-controlled parallel walls [45, 51] exhibit similar
parallel flows.

We use $L=200d_{0}$, $s=1.4$, and $\eta/\sqrt{m_{0}k}=1$, where $m_{0}$ is the
mass of a particle with diameter $d_{0}$. The critical packing fraction is
estimated as $\phi_{\rm c}\simeq 0.843$ in Appendix A, below which the
particles cannot form the jammed solid-like state. The parallel flow for
$\phi_{0}<\phi_{\rm c}$ has been studied in Refs. [52, 53, 54], but we focus
on dense systems with $\phi_{0}\geq\phi_{\rm c}$. The time evolution equation
is solved numerically using LAMMPS (the open source molecular dynamics program
from Sandia National Laboratories) with the time step $\Delta
t=0.05\sqrt{m_{0}/k}$ [55, 56]. The units of the velocity, the external force,
and the mass flux are $v_{0}=d_{0}/\sqrt{m_{0}/k}$, $f_{0}=kd_{0}/m_{0}$, and
$Q_{0}=\sqrt{m_{0}k}$, respectively.

##  III Flow near jamming

Figure 2 displays the velocity profile $u(z)$ for various values of $f$ with
$\phi_{0}=0.850$ and $H=400d_{0}$. Here, $u(z)$ is obtained as the average of
the particle velocity $v_{i,x}=\frac{d}{dt}x_{i}$ at $z$. The velocity profile
is symmetric about $z=0$, and $u(z)=0$ at the rough plates ($z=\pm H/2$). For
$f=1.4\times 10^{-6}f_{0}$ and $f=1.3\times 10^{-6}f_{0}$, a plug flow with a
constant velocity is visible near $z=0$. Near the plates, fluidized regions
with nonzero shear rate are observed, where the velocity profile appears
parabolic. For $f=1.2\times 10^{-6}f_{0}$ and $f=1.0\times 10^{-6}f_{0}$, the
velocity $u(z)$ is zero in the entire system. This indicates that the
particles are jammed for small $f$.

Figure 2: Velocity profile $u(z)$ with $\phi_{0}=0.850$ and $H=400d_{0}$ for
various values of $f$.

In Fig. 3, the mass flux given by $Q=\sum_{i}m_{i}v_{i,x}/L$ is plotted
against $f$ for various values of $\phi_{0}$ and $H$. The mass flux $Q$
increases with $H$ and decreases as $\phi_{0}$ increases. For any $\phi_{0}$
and $H$, $Q$ increases from $0$ as $f$ exceeds a threshold $f_{\rm c}$. The
critical force $f_{\rm c}$ depending on $\phi_{0}$ and $H$ distinguishes the
jammed state with $Q=0$ from the unjammed state with $Q>0$.

Figure 3:  Mass flux $Q$ against $f$ for various values of $\phi_{0}$ with
$H=300d_{0}$ (filled symbols) and $H=400d_{0}$ (open symbols).

In Fig. 5, the critical force $f_{\rm c}$ is plotted against $\phi_{0}$ for
various values of $H$. We estimate $f_{\rm c}$ as the force $f$ at which $Q$
exceeds $Q_{\rm th}=1.0\times 10^{-4}Q_{0}(H/d_{0})^{2}$. Note that $f_{\rm
c}$ is unchanged if we choose a smaller $Q_{\rm th}$. The critical force
$f_{\rm c}$ increases as $\phi_{0}$ exceeds $\phi_{\rm c}\simeq 0.843$ and
decreases as $H$ increases. This indicates that the system jams as $H$
decreases or $\phi_{0}$ increases.

Figure 4:  Critical force $f_{c}$ against $\phi_{0}$ for various values of
$H$.

##  IV Scaling laws near jamming

###  IV.1 Continuum equations

In the steady flow shown in Fig. 2, the pressure $p(z)$, the shear stress
$\sigma(z)$, and the packing fraction $\phi(z)$ at $z$ satisfy the momentum
conservation [57, 52]

$\displaystyle\dfrac{d\sigma(z)}{dz}=-\rho_{s}\phi(z)f,$ (5)
$\displaystyle\dfrac{dp(z)}{dz}=0.$ (6)

We adopt the $\mu(I)$-rheology [25, 39] as the constitutive equation, where
the local inertia number $I(z)=\dot{\gamma}(z)d/\sqrt{p(z)/\rho_{s}}$ and bulk
friction coefficient $\mu(z)=|\sigma(z)|/p(z)$ satisfy

$\displaystyle I(z)=\begin{cases}\mathcal{I}(\mu(z))&(\mu(z)>\mu_{\rm s})\\\
0&(\mu(z)\leq\mu_{\rm s}),\end{cases}$ (7)

with

$\displaystyle\mathcal{I}(\mu)=I_{0}\dfrac{\mu-\mu_{\rm s}}{\mu_{2}-\mu},$ (8)

and the shear rate $\dot{\gamma}(z)=|du(z)/dz|$. Here, $\mu_{\rm s}$ is the
static yield coefficient, $\mu_{2}$ is the maximum value of $\mu$, and $I_{0}$
is a dimensionless parameter characterizing the non-linear response. This
constitutive equation indicates that the yield stress depends on the pressure,
as $\sigma_{\rm Y}=\mu_{\rm s}p$. We assume that the local pressure $p(z)$
depends on the shear rate $\dot{\gamma}(z)$ and packing fraction $\phi(z)$:
[10, 14, 15, 16, 58]

$\displaystyle p(z)=\mathcal{P}(\dot{\gamma}(z),\phi(z)).$ (9)

The function $\mathcal{P}(\dot{\gamma},\phi)$ satisfies

$\displaystyle\lim_{\dot{\gamma}\to
0}\mathcal{P}(\dot{\gamma},\phi)=\Pi(\phi),$ (10)

where $\Pi(\phi)$ is the pressure in the quasi-static limit, $\dot{\gamma}\to
0$. For particles with the linear repulsive interaction (Eq. (4)), $\Pi(\phi)$
obeys a scaling law [5]

$\displaystyle\lim_{\phi\to\phi_{\rm c}}\Pi(\phi)=B|\phi-\phi_{\rm c}|$ (11)

for $\phi\geq\phi_{\rm c}$ with a constant $B$, which is numerically confirmed
in our system as shown in Appendix A. Here, we have assumed that $\sigma(z)$
and $p(z)$ are determined by the local quantities $\dot{\gamma}(z)$ and
$\phi(z)$, as the non-local effect is irrelevant for sufficiently large
systems [45]. In fact, we have verified that the same analytical results are
obtained for $H/d\gg 1$ even when we apply the non-local constitutive equation
used in Ref. [45]. We have also verified that different forms of
$\mathcal{I}(\mu)$ proposed in Refs. [40, 41, 42, 43, 44] give identical
scaling laws for $f_{\rm c}$ and $Q$ shown below.

The mean packing fraction $\phi_{0}$ satisfies

$\displaystyle\phi_{0}=\int_{-H/2}^{H/2}dz\ \phi(z)/H.$ (12)

The boundary condition is given by

$\displaystyle u(z=\pm H/2)=0.$ (13)

The system is symmetric with respect to $z=0$:

$\displaystyle u(z)=u(-z),$ (14) $\displaystyle\sigma(z)=-\sigma(-z).$ (15)

The mass flux $Q$ is given by

$\displaystyle Q=\int_{-H/2}^{H/2}dz\ \rho_{s}\phi(z)u(z).$ (16)

###  IV.2 Critical force

First, we consider the jammed state with $\dot{\gamma}(z)=0$ for $f<f_{\rm
c}$. From Eqs. (6), (9), (10), and (12) with $\dot{\gamma}(z)=0$, $\phi(z)$
and $p(z)$ are given by

$\displaystyle\phi(z)=\phi_{0},$ (17) $\displaystyle p(z)=\Pi(\phi_{0}),$ (18)

respectively. Substituting Eq. (17) into Eq. (5) and using Eq. (15), we obtain

$\displaystyle\sigma(z)=-\rho_{s}\phi_{0}fz.$ (19)

From Eqs. (18) and (19), the bulk friction $\mu(z)=|\sigma(z)|/p(z)$ is given
by

$\displaystyle\mu(z)=\rho_{s}\phi_{0}f|z|/\Pi(\phi_{0}).$ (20)

According to the constitutive equation, Eq. (7), the jammed state with
$\dot{\gamma}(z)=0$ is realized when $\mu(z)\leq\mu_{\rm s}$ is satisfied for
$-H/2\leq z\leq H/2$. Using Eq. (20), this condition is replaced with $f\leq
f_{\rm c}$, where

$\displaystyle f_{\rm c}=\dfrac{2\mu_{\rm s}\Pi(\phi_{0})}{\rho_{\rm
s}\phi_{0}H}.$ (21)

For $\phi_{0}-\phi_{\rm c}\ll 1$, Eqs. (11) and (21) give the scaling law

$\displaystyle f_{\rm
c}H=\left(f_{0}d_{0}\right)\mathcal{F}\left(\phi_{0}-\phi_{\rm c}\right),$
(22)

where $\mathcal{F}(\xi)$ denotes a scaling function. In our model,
$\mathcal{F}(\xi)$ is given by $\mathcal{F}(\xi)=2\mu_{\rm s}B\xi/(\rho_{\rm
s}f_{0}\phi_{\rm c}d_{0})$. In Fig. 5, we present the scaling plot based on
Eq. (22), which demonstrates the collapse of all data points across various
values of $\phi_{0}$ and $H$. This consolidation confirms the validity of the
scaling law expressed in Eq. (22). where the all the data with various
$\phi_{0}$ and $H$ are nicely collapsed. This confirms the validity of the
scaling law, Eq. (22).

Figure 5: Scaled plot of $f_{\rm c}$ based on Eq. (22) for various values of
$H$.

###  IV.3 Mass flux

Next, we focus on the critical behavior of the mass flux $Q$ near $f_{\rm c}$.
By introducing a reduced force

$\displaystyle\epsilon=\dfrac{f-f_{c}}{f_{c}},$ (23)

we perturbatively obtain $\phi(z)$, $u(z)$, and $Q$ for $\epsilon\ll 1$ from
Eqs. (5)–(15). The details of the perturbation theory are shown in Appendix B.
In the lowest order approximation, $\phi(z)$ and $p(z)$ are the same as the
solutions for $\dot{\gamma}=0$ given by Eqs. (17) and (18). The local shear
rate $\dot{\gamma}(z)$ is given by

$\displaystyle\dot{\gamma}(z)=\begin{cases}\dfrac{\mu_{\rm
s}I_{0}\sqrt{\Pi(\phi_{0})/\rho_{\rm s}}}{(\mu_{\rm
s}-\mu_{2})d}\epsilon\zeta(z)&(|z|>H/2-\ell)\\\ 0&(|z|\leq
H/2-\ell),\end{cases}$ (24)

with $\zeta(z)=(|z|-H/2+\ell)/\ell$. The area for $|z|>H/2-\ell$ represents
the fluidized region with size $\ell$ as shown in Fig. 1, and $\zeta(z)$
represents the normalized position in the fluidized region. Size $\ell$ of the
fluidized region is given by $\ell=H\epsilon/2$ for $\epsilon\ll 1$.
Integrating the shear rate $\dot{\gamma}(z)$ with respect to $z$ and using the
boundary condition, Eq. (13), we obtain the velocity profile for $\epsilon\ll
1$ as

$\displaystyle u(z)=\begin{cases}U_{\rm
M}\left|1-\\{\zeta(z)\\}^{2}\right|&(|z|>H/2-\ell)\\\ U_{\rm M}&(|z|\leq
H/2-\ell),\end{cases}$ (25)

with the maximum velocity

$\displaystyle U_{\rm M}=\dfrac{\mu_{\rm s}I_{0}}{4(\mu_{\rm
s}-\mu_{2})}\dfrac{\sqrt{\Pi(\phi_{0})/\rho_{\rm s}}H}{d}\epsilon^{2}.$ (26)

In the fluidized region for $|z|>H/2-\ell$, the granular materials exhibit a
parabolic velocity profile, while a plug flow appears for $|z|\leq H/2-\ell$.
This velocity profile is similar to that in yield stress fluids [59, 60, 61].
The velocity profile given by Eq. (25) is consistent with that in Fig. 2. A
similar plug flow is obtained for granular materials with $\phi_{0}<\phi_{c}$,
but the critical force does not exist in relatively dilute systems [52, 53,
54].

The mass flux $Q$ is obtained by substituting the solutions of $\phi(z)$ and
$u(z)$ into Eq. (16) as

$\displaystyle Q=\rho_{s}\phi_{0}HU_{\rm M}$ (27)

for $\epsilon\ll 1$. Substituting Eqs. (11), (21), (23), and (26) into this
equation, we obtain a critical scaling law near $\phi_{\rm c}$ for the mass
flux $Q$ as

$\displaystyle\dfrac{Q}{H^{2}\sqrt{\phi_{0}-\phi_{\rm
c}}}=(Q_{0}/d_{0}^{2})\mathcal{Q}\left(\dfrac{fH/(f_{0}d_{0})}{\phi_{0}-\phi_{\rm
c}}\right),$ (28)

with a scaling function $\mathcal{Q}(\xi)$. In our analytical model,
$\mathcal{Q}(\xi)$ is given by $\mathcal{Q}(\xi)=\alpha(\beta\xi-1)^{2}$, with
$\alpha=\dfrac{\mu_{\rm s}I_{0}\phi_{0}\sqrt{\rho_{\rm
s}B}d_{0}^{2}}{4(\mu_{\rm s}-\mu_{2})dQ_{0}}$ and $\beta=\dfrac{\rho_{\rm
s}\phi_{0}f_{0}d_{0}}{2\mu_{\rm s}B}$.

In Fig. 6, we demonstrate the scaling plot of the mass flux $Q$ based on Eq.
(28). Although the system is inhomogeneous, with the shear rate depending on
$z$ as shown in Fig. 2, we find excellent consolidation of data in Fig. 6,
which mirrors critical scaling laws observed in homogeneous systems near
jamming [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 24].

Figure 6: Scaling plot of the mass flux $Q$ based on Eq. (28) for various
values of $\phi_{0}$ and $H$. The legends represent $(\phi_{0},~{}H)$.

##  V Conclusion and discussion

We have numerically studied the critical behavior of frictionless granular
materials between rough parallel plates under an external force. Jamming is
observed as the transition from the jammed state with the mass flux $Q=0$ to
the unjammed state with $Q>0$ as $f$ exceeds the critical force $f_{\rm c}$.
Based on a continuum model, we have analytically obtained the velocity profile
$u(z)$, which gives the scaling laws for $f_{\rm c}$ and $Q$. The scaling laws
have been numerically verified. This indicates that critical behaviors similar
to those in homogeneous systems can be observed in the inhomogeneous system.

In our analysis, we have considered the $z$-dependence of the packing fraction
$\phi(z)$, but it is represented as a constant in the lowest order
approximation. In previous studies based on continuum models [34, 45, 31, 51],
the packing fraction is assumed as uniform, but we have confirmed that this
assumption is justified for $\epsilon=(f-f_{\rm c})/f_{\rm c}\ll 1$ with
$\phi_{0}\geq\phi_{\rm c}$. For relatively dilute systems with
$\phi_{0}<\phi_{\rm c}$, the critical force $f_{\rm c}$ is zero, and the
$z$-dependence of the packing fraction $\phi(z)$ is essential [52, 53, 54].

We have considered frictionless particles with a repulsive interaction.
However, we expect that critical behaviors are also observed in particles with
friction or attractive interactions, as the $\mu(I)$ rheology is applicable in
these systems [34, 62, 63]. An extension of our theory to these systems will
be our future work.

###### Acknowledgements.

The authors thank R. Kuroda, K. Miyazaki, T. Uneyama, N. Oyama, H. Hayakawa,
K. Saitoh, S. Takada, and T. Barker for fruitful discussions. M.O. is
partially supported by JSPS KAKENHI (Grant Nos. JP21H01006 and JP23K03248).
K.Y. is partially supported by JSPS Fellows (Grant No. 21J13720). K.Y. would
like to thank the Research center for Computational Science, Okazaki, Japan
for making its supercomputer system available (Project: 22-IMS-C267 and
23-IMS-C126). We would like to thank Editage (www.editage.jp) for English
language editing.

The analytical calculation is performed by K.H and M.O. The numerical
simulations using LAMMPS are performed by K.Y. The first draft of the
manuscript was written by M.O.

##  Appendix A Pressure in the jammed state

In this section, we present the numerical results for the pressure $P$, which
is given by $P=\sum_{i}\sum_{j>i}(x_{ij}F_{ij,x}+y_{ij}F_{ij,y})/(2LH)$ with
$\bm{r}_{ij}=(x_{ij},y_{ij})$ and $\bm{F}_{ij}=(F_{ij,x},F_{ij,y})$ . In Fig.
7, we plot the pressure $P$ against $\phi_{0}-\phi_{\rm c}$ for $f=0$. We
estimate the critical packing fraction $\phi_{\rm c}$ as $\phi_{0}$ where $P$
exceeds a threshold $P_{\rm th}=1.0\times 10^{-5}f_{0}/d_{0}$ because the
nonzero pressure is considered the indicator of jamming [5]. The critical
fraction $\phi_{\rm c}$ has a slight dependence on $H$ as $\phi_{\rm
c}=0.8425,0.8430,$ and $0.8430$ for $H/d_{0}=200,300$, and $400$,
respectively. The pressure $P$ is almost proportional to $\phi_{0}-\phi_{\rm
c}$ in Fig. 7, which is consistent with Eq. (10).

Figure 7: Pressure $P$ at the rough plates with $f=0$ against
$\phi_{0}-\phi_{\rm c}$ for various values of $H$. The dashed line represents
Eq. (11) with $B=0.4$.

##  Appendix B Details of analysis

In this section, we demonstrate the details of our perturbation theory. From
Eqs. (5) and (6) with Eq. (15), the stress and pressure are given by

$\displaystyle\sigma(z)=-\rho_{s}f\int_{0}^{z}dz^{\prime}\ \phi(z^{\prime}),$
(29) $\displaystyle p(z)=P_{\rm b},$ (30)

where $P_{\rm b}$ is the pressure at the plates. Near $z=0$, there is a region
where $|\sigma|$ is less than $\mu_{\rm s}P_{\rm b}$. According to Eq. (7),
the shear rate $\dot{\gamma}$ is zero in this region. The area for $|z|\leq
H/2-\ell$ corresponds to the plug region with $\dot{\gamma}=0$, while the area
for $|z|>H/2-\ell$ represents the fluidized region with $\dot{\gamma}>0$ with
the size $\ell$.

With the introduction of the nondimensionalized position

$\displaystyle\zeta(z)=(|z|-H/2+\ell)/\ell,$ (31)

we represent the shear rate $\dot{\gamma}(z)$ as

$\displaystyle\dot{\gamma}(z)=\begin{cases}\Gamma(\zeta(z))&(|z|>H/2-\ell)\\\
0&(|z|\leq H/2-\ell),\end{cases}$ (32)

where $\Gamma(\zeta)$ is the shear rate in the fluidized region satisfying

$\displaystyle\lim_{\zeta\to 0}\Gamma(\zeta)=0.$ (33)

According to Eqs. (9), (10), (30), and (32), the packing fraction $\phi(z)$ is
constant in the plug region. Hence, $\phi(z)$ is represented as

$\displaystyle\phi(z)=\begin{cases}\phi_{\rm
S}\left\\{1-\mathcal{F}(\zeta(z))\right\\}&(|z|>H/2-\ell)\\\ \phi_{\rm
S}&(|z|\leq H/2-\ell)\end{cases}$ (34)

with the packing fraction $\phi_{\rm S}$ in the plug region and a function
$\mathcal{F}(\zeta)$ satisfying

$\displaystyle\lim_{\zeta\to 0}\mathcal{F}(\zeta)=0.$ (35)

At the boundary of the plug region, $z=H/2-\ell$, the shear stress $\sigma$
satisfy $|\sigma|=\mu_{s}P_{\rm b}$. Using Eqs. (29), (30), and (34), this
condition is replaced by

$\displaystyle\mu_{s}P_{\rm b}=\rho_{s}\phi_{\rm
S}f\left(\dfrac{H}{2}-\ell\right).$ (36)

Substituting (34) into Eq. (12), we obtain

$\displaystyle\phi_{0}=\phi_{\rm
S}\left\\{1-\dfrac{2\ell}{H}\int_{0}^{1}d\zeta\mathcal{F}(\zeta)\right\\}.$
(37)

From Eqs. (29), (30), (34), and (36), we derive the bulk friction coefficient
$\mu=|\sigma|/p$ at $\zeta\geq 0$ as

$\displaystyle\mu(\zeta)-\mu_{s}=\mu_{s}\dfrac{\ell}{H/2-\ell}\left\\{\zeta-\int_{0}^{\zeta}d\zeta^{\prime}\
\mathcal{F}(\zeta^{\prime})\right\\}.$ (38)

Substituting Eqs. (30) and (32) into Eq. (7) with
$I=\dot{\gamma}d/\sqrt{p/\rho_{\rm s}}$, we obtain

$\displaystyle\dfrac{\Gamma(\zeta)d}{\sqrt{P_{\rm
b}/\rho_{s}}}=\mathcal{I}\left(\mu(\zeta)-\mu_{s}\right).$ (39)

Using Eqs. (9) and (10) with Eqs. (30), (32), and (34), we find

$\displaystyle P_{\rm b}=\Pi(\phi_{\rm S}),$ (40) $\displaystyle
0=\Delta\left(\phi_{\rm
S}\left\\{1-\mathcal{F}(\zeta)\right\\},\Gamma(\zeta)\right),$ (41)

with $\Delta(\phi,\dot{\gamma})=P(\phi,\dot{\gamma})-\Pi(\phi)$, which
satisfies $\lim_{\dot{\gamma}\to 0}\Delta(\phi,\dot{\gamma})=0$. From Eqs.
(36)-(41), the variables $\ell$, $\phi_{\rm S}$, $P_{\rm b}$, $\Gamma(\zeta)$,
and $\mathcal{F}(\zeta)$ are determined.

For $f\simeq f_{\rm c}$ with $\epsilon\ll 1$, $\ell$, $\phi_{\rm S}$, $P_{\rm
b}$, $\Gamma(\zeta)$, and $\mathcal{F}(\zeta)$ are expanded as

$\displaystyle\ell$ $\displaystyle=\epsilon\ell^{(1)}+O(\epsilon^{2}),$ (42)
$\displaystyle\phi_{\rm S}$ $\displaystyle=\phi_{0}+\epsilon\phi_{\rm
S}^{(1)}+O(\epsilon^{2}),$ (43) $\displaystyle P_{\rm b}$
$\displaystyle=\Pi(\phi_{0})+\epsilon P_{\rm b}^{(1)}+O(\epsilon^{2}),$ (44)
$\displaystyle\Gamma(\zeta)$
$\displaystyle=\epsilon\Gamma^{(1)}(\zeta)+O(\epsilon^{2}),$ (45)
$\displaystyle\mathcal{F}(\zeta)$
$\displaystyle=\epsilon\mathcal{F}^{(1)}(\zeta)+O(\epsilon^{2}).$ (46)

Substituting these equations into Eqs. (36)-(41) and extract all terms
proportional to $\epsilon$, we obtain

$\displaystyle\ell^{(1)}$ $\displaystyle=H/2,$ (47) $\displaystyle\phi_{\rm
S}^{(1)}$ $\displaystyle=0,$ (48) $\displaystyle P_{\rm b}^{(1)}$
$\displaystyle=0,$ (49) $\displaystyle\Gamma^{(1)}(\zeta)$
$\displaystyle=\dfrac{\mu_{\rm s}I_{0}\sqrt{\Pi(\phi_{0})}/\rho_{\rm
s}}{(\mu_{2}-\mu_{\rm s})d}\zeta,$ (50)
$\displaystyle\mathcal{F}^{(1)}(\zeta)$
$\displaystyle=\dfrac{C_{1}}{C_{2}+C_{3}}\Gamma^{(1)}(\zeta),$ (51)

with

$\displaystyle C_{1}$
$\displaystyle=\left.\dfrac{\partial\Delta(\phi,\dot{\gamma})}{\partial\dot{\gamma}}\right|_{(\phi,\dot{\gamma})=(\phi_{0},0)},$
(52) $\displaystyle C_{2}$
$\displaystyle=\phi_{0}\left.\dfrac{d\Pi(\phi)}{d\phi}\right|_{\phi=\phi_{0}},$
(53) $\displaystyle C_{3}$
$\displaystyle=\left.\phi_{0}\dfrac{\partial\Delta(\phi,\dot{\gamma})}{\partial\phi}\right|_{(\phi,\dot{\gamma})=(\phi_{0},0)}.$
(54)

Substituting these equations into Eqs. (42)-(46), we obtain

$\displaystyle\ell$ $\displaystyle=\dfrac{H}{2}\epsilon+O(\epsilon^{2}),$ (55)
$\displaystyle\phi(z)$ $\displaystyle=\phi_{0}+O(\epsilon),$ (56)
$\displaystyle p(z)$ $\displaystyle=\Pi(\phi_{0})+O(\epsilon^{2}),$ (57)
$\displaystyle\dot{\gamma}(z)$ $\displaystyle=\begin{cases}\dfrac{\mu_{\rm
s}I_{0}\sqrt{\Pi(\phi_{0})/\rho_{\rm s}}}{(\mu_{\rm
s}-\mu_{2})d}\epsilon\zeta(z)+O(\epsilon^{2})&(|z|>H/2-\ell)\\\ 0&(|z|\leq
H/2-\ell).\end{cases}$ (58)

Neglecting the higher order terms in Eqs. (55)-(58), we obtain $\phi(z)$,
$p(z)$, and $\dot{\gamma}(z)$ as Eqs. (17), (18), and (24), respectively.

Integrating $\dot{\gamma}(z)$ in Eq. (58) with respect to $z$, we obtain

$\displaystyle u(z)=\begin{cases}\tilde{U}_{\rm
M}\epsilon^{2}\left|1-\\{\zeta(z)\\}^{2}\right|+O(\epsilon^{3})&(|z|>H/2-\ell)\\\
\tilde{U}_{\rm M}\epsilon^{2}+O(\epsilon^{3})&(|z|\leq H/2-\ell)\end{cases}$
(59)

with

$\displaystyle\tilde{U}_{\rm M}=\dfrac{\mu_{\rm s}I_{0}}{4(\mu_{\rm
s}-\mu_{2})}\dfrac{\sqrt{\Pi(\phi_{0})/\rho_{\rm s}}H}{d}.$ (60)

Neglecting $O(\epsilon^{3})$ in Eq. (59), we obtain Eq. (25). Subsituting Eqs.
(56), (59), and (55) into Eq. (16), we obtain

$\displaystyle Q=\rho_{\rm s}\phi_{0}\tilde{U}_{\rm
M}H\epsilon^{2}+O(\epsilon^{3}).$ (61)

Neglecting $O(\epsilon^{3})$ in Eq. (61), we obtain Eq. (27).

## References

  * van Hecke [2010] M. van Hecke, Jamming of soft particles: geometry, mechanics, scaling and isostaticity, J. Phys. Condens. Matter 22, 033101 (2010).
  * Behringer and Chakraborty [2019] R. P. Behringer and B. Chakraborty, The physics of jamming for granular materials: a review., Rep. Prog. Phys. 82, 012601 (2019).
  * Bonn  _et al._ [2017] D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Yield stress materials in soft condensed matter, Rev. Mod. Phys. 89, 035005 (2017).
  * Liu and Nagel [1998] A. J. Liu and S. R. Nagel, Jamming is not just cool any more, Nature 396, 21 (1998).
  * O’Hern  _et al._ [2002] C. S. O’Hern, S. A. Langer, A. J. Liu, and S. R. Nagel, Random packings of frictionless particles, Phys. Rev. Lett. 88, 075507 (2002).
  * O’Hern  _et al._ [2003] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Jamming at zero temperature and zero applied stress: The epitome of disorder, Phys. Rev. E 68, 011306 (2003).
  * Wyart [2005] M. Wyart, On the rigidity of amorphous solids, Ann. Phys. 30, 1 (2005).
  * Zaccone and Scossa-Romano [2011] A. Zaccone and E. Scossa-Romano, Approximate analytical description of the nonaffine response of amorphous solids, Phys. Rev. B 83, 184205 (2011).
  * Wyart  _et al._ [2005] M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Witten, Effects of compression on the vibrational modes of marginally jammed solids, Phys. Rev. E 72, 051306 (2005).
  * Hatano  _et al._ [2007] T. Hatano, M. Otsuki, and S. ichi Sasa, Criticality and scaling relations in a sheared granular material, J. Phys. Soc. Japan 76, 023001 (2007).
  * Ikeda  _et al._ [2012] A. Ikeda, L. Berthier, and P. Sollich, Unified study of glass and jamming rheology in soft particle systems, Phys. Rev. Lett. 109, 018301 (2012).
  * Kawasaki  _et al._ [2015] T. Kawasaki, D. Coslovich, A. Ikeda, and L. Berthier, Diverging viscosity and soft granular rheology in non-brownian suspensions, Phys. Rev. E 91, 012203 (2015).
  * Olsson and Teitel [2007] P. Olsson and S. Teitel, Critical scaling of shear viscosity at the jamming transition, Phys. Rev. Lett. 99, 178001 (2007).
  * Hatano [2008] T. Hatano, Scaling properties of granular rheology near the jamming transition, J. Phys. Soc. Japan 77, 123002 (2008).
  * Otsuki and Hayakawa [2009a] M. Otsuki and H. Hayakawa, Universal scaling for the jamming transition, Prog. Theor. Phys. 121, 647 (2009a).
  * Otsuki and Hayakawa [2009b] M. Otsuki and H. Hayakawa, Critical behaviors of sheared frictionless granular materials near the jamming transition, Phys. Rev. E 80, 011308 (2009b).
  * Tighe  _et al._ [2010] B. P. Tighe, E. Woldhuis, J. J. C. Remmers, W. van Saarloos, and M. van Hecke, Model for the scaling of stresses and fluctuations in flows near jamming, Phys. Rev. Lett. 105, 088303 (2010).
  * Nordstrom  _et al._ [2010] K. N. Nordstrom, E. Verneuil, P. E. Arratia, A. Basu, Z. Zhang, A. G. Yodh, J. P. Gollub, and D. J. Durian, Microfluidic rheology of soft colloids above and below jamming, Phys. Rev. Lett. 105, 175701 (2010).
  * Tighe [2011] B. P. Tighe, Relaxations and rheology near jamming, Phys. Rev. Lett. 107, 158303 (2011).
  * Otsuki and Hayakawa [2014] M. Otsuki and H. Hayakawa, Avalanche contribution to shear modulus of granular materials, Phys. Rev. E 90, 042202 (2014).
  * Otsuki and Hayakawa [2017] M. Otsuki and H. Hayakawa, Discontinuous change of shear modulus for frictional jammed granular materials, Phys. Rev. E 95, 062902 (2017).
  * Saitoh and Kawasaki [2020] K. Saitoh and T. Kawasaki, Critical scaling of diffusion coefficients and size of rigid clusters of soft athermal particles under shear, Front. Phys. 8, 99 (2020).
  * Otsuki and Hayakawa [2012] M. Otsuki and H. Hayakawa, Critical scaling of a jammed system after a quench of temperature, Phys. Rev. E 86, 031505 (2012).
  * Saitoh  _et al._ [2020] K. Saitoh, T. Hatano, A. Ikeda, and B. P. Tighe, Stress relaxation above and below the jamming transition, Phys. Rev. Lett. 124, 118001 (2020).
  * MiDi [2004] G. D. R. MiDi, On dense granular flows, Eur. Phys. J. E 14, 341 (2004).
  * Zuriguel [2014] I. Zuriguel, Invited review: Clogging of granular materials in bottlenecks, Pap. Phys. 6, 060014 (2014).
  * Cai  _et al._ [2021] G. Cai, A. B. Harada, and K. Nordstrom, Mesoscale metrics on approach to the clogging point, Granul. Matter 23, 69 (2021).
  * Pouliquen [1999] O. Pouliquen, On the shape of granular fronts down rough inclined planes, Phys. Fluids 11, 1956 (1999).
  * Silbert  _et al._ [2003] L. E. Silbert, J. W. Landry, and G. S. Grest, Granular flow down a rough inclined plane: Transition between thin and thick piles, Phys. Fluids 15, 1 (2003).
  * Forterre and Pouliquen [2002] Y. Forterre and O. Pouliquen, Stability analysis of rapid granular chute flows: formation of longitudinal vortices, J. Fluid Mech. 467, 361 (2002).
  * Liu and Henann [2017] D. Liu and D. L. Henann, Non-local continuum modelling of steady, dense granular heap flows, J. Fluid Mech. 831, 212 (2017).
  * Lemieux and Durian [2000] P.-A. Lemieux and D. J. Durian, From avalanches to fluid flow: A continuous picture of grain dynamics down a heap, Phys. Rev. Lett. 85, 4273 (2000).
  * Komatsu  _et al._ [2001] T. S. Komatsu, S. Inagaki, N. Nakagawa, and S. Nasuno, Creep motion in a granular pile exhibiting steady surface flow, Phys. Rev. Lett. 86, 1757 (2001).
  * Jop  _et al._ [2006] P. Jop, Y. Forterre, and O. Pouliquen, A constitutive law for dense granular flows, Nature 441, 727 (2006).
  * Parker  _et al._ [1997] D. Parker, A. Dijkstra, T. Martin, and J. Seville, Positron emission particle tracking studies of spherical particle motion in rotating drums, Chem. Eng. Sci. 52, 2011 (1997).
  * Pignatel  _et al._ [2012] F. Pignatel, C. Asselin, L. Krieger, I. C. Christov, J. M. Ottino, and R. M. Lueptow, Parameters and scalings for dry and immersed granular flowing layers in rotating tumblers, Phys. Rev. E 86, 011304 (2012).
  * Orpe and Khakhar [2001] A. V. Orpe and D. V. Khakhar, Scaling relations for granular flow in quasi-two-dimensional rotating cylinders, Phys. Rev. E 64, 031302 (2001).
  * Zheng  _et al._ [2019] Q. Zheng, L. Bai, L. Yang, and A. Yu, 110th anniversary: continuum modeling of granular mixing in a rotating drum, Ind. Eng. Chem. 58, 19251 (2019).
  * da Cruz  _et al._ [2005] F. da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir, Rheophysics of dense granular materials: Discrete simulation of plane shear flows, Phys. Rev. E 72, 021309 (2005).
  * Hatano [2007] T. Hatano, Power-law friction in closely packed granular materials, Phys. Rev. E 75, 060301 (2007).
  * Peyneau and Roux [2008] P.-E. Peyneau and J.-N. Roux, Frictionless bead packs have macroscopic friction, but no dilatancy, Phys. Rev. E 78, 011307 (2008).
  * Émilien Azéma  _et al._ [2018] Émilien Azéma, F. Radjaï, and J.-N. Roux, Inertial shear flow of assemblies of frictionless polygons: Rheology and microstructure, Eur. Phys. J. E 41, 2 (2018).
  * Man  _et al._ [2023] T. Man, P. Zhang, Z. Ge, S. A. Galindo-Torres, and K. M. Hill, Friction-dependent rheology of dry granular systems, Acta Mech. Sin. 39, 722191 (2023).
  * Bouzid  _et al._ [2013] M. Bouzid, M. Trulsson, P. Claudin, E. Clément, and B. Andreotti, Nonlocal rheology of granular flows across yield conditions, Phys. Rev. Lett. 111, 238301 (2013).
  * Kamrin and Koval [2012] K. Kamrin and G. Koval, Nonlocal constitutive relation for steady granular flow, Phys. Rev. Lett. 108, 178301 (2012).
  * Henann and Kamrin [2013] D. L. Henann and K. Kamrin, A predictive, size-dependent continuum model for dense granular flows, Proc. Natl. Acad. Sci. U.S.A 110, 6730 (2013).
  * Bouzid  _et al._ [2015] M. Bouzid, A. Izzet, M. Trulsson, E. Clément, P. Claudin, and B. Andreotti, Non-local rheology in dense granular flows, Eur. Phys. J. E 38, 125 (2015).
  * Zhang and Kamrin [2017] Q. Zhang and K. Kamrin, Microscopic description of the granular fluidity field in nonlocal flow modeling, Phys. Rev. Lett. 118, 058001 (2017).
  * Kim and Kamrin [2020] S. Kim and K. Kamrin, Power-law scaling in granular rheology across flow geometries, Phys. Rev. Lett. 125, 088002 (2020).
  * Kim and Kamrin [2023] S. Kim and K. Kamrin, A second-order non-local model for granular flows, Front. Phys. 11, 1092233 (2023).
  * Liu and Henann [2018] D. Liu and D. L. Henann, Size-dependence of the flow threshold in dense granular materials, Soft Matter 14, 5294 (2018).
  * Barker  _et al._ [2022] T. Barker, C. Zhu, and J. Sun, Exact solutions for steady granular flow in vertical chutes and pipes, J. Fluid Mech. 930, A21 (2022).
  * Islam  _et al._ [2022] M. U. Islam, J. T. Jenkins, and S. L. Das, Extended kinetic theory for granular flow in a vertical chute, J. Fluid Mech. 950, A13 (2022).
  * Islam  _et al._ [2023] M. U. Islam, J. T. Jenkins, and S. L. Das, Granular flow through a vertical axisymmetric pipe, Phys. Rev. Fluids 8, L072301 (2023).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117, 1 (1995).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, _et al._ , Lammps-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun. 271, 108171 (2022).
  * Evans and Morriss [2008] D. J. Evans and G. Morriss,  _Statistical Mechanics of Nonequilibrium Liquids_ (Cambridge University Press, 2008).
  * Saitoh and Mizuno [2016] K. Saitoh and H. Mizuno, Anomalous energy cascades in dense granular materials yielding under simple shear deformations, Soft Matter 12, 1360 (2016).
  * Oldroyd [1947] J. G. Oldroyd, Two-dimensional plastic flow of a bingham solid, Math. Proc. Camb. Philos. Soc. 43, 383 (1947).
  * Frigaard  _et al._ [1994] I. A. Frigaard, S. D. Howison, and I. J. Sobey, On the stability of poiseuille flow of a bingham fluid, J. Fluid Mech. 263, 133 (1994).
  * Khatib and Wilson [2001] M. A. Khatib and S. Wilson, The development of poiseuille flow of a yield-stress fluid, J. Fluid Mech. 100, 1 (2001).
  * Vo  _et al._ [2020] T. T. Vo, S. Nezamabadi, P. Mutabaruka, J.-Y. Delenne, and F. Radjai, Additive rheology of complex granular flows, Nat. Commun. 11, 1476 (2020).
  * Roy  _et al._ [2017] S. Roy, S. Luding, and T. Weinhart, A general(ized) local rheology for wet granular materials, New J. Phys. 19, 043014 (2017).


# On $3-2-1$ values of finite multiple harmonic $q$-series at roots of unity

Kh. Hessami Pilehrood  The Fields Institute for Research in Mathematical
Sciences, 222 College St, Toronto, Ontario M5T 3J1 Canada<EMAIL_ADDRESS>,  T. Hessami Pilehrood  The Fields Institute for Research in Mathematical
Sciences, 222 College St, Toronto, Ontario M5T 3J1 Canada<EMAIL_ADDRESS>and  R. Tauraso  Dipartimento di Matematica, Università di Roma “Tor Vergata”,
via della Ricerca Scientifica, 00133 Roma, Italy<EMAIL_ADDRESS>
###### Abstract.

We mainly answer two open questions about finite multiple harmonic $q$-series
on 3-2-1 indices at roots of unity, posed recently by H. Bachmann, Y.
Takeyama, and K. Tasaka. Two conjectures regarding cyclic sums which
generalize the given results are also provided.

###### Key words and phrases:

Multiple harmonic sums, $q$-analogs, roots of unity

###### 2010 Mathematics Subject Classification:

11M32, 11M35, 05A15, 30B10, 30D05

##  1\. Introduction.

For two $r$-tuples of non-negative integers ${\bf s}=(s_{1},\ldots,s_{r})$ and
${\bf t}=(t_{1},\ldots,t_{r})$ and a positive integer $n$, with a complex
number $q$ satisfying $q^{m}\neq 1$ for $n>m>0$, we define two classes of
multiple $q$-harmonic sums

$\displaystyle H_{n}({\bf s};{\bf
t};q)=H_{n}(s_{1},\ldots,s_{r};t_{1},\ldots,t_{r};q)$
$\displaystyle=\sum_{1\leq k_{1}<\cdots<k_{r}\leq
n}\frac{q^{k_{1}t_{1}+\cdots+k_{r}t_{r}}}{[k_{1}]_{q}^{s_{1}}\cdots[k_{r}]_{q}^{s_{r}}},$
$\displaystyle H_{n}^{\star}({\bf s};{\bf
t};q)=H_{n}^{\star}(s_{1},\ldots,s_{r};t_{1},\ldots,t_{r};q)$
$\displaystyle=\sum_{1\leq k_{1}\leq\cdots\leq k_{r}\leq
n}\frac{q^{k_{1}t_{1}+\cdots+k_{r}t_{r}}}{[k_{1}]_{q}^{s_{1}}\cdots[k_{r}]_{q}^{s_{r}}},$

where

$[n]_{q}=\frac{1-q^{n}}{1-q}=1=q+\cdots+q^{n-1}$

is a $q$-analog of positive integer $n$. By convention, we put
$H_{n}(\emptyset)=H_{n}^{\star}(\emptyset)=1$, and $H_{n}({\bf s};{\bf
t};q)=0$ if $n<r$. The number $w({\bf s})=\sum_{j=1}^{r}s_{j}$ is called the
weight of the multiple harmonic sum.

For a primitive $n$-th root of unity $\zeta_{n}$, following work [2], we adopt
the notation

$\displaystyle z_{n}({\bf s};\zeta_{n})$ $\displaystyle=H_{n-1}({\bf s};{\bf
s}-\\{1\\}^{r};\zeta_{n}),$ $\displaystyle z_{n}^{\star}({\bf s};\zeta_{n})$
$\displaystyle=H_{n-1}^{\star}({\bf s};{\bf s}-\\{1\\}^{r};\zeta_{n}),$

where $\\{a\\}^{r}$ denotes an $r$-tuple with $r$ consecutive copies of the
letter $a$ (note that we reversed the order of summation for convenience in
our settings).

In [2], Bachmann, Takeyama, and Tasaka studied special values of
$z_{n}(\\{k\\}^{r};\zeta_{n})$ and in particular for $k=1,2,3$, showed that

$\displaystyle z_{n}(\\{1\\}^{r};\zeta_{n})$
$\displaystyle=\frac{1}{n}\binom{n}{r+1}(1-\zeta_{n})^{r},$ (1) $\displaystyle
z_{n}(\\{2\\}^{r};\zeta_{n})$
$\displaystyle=\frac{(-1)^{r}}{n(r+1)}\binom{n+r}{2r+1}(1-\zeta_{n})^{2r},$
(2) $\displaystyle z_{n}(\\{3\\}^{r};\zeta_{n})$
$\displaystyle=\frac{1}{n^{2}(r+1)}\left(\binom{n+2r+1}{3r+2}+(-1)^{r}\binom{n+r}{3r+2}\right)(1-\zeta_{n})^{3r}.$
(3)

The authors of [2] also formulated two open questions for finite multiple
harmonic $q$-series $z_{n}$ on 3-2-1 indices, namely,

$\begin{split}z_{n}(\\{1\\}^{a},2,\\{1\\}^{b};\zeta_{n})+z_{n}(\\{1\\}^{b},2,\\{1\\}^{a};\zeta_{n})&\stackrel{{\scriptstyle?}}{{=}}-\frac{1}{n}\binom{n+1}{a+b+3}(1-\zeta_{n})^{a+b+2},\\\\[3.0pt]
z_{n}(\\{2\\}^{a},3,\\{2\\}^{b};\zeta_{n})+z_{n}(\\{2\\}^{b},3,\\{2\\}^{a};\zeta_{n})&\stackrel{{\scriptstyle?}}{{=}}\frac{(-1)^{a+b}}{n(a+b+2)}\binom{n+a+b+1}{2(a+b)+3}(1-\zeta_{n})^{2(a+b)+3}.\end{split}$
(4)

In this paper, we prove the above relations and obtain related formulas for
corresponding values of $\xi({\bf s})$, which are defined as the limit values
(see [1, Thm. 1.2])

$\xi({\bf s})=\lim_{n\to\infty}z_{n}({\bf s};e^{\frac{2\pi i}{n}}).$

Note that when $n$ is a prime, formulas (4) as well as (2) and (3) follow from
our results on $q$-congruences for multiple $q$-harmonic sums [3, Thm. 4.1,
Thm. 5.1, Thm. 6.1, and Thm. 8.3], while formula (1) follows from [4, Cor.
2.2]. The methods of our paper [3] can be easily adjusted to prove (4) for
arbitrary positive integer $n$.

######  Theorem 1.1.

For all non-negative integers $a,b$ and any $n$-th primitive root of unity
$\zeta_{n}$,

$z_{n}(\\{2\\}^{a},3,\\{2\\}^{b};\zeta_{n})+z_{n}(\\{2\\}^{b},3,\\{2\\}^{a};\zeta_{n})=\frac{(-1)^{a+b}}{n(a+b+2)}\binom{n+a+b+1}{2(a+b)+3}(1-\zeta_{n})^{2(a+b)+3}.$

######  Theorem 1.2.

For all non-negative integers $a,b$ and any $n$-th primitive root of unity
$\zeta_{n}$,

$z_{n}(\\{1\\}^{a},2,\\{1\\}^{b};\zeta_{n})+z_{n}(\\{1\\}^{b},2,\\{1\\}^{a};\zeta_{n})=-\frac{1}{n}\binom{n+1}{a+b+3}(1-\zeta_{n})^{a+b+2}.$

The complex numbers $\xi({\bf s})$ are of interest in view of their
connections to the finite and symmetric multiple zeta values as was shown in
[1]. After letting $\zeta_{n}=e^{\frac{2\pi i}{n}}$ in Theorem 1.1 and Theorem
1.2, and by noting that for $j,k\in\mathbb{N}$,

$\lim_{n\to\infty}\binom{n+j}{k}\frac{k!}{n^{k}}=1,\quad\text{and}\quad\lim_{n\to\infty}n(1-e^{\frac{2\pi
i}{n}})=-2\pi i,$

we obtain the following corollary.

######  Corollary 1.1.

For all non-negative integers $a,b$,

$\xi(\\{1\\}^{a},2,\\{1\\}^{b})+\xi(\\{1\\}^{b},2,\\{1\\}^{a})=-\frac{(-2\pi
i)^{a+b+2}}{(a+b+3)!},$

and

$\xi(\\{2\\}^{a},3,\\{2\\}^{b})+\xi(\\{2\\}^{b},3,\\{2\\}^{a})=0.$

Note that the last relation can also be readily obtained from the definition
of the symmetric multiple zeta values (see, for example, [1, Def. 2.5]).

Finally, we put forward the following conjectures regarding cyclic sums of
multiple $q$-harmonic sums $z_{n}$ at roots of unity, which generalize both of
the theorems above.

######  Conjecture 1.1 (Cyclic-sum).

Let $d_{0},d_{1},\dots,d_{t}$ be non-negative integers. Then

  * (i)

For every integer $n>r$, where $r=\sum_{j=0}^{t}d_{j}+2t$, and any primitive
root of unity $\zeta_{n}$,

$\sum_{j=0}^{t}z_{n}\left(\\{1\\}^{d_{j}},2,\\{1\\}^{d_{j+1}},2,\dots,2,\\{1\\}^{d_{j+t}}\right)=\frac{(-1)^{t}}{n}\binom{n+t}{r+1}(1-\zeta_{n})^{r}.$

  * (ii)

For every integer $n>r$, where $r=\sum_{j=0}^{t}2d_{j}+3t$, and any primitive
root of unity $\zeta_{n}$,

$\sum_{j=0}^{t}z_{n}\left(\\{2\\}^{d_{j}},3,\\{2\\}^{d_{j+1}},3,\dots,3,\\{2\\}^{d_{j+t}}\right)\in(1-\zeta_{n})^{r}{\mathbb{Q}}.$

In both sums above it is understood that $d_{j}=d_{k}$ if $j\equiv k$ modulo
$t+1$.

Note that the case $t=1$ follows from Theorem 1.1 and Theorem 1.2. The case of
arbitrary $t$ when all $d_{j}$ are zeros follows from (2) and (3).

##  2\. Proof of Theorem 1.1.

Let $\overline{\bf s}=(s_{r},s_{r-1},\ldots,s_{1})$ denote the reverse of
${\bf s}=(s_{1},\ldots,s_{r-1},s_{r})$. Then we have the following relations.

######  Lemma 2.1.

Let ${\bf s}=(s_{1},\ldots,s_{r})$ and ${\bf t}=(t_{1},\ldots,t_{r})$ be two
$r$-tuples of non-negative integers, and $\zeta_{n}$ be an $n$-th primitive
root of unity. Then

$\displaystyle H_{n-1}({\bf s};{\bf t};\zeta_{n})$ $\displaystyle=(-1)^{w({\bf
s})}H_{n-1}(\overline{\bf s};\overline{\bf s}-\overline{\bf t};\zeta_{n}),$
$\displaystyle H_{n-1}^{\star}({\bf s};{\bf t};\zeta_{n})$
$\displaystyle=(-1)^{w({\bf s})}H_{n-1}^{\star}(\overline{\bf s};\overline{\bf
s}-\overline{\bf t};\zeta_{n}),$

and in particular,

$\displaystyle z_{n}({\bf s};\zeta_{n})$ $\displaystyle=(-1)^{w({\bf
s})}H_{n-1}(\overline{\bf s};\\{1\\}^{r};\zeta_{n}),$ (5) $\displaystyle
z_{n}^{\star}({\bf s};\zeta_{n})$ $\displaystyle=(-1)^{w({\bf
s})}H_{n-1}^{\star}(\overline{\bf s};\\{1\\}^{r};\zeta_{n}).$

###### Proof.

Replacing each $k_{i}$ by $n-k_{i}$ and reversing the order of summation, we
get

$\begin{split}H_{n-1}({\bf s};{\bf
t};\zeta_{n})&=\sum_{0<k_{1}<\cdots<k_{r}<n}\frac{\zeta_{n}^{t_{1}k_{1}+\ldots+t_{r}k_{r}}}{[k_{1}]_{\zeta_{n}}^{s_{1}}\dots[k_{r}]_{\zeta_{n}}^{s_{r}}}=\sum_{0<n-k_{1}<\cdots<n-k_{r}<n}\frac{\zeta_{n}^{t_{1}(n-k_{1})+\ldots+t_{r}(n-k_{r})}}{[n-k_{1}]_{\zeta_{n}}^{s_{1}}\dots[n-k_{r}]_{\zeta_{n}}^{s_{r}}}\\\
&=\sum_{0<k_{r}<k_{r-1}\cdots<k_{1}<n}\frac{\zeta_{n}^{-t_{1}k_{1}-\ldots-
t_{r}k_{r}}}{[k_{1}]_{\zeta_{n}}^{s_{1}}\dots[k_{r}]_{\zeta_{n}}^{s_{r}}}\times(-1)^{w({\bf
s})}\zeta_{n}^{k_{1}s_{1}+\ldots+k_{r}s_{r}}\\\\[3.0pt] &=(-1)^{w({\bf
s})}H_{n-1}(\overline{\bf s};\overline{\bf s}-\overline{\bf
t};\zeta_{n}),\end{split}$

where we used the identity

$[n-k_{i}]_{\zeta_{n}}=\frac{1-\zeta_{n}^{n-k_{i}}}{1-\zeta_{n}}=\frac{1-\zeta_{n}^{-k_{i}}}{1-\zeta_{n}}=-\zeta_{n}^{-k_{i}}[k_{i}]_{\zeta_{n}}.$
(6)

Setting ${\bf t}={\bf s}-\\{1\\}^{r}$, we get (5). The proofs for the multiple
harmonic star sums are similar. ∎

Proof of Theorem 1.1. We have

$\begin{split}&\qquad
z_{n}(\\{2\\}^{a},3,\\{2\\}^{b};\zeta_{n})+(1-\zeta_{n})z_{n}(\\{2\\}^{a+b+1})\\\\[3.0pt]
&=\sum_{0<k_{1}<\cdots<k_{a}}\frac{\zeta_{n}^{k_{1}+\ldots+k_{a}}}{[k_{1}]_{\zeta_{n}}^{2}\cdots[k_{a}]_{\zeta_{n}}^{2}}\sum_{k_{a}<k_{a+1}<k_{a+2}}\left(\frac{\zeta_{n}^{2k_{a+1}}}{[k_{a+1}]_{\zeta_{n}}^{3}}+\frac{(1-\zeta_{n})\zeta_{n}^{k_{a+1}}}{[k_{a+1}]_{\zeta_{n}}^{2}}\right)\\\\[3.0pt]
&\times\sum_{k_{a+1}<k_{a+2}<\cdots<k_{a+b+1}\leq
n}\frac{\zeta_{n}^{k_{a+2}+\ldots+k_{a+b+1}}}{[k_{a+2}]_{\zeta_{n}}^{2}\cdots[k_{a+b+1}]_{\zeta_{n}}^{2}}\\\\[3.0pt]
&=H_{n-1}(\\{2\\}^{a},3,\\{2\\}^{b};\\{1\\}^{a+b+1};\zeta_{n})=-z_{n}(\\{2\\}^{b},3,\\{2\\}^{a};\zeta_{n}),\end{split}$

where in the last equality we used (5). Hence

$z_{n}(\\{2\\}^{a},3,\\{2\\}^{b};\zeta_{n})+z_{n}(\\{2\\}^{b},3,\\{2\\}^{a};\zeta_{n})=-(1-\zeta_{n})z_{n}(\\{2\\}^{a+b+1}),$

which, by (2), implies the theorem. ∎

##  3\. Proof of Theorem 1.2.

The $q$-binomial coefficient, or Gaussian coefficient, when $q$ is specified
to a primitive root of unity has the following properties.

######  Lemma 3.1.

Let $n>1$ be a positive integer. Then for any primitive $n$-th root of unity
$\zeta_{n}$ and $1\leq k<n$,

$\genfrac{[}{]}{0.0pt}{}{n-1}{k}_{\zeta_{n}}=(-1)^{k}\zeta_{n}^{-\binom{k+1}{2}}.$

###### Proof.

We have

$\genfrac{[}{]}{0.0pt}{}{n-1}{k}_{\zeta_{n}}=\prod_{j=1}^{k}\frac{[n-j]_{\zeta_{n}}}{[j]_{\zeta_{n}}}=\prod_{j=1}^{k}\frac{1-\zeta_{n}^{n-j}}{1-\zeta_{n}^{j}}=\prod_{j=1}^{k}\frac{1-\zeta_{n}^{-j}}{1-\zeta_{n}^{j}}=\prod_{j=1}^{k}(-\zeta_{n}^{-j})=(-1)^{k}\zeta_{n}^{-\binom{k+1}{2}}.$

∎

The proof of Theorem 1.2 is based on the following multiple $q$-binomial
identity.

Theorem A ([3], Thm. 8.1) Let $n,s_{1},\ldots,s_{r}$ be positive integers.
Then

$\sum_{k=1}^{n}\genfrac{[}{]}{0.0pt}{}{n}{k}_{q}(-1)^{k}q^{\binom{k+1}{2}}\sum_{1\leq
k_{1}<k_{2}<\ldots<k_{r}=k}\,\prod_{i=1}^{r}\frac{q^{(s_{i}-1)k_{i}}}{[k_{i}]_{q}^{s_{i}}}=(-1)^{r}\underset{j_{i}<j_{i+1},\,i\in
I}{\sum_{1\leq j_{1}\leq j_{2}\leq\ldots\leq j_{w}\leq
n}}\,\prod_{i=1}^{w}\frac{q^{j_{i}}}{[j_{i}]_{q}},$

where $w=w({\bf s})=\sum_{i=1}^{r}s_{i}$,
$I=\\{s_{1},s_{1}+s_{2},\ldots,s_{1}+s_{2}+\dots+s_{r-1}\\}$, and the sum on
the right is taken over all integers $j_{1},\ldots,j_{w}$ satisfying the
conditions $1\leq j_{i}\leq n$, $j_{i}<j_{i+1}$ for $i\in I$, and $j_{i}\leq
j_{i+1}$ otherwise.

From Theorem A we get a kind of duality for finite multiple $q$-harmonic sums
$z_{n}$ at roots of unity.

######  Theorem 3.1.

Let $n,s_{1},\ldots,s_{r}$ be positive integers. Then

$z_{n}({\bf s};\zeta_{n})=(-1)^{r}\underset{j_{i}<j_{i+1},\,i\in
I}{\sum_{1\leq j_{1}\leq j_{2}\leq\ldots\leq
j_{w}<n}}\,\prod_{i=1}^{w}\frac{\zeta_{n}^{j_{i}}}{[j_{i}]_{\zeta_{n}}},$ (7)

where $w=w({\bf s})=\sum_{i=1}^{r}s_{i}$,
$I=\\{s_{1},s_{1}+s_{2},\ldots,s_{1}+s_{2}+\dots+s_{r-1}\\}$, and the sum on
the right is taken over all integers $j_{1},\ldots,j_{w}$ satisfying the
conditions $1\leq j_{i}\leq n$, $j_{i}<j_{i+1}$ for $i\in I$, and $j_{i}\leq
j_{i+1}$ otherwise.

###### Proof.

To get (7), we replace $n$ by $n-1$, and $q$ by a primitive root of unity
$\zeta_{n}$ in Theorem A, and apply Lemma 3.1. ∎

Proof of Theorem 1.2. Let ${\bf s}=(\\{1\\}^{a},2,\\{1\\}^{b})$ and $w=w({\bf
s})=a+b+2$. Applying Theorem 3.1 and noticing that
$I=\\{1,2,\ldots,a,a+2,a+3,\ldots,a+b+1\\}$, we get

$\begin{split}-z_{n}(\\{1\\}^{a},\,&2,\\{1\\}^{b};\zeta_{n})=(-1)^{a+b}\sum_{1\leq
j_{1}<j_{2}<\ldots<j_{a+1}\leq
j_{a+2}<j_{a+3}<\ldots<j_{w}<n}\,\prod_{i=1}^{w}\frac{\zeta_{n}^{j_{i}}}{[j_{i}]_{\zeta_{n}}}\\\
&=(-1)^{a+b}\sum_{0<n-j_{1}<\ldots<n-j_{a+1}\leq
n-j_{a+2}<n-j_{a+3}<\ldots<n-j_{w}<n}\,\prod_{i=1}^{w}\frac{\zeta_{n}^{n-j_{i}}}{[n-j_{i}]_{\zeta_{n}}}.\end{split}$

Applying identity (6), we get

$\begin{split}-z_{n}(&\\{1\\}^{a},2,\\{1\\}^{b};\zeta_{n})=\sum_{n>j_{1}>\ldots>j_{a+1}\geq
j_{a+2}>j_{a+3}>\ldots>j_{w}\geq
1}\frac{1}{[j_{1}]_{\zeta_{n}}\cdots[j_{w}]_{\zeta_{n}}}\\\
&=\sum_{n>j_{1}>\ldots>j_{a+1}>j_{a+3}>\ldots>j_{w}\geq
1}\frac{1}{[j_{1}]_{\zeta_{n}}\cdots[j_{a}]_{\zeta_{n}}[j_{a+1}]_{\zeta_{n}}^{2}[j_{a+3}]_{\zeta_{n}}\cdots[j_{w}]_{\zeta_{n}}}+z_{n}(\\{1\\}^{w};\zeta_{n}).\end{split}$

Noticing that

$\frac{1}{[j_{a+1}]_{\zeta_{n}}^{2}}=\frac{(1-\zeta_{n})[j_{a+1}]_{\zeta_{n}}+\zeta_{n}^{j_{a+1}}}{[j_{a+1}]_{\zeta_{n}}^{2}}=\frac{1-\zeta_{n}}{[j_{a+1}]_{\zeta_{n}}}+\frac{\zeta_{n}^{j_{a+1}}}{[j_{a+1}]_{\zeta_{n}}^{2}},$

we obtain

$-z_{n}(\\{1\\}^{a},2,\\{1\\}^{b};\zeta_{n})=(1-\zeta_{n})z_{n}(\\{1\\}^{w-1};\zeta_{n})+z_{n}(\\{1\\}^{b},2,\\{1\\}^{a};\zeta_{n})+z_{n}(\\{1\\}^{w};\zeta_{n}).$

Therefore, by (1),

$\begin{split}z_{n}(\\{1\\}^{a},2,\\{1\\}^{b};\zeta_{n})&+z_{n}(\\{1\\}^{b},2,\\{1\\}^{a};\zeta_{n})=-(1-\zeta_{n})z_{n}(\\{1\\}^{w-1};\zeta_{n})-z_{n}(\\{1\\}^{w};\zeta_{n})\\\
&=-\frac{1}{n}(1-\zeta_{n})\binom{n}{w}(1-\zeta_{n})^{w-1}-\frac{1}{n}\binom{n}{w+1}(1-\zeta_{n})^{w}\\\
&=-\frac{1}{n}\binom{n+1}{w+1}(1-\zeta_{n})^{w}.\end{split}$

∎

## References

  * [1] H. Bachmann, Y. Takeyama, K. Tasaka, Cyclotomic analogues of finite multiple zeta values, Compositio Math. 154 (2018), 2701-2721; arXiv:1707.05008v2 [math.NT].
  * [2] H. Bachmann, Y. Takeyama, K. Tasaka, Special values of finite multiple harmonic q-series at roots of unity, (2020), IRMA Lectures in Mathematics and Theoretical Physics 31,“Algebraic Combinatorics, Resurgence, Moulds and Applications (CARMA)”, Vol. 2, 1-18; arXiv:1807.00411 [math.NT].
  * [3] Kh. Hessami Pilehrood, T. Hessami Pilehrood, R. Tauraso, Some $q$-congruences for homogeneous and quasi-homogeneous multiple $q$-harmonic sums, Ramanujan J. 43 (2017), 113-139.
  * [4] J. Zhao, _On $q$-analog of Wolstenholme type congruences for multiple harmonic sums_, Integers 13 (2013), A23.


# Instantaneous Clear Sky Radiative Forcings of Halogenated Gases

W. A. van Wijngaarden  Department of Physics and Astronomy, York University,
Canada<EMAIL_ADDRESS>W. Happer  Department of Physics, Princeton University,
USA<EMAIL_ADDRESS>
The clear sky instantaneous radiative forcings of the 14 halogenated gases
previously shown to have the largest contribution to global warming, were
found. The calculation used the absorption cross sections for the halogenated
gases which are assumed to be independent of temperature as well as over 1/3
million line strengths for the five main naturally occurring greenhouse gases:
H2O, CO2, O3, CH4 and N2O, from the Hitran database. The total radiative
forcing of the halogenated gases at their 2020 concentrations is 0.52 (0.67)
W/m2 at the tropopause (mesopause). Over half of this forcing is due to CFC11
and CFC12 whose concentrations are steadily declining as a result of the
Montreal Protocol. The rate of total forcing change for all 14 halogenated
gases is 1.5 (2.2) mW/m2/year at the tropopause (mesopause). The calculations
assumed a constant altitude concentration for all halogenated gases except
CFC11, CFC12 and SF6. Using the observed altitude dependence for those 3
molecules reduced their radiative forcings by about 10%. The global warming
potential (GWP) values were found to be comparable to those given by the
Intergovernmental Panel on Climate Change. The contribution of a gas to global
warming can be estimated using the forcing power per molecule defined as the
derivative of its radiative forcing with respect to its column density. For
the present atmosphere, the per-molecule forcing powers of halogenated gases
are orders of magnitude larger than those for the 5 naturally occuring
greenhouse gases because the latter have much higher concentrations and are
strongly saturated. However, the rates of concentration increase of the 5 main
greenhouse gases are orders of magnitude greater than that of any halogenated
gas. Assuming the temperature increase caused by each gas is proportional to
its radiative forcing increase, the 14 halogenated gases are responsible for
only 2% of the total global warming.

##  1 Introduction

Halogenated molecules frequently contain one or more carbon atoms as well as
various halogen atoms and hydrogen. Chlorofluorocarbons (CFCs) are nontoxic
and nonflammable. CFC11 and CFC12, commercially known as freon 11 and freon 12
respectively, were discovered about a century ago and widely used in
refrigeration. When released into the atmosphere, these molecules can travel
to the stratosphere where they break down releasing chlorine atoms. In the
1970s, it was found that chlorine catalyzed the destruction of the
stratospheric ozone layer [1, 2]. An international agreement known as the
Montreal Protocol was created in 1987 to eliminate the use of CFCs [3].
Subsequent amendments seek to phase out hydrochlorofluorocarbons (HCFCs) and
hydrofluorocarbons (HFCs). These molecules were developed to replace CFCs.
They have much less ozone destroying ability but along with the CFCs are
greenhouse gases as pointed out by the Integovernmental Panel on Climate
Change (IPCC) [4].

There have been a number of calculations of the radiative forcings and global
warming potentials for various halogenated molecules [5, 6, 7]. These
calculations used various radiative transfer models along with results from
the Hitran database for the absorption cross sections for the halogenated
gases as well as line strengths for greenhouse gases such as carbon dioxide
[8]. The calculations use temperature profiles and altitudinal profiles of
greenhouse gases such as water vapor at various locations on the Earth.
Unfortunately, these profiles are not always readily available which would
facilitate comparison of different calculations.

The instantaneous radiative forcing was calculated for 14 molecules found to
be the strongest greenhouse gases of the halogenated molecules. The
calculations were done for a clear sky and use the standard US temperature
profile for the midlatitude [9]. This expands our previous work that found the
radiative forcings of the 5 main naturally occurring greenhouse gases H2O,
CO2, O3, CH4 and N2O [10]. Our calculations use the observed altitudinal
profiles of these gases [11]. The radiative forcings of the halogenated gases
using their 2020 concentrations, were found at the tropopause and top of the
atmosphere or mesopause. The rate of change of the total halogenated gas
forcing was found using the observed concentration change rates of the various
gases. The warming effectiveness of each molecule was found by taking the
derivative of its radiative forcing with respect to its column density to
obtain the per-molecule forcing power. Finally, the global warming potentials
were obtained and compared to values in the literature.

##  2 Concentration of Greenhouse Gases

Careful measurements of surface greenhouse gas concentrations have been
ongoing for more than half a century. The pioneering work of C. Keeling
studying carbon dioxide concentrations at Mauna Loa has been extended to
include a number of halogenated species as shown in Fig. 1 [12]. The graph
clearly shows the effectiveness of the Montreal Protocol as the concentrations
of the two most abundant halogenated gases, CFC11 and CFC12 have decreased by
about 20% and 10% respectively, since their peak concentrations in the 1990s.
All of the halogenated molecules have an entirely anthropogenic origin with
the exception of CF4 whose preindustrial level is estimated to be 34 ppt [12].

Figure 1: Measurements of Halogenated Gas Concentrations. The data for CFC114
are taken from Vollmer et al [13] while those for HFC23 are from Stanley et al
[14]. The data for all other gases were obtained by NOAA’s Global Greenhouse
Gas Reference Network at Mauna Loa, Hawaii [12]. Note the concentrations of
several gases have been multiplied by 10 in order to be visible on the graph
as indicated in the legend.

Radiation transfer in the cloud-free atmosphere of the Earth is controlled by
the temperature $T=T(z)$ at the altitude $z$ and the number densities,
$N^{\\{i\\}}=N^{\\{i\\}}(z)$ of the $i$th type of molecule. Representative
midlatitude altitude profiles of temperature [9], and greenhouse gas
concentrations [11], are shown in Fig. 2. Altitude profiles directly measured
by radiosondes in ascending balloons are always more complicated than those of
Fig. 2, which can be thought of as appropriate average profiles for the year
2020.

The standard US temperature profile for the midlatitude, plotted as the solid
blue line in the left panel of Fig. 2, consists of 5 linear segments. The
average surface temperature of the Earth is 288.7 K. This decreases to 217.2 K
at the tropopause altitude, $z_{\rm tp}$, of 11 km and remains constant until
20 km. The temperature in the stratosphere increases to 229.2 K at 32 km with
an additional increase to 271.2 K at 47 km as a result of absorption of
incoming solar UV radiation by ozone molecules. Finally, the temperature
decreases in the mesosphere to 187.5 K at the mesopause altitude, $z_{\rm
mp}$, of 86 km.

Figure 2: Left. A midlatitude atmospheric temperature profile, $T=T(z)$. The
Earth’s mean surface temperature $T(0)=288.7$ K. Right. Standard observed
concentrations, $C^{\\{i\\}}_{\rm sd}$ for greenhouse molecules versus
altitude $z$.

The concentrations for the ith greenhouse gas, $C^{\\{i\\}}_{\rm sd}$, based
on observations [11], are shown as functions of altitude on the right of Fig.
2 for the 5 main greenhouse gases. The 2020 surface concentrations are $7,750$
ppm of H2O, $1.9$ ppm of CH4, $0.33$ ppm of N2O and 10 ppt for SF6 [12]. The
O3 concentration peaks at $7.8$ ppm at an altitude of 35 km, while the
concentrations of CO2 and CF4 were 413 ppm and 86 ppt, respectively, at all
altitudes.

The dependence of the concentrations relative to the surface of CFC11 and
CFC12 with altitude is shown in Fig. 3. The concentrations of other
halogenated gases undoubtedly also depend on altitude but unfortunately few
observations exist. This study considered a constant altitude dependence of
the halogenated gas concentrations other than for CFC11, CFC12 and SF6.

Integrating the number density, $N^{\\{i\\}}$, over an atmospheric column
having a cross sectional area of 1 cm2 yields the column number density of the
$i$th type of molecule $\hat{N}^{\\{i\\}}_{\rm sd}$ listed in Table 2.

Figure 3: Observed concentrations of CFC11 and CFC12 relative to the surface
were found using data taken by the ACE Fourier Transform Spectrometer in 2009
between latitudes of 30o and 45o N [15].

##  3 Calculation of Radiative Forcing

Our previous publications describe in detail how to calculate radiation
transport from the Earth’s surface through the atmosphere [10, 16] and is here
only briefly discussed. Radiation transport is governed by the Schwarzschild
equation [17] in cloud-free air where scattering is negligible.

$\cos\theta\frac{\partial\tilde{I}}{\partial\tau}=\tilde{B}-\tilde{I}$ (1)

Here ${\tilde{I}}={\tilde{I}}(\nu,z,\theta)$ is the spectral intensity of a
pencil of radiation of frequency between $\nu$ and $\nu+d\nu$ at altitude $z$.
The pencil makes an angle $\theta$ to the vertical. The Planck intensity is
given by

$\tilde{B}(\nu,T)=\frac{2hc^{2}\nu^{3}}{e^{\nu c\,h/(k_{\rm B}T)}-1}$ (2)

The radiation frequency, $\nu$ has units of cm-1, $h$ is Planck’s constant,
$k_{B}$ is Boltzmann’s constant and $c$ is the speed of light.

The optical depth is defined by

$\tau(z,\nu)=\int_{0}^{z}dz^{\prime}\kappa(z^{\prime},\nu),$ (3)

where the net attenuation coefficient due to molecules absorbing and
reemitting light of frequency $\nu$ at altitude $z$ is given by

$\kappa(z,\nu)=\sum_{i}N^{\\{i\\}}(z)\sigma^{\\{i\\}}(z,\nu).$ (4)

Molecule | Temperature | Pressure | $\nu_{min}$ | $\nu_{max}$ | Cross Section
---|---|---|---|---|---
| (K) | (mbar) | (cm-1) | (cm-1) | ($10^{-17}$ cm2)
Chlorofluorocarbons |  |  |  |  |
CFC10 (CCl4) | 295.7 | 760.0 | 700 | 860 | 6.73
CFC11 | 293.0 | 760.0 | 710 | 1290 | 9.95
| 273.0 | 7.5 | 710 | 1290 | 9.96
| 232.6 | 7.5 | 710 | 1290 | 9.96
| 216.7 | 49.9 | 710 | 1290 | 9.96
| 216.6 | 201 | 710 | 1290 | 9.96
| 191.7 | 7.5 | 710 | 1290 | 9.97
CFC12 | 293.7 | 761.0 | 800 | 1270 | 13.5
CFC113 | 298.1 | 760.0 | 620 | 5000 | 14.6
CFC114 | 298.1 | 760.0 | 600 | 5000 | 17.4
Hydrochlorofluorocarbons |  |  |  |  |
HCFC22 | 294.8 | 761.6 | 730 | 1380 | 10.5
HCFC141b | 298.1 | 760.0 | 550 | 6500 | 8.57
HCFC142b | 298.1 | 760.0 | 600 | 6500 | 11.1
Hydrofluorocarbons |  |  |  |  |
HFC23 | 293.7 | 762.4 | 950 | 1500 | 12.3
HFC125 | 298.1 | 760.0 | 510 | 6500 | 17.5
HFC134a | 295.8 | 760.8 | 750 | 1600 | 13.2
HFC143a | 298.1 | 760.0 | 500 | 6500 | 14.1
Fully Fluorinated Species |  |  |  |  |
SF6 | 189.2 | 50.4 | 780 | 1100 | 19.9
| 298.1 | 760.0 | 560 | 6500 | 21.2
| 323.1 | 760.0 | 560 | 6500 | 21.6
CF4 | 295.0 | 760.5 | 1190 | 1336 | 19.3
Table 1: Hitran Absorption Cross Section Datasets Used in This Work.

Here $N^{\\{i\\}}(z)$ is the density of a greenhouse gas molecule of type $i$
and $\sigma^{\\{i\\}}=\sigma^{\\{i\\}}(z,\nu)$ is its absorption cross section
for radiation of frequency $\nu$ at the altitude $z$. The cross section can
depend strongly on altitude because temperature and pressure are functions of
altitude. Temperature controls the distribution of the molecules between
translational, rotational and vibrational states. Pressure, together with
temperature, determines the width of the molecular resonance lines. The
halogenated gas absorption cross sections were obtained from the Hitran
website [8] and are listed in Table 1. Halogens are very electronegative which
implies that vibrations of halogenated molecules can have unusually large
transition electric dipole moments, and therefore large absorption cross
sections for thermal radiation. The uncertainty of the Hitran cross sections
is given to be about 8% [7]. Integrating over frequency gave the total cross
sections listed in the final column of Table 1. The total cross sections are
believed to be independent of temperature [18]. Table 1 shows the Hitran
absorption cross sections for CFC11 and SF6 vary less than 0.1% and 2%
respectively, over a considerable range of temperatures and pressures. This
will be discussed more in Section 6.

The cross section can be written as the sum of partial cross sections
$\sigma_{ul}$, corresponding to each Bohr transition frequency $\nu_{ul}$,

$\sigma^{\\{i\\}}=\sum_{ul}\sigma_{ul}.$ (5)

The partial cross section, $\sigma_{ul}$, is the product of a lineshape
function, $G_{ul}=G_{ul}(\nu,\tau)$, and a line strength, $S_{ul}=S_{ul}(T)$,

$\sigma_{ul}=G_{ul}S_{ul}.$ (6)

The absorption of the 5 main greenhouse gases, H2O, CO2, O3, CH4 and N2O were
obtained using the Hitran line strengths [10].

The optical depth from the surface to the top of the radiative atmosphere, the
altitude $z_{\rm mp}$ of the mesopause, is

$\tau_{\infty}=\tau_{\rm mp}=\int_{0}^{z_{\rm
mp}}dz^{\prime}\kappa(z^{\prime},\nu).$ (7)

As indicated by the notation (7), we have assumed that the optical depth
$\tau_{\rm mp}$ at the mesopause altitude $z_{\rm mp}$ differs negligibly from
the optical depth $\tau_{\infty}$ at infinite altitude since there is so
little opacity of the atmosphere above the mesopause.

The Schwarzschild equation (1) can be solved to find the intensity [10]

$\displaystyle\hbox{For $\varsigma>0$}:\quad\tilde{I}(\tau,\varsigma)$
$\displaystyle=$
$\displaystyle+\varsigma\int_{0}^{\tau}d\tau^{\prime}e^{-\varsigma(\tau-\tau^{\prime})}\tilde{B}(\tau^{\prime})+e^{-\varsigma\tau}\tilde{I}(0,\varsigma)$
(8) $\displaystyle\hbox{For $\varsigma<0$}:\quad\tilde{I}(\tau,\varsigma)$
$\displaystyle=$
$\displaystyle-\varsigma\int_{\tau}^{\tau_{\infty}}d\tau^{\prime}e^{-\varsigma(\tau-\tau^{\prime})}\tilde{B}(\tau^{\prime})$
(9)

where $\varsigma=\sec\theta$.

The upwards flux defined by

$\tilde{Z}=\int_{4\pi}d\Omega\,\cos\theta\,\tilde{I}.$ (10)

can be rewritten after substituting (8) and (9) into (10) to give

$\displaystyle{{\tilde{Z}}\over{2\pi}}$ $\displaystyle=$
$\displaystyle\int_{0}^{\tau}d\tau^{\prime}E_{2}(\tau-\tau^{\prime})\tilde{B}(\tau^{\prime})-\int_{\tau}^{\tau_{\infty}}d\tau^{\prime}E_{2}(\tau^{\prime}-\tau)\tilde{B}(\tau^{\prime})+\epsilon_{s}\tilde{B}_{s}E_{3}(\tau)$
(11) $\displaystyle=$
$\displaystyle-\int_{0}^{\tau_{\infty}}d\tau^{\prime}E_{3}(|\tau-\tau^{\prime}|){{\partial\tilde{B}(\tau^{\prime})}\over{\partial\tau^{\prime}}}+\tilde{B}(\tau_{\infty})E_{3}(\tau_{\infty}-\tau).$

Here, $\tilde{B}_{s}=\tilde{B}(T_{s})$, is the Planck intensity evaluated at
the surface temperature $T_{s}$. The exponential-integral functions,
$E_{n}(\tau)$, that account for slant paths of radiation between different
altitudes are defined for integers $n=1,2,3,\ldots$ by

$E_{n}(\tau)=\int_{1}^{\infty}d\varsigma\,\varsigma^{-n}\,e^{-\varsigma\tau}.$
(12)

Figure 4: Effects of different concentrations of CFC12 on the spectral flux at
the mesopause altitude, $z_{\rm mp}$. The smooth blue line is the spectral
flux from a surface at the temperature of $T_{s}$ = 288.7 K without any
greenhouse gases present. Both the black and red lines show the flux for an
atmosphere containing the 5 main greenhouse gases at their concentrations of
Fig. 2 as well as CFC12. The effect of CFC12 is only apparent if its
concentration is 10 times its observed 2020 value, 499 ppt, as shown by the
red curve.

The spectral forcing, $\tilde{F}$, is defined as the difference between the
spectral flux $\pi\tilde{B}_{s}$ through a transparent atmosphere from a
surface with temperature $T_{s}$, and the spectral flux $\tilde{Z}$ for an
atmosphere having greenhouse gases,

$\tilde{F}=\pi\tilde{B}_{s}-\tilde{Z}\\\ .$ (13)

The spectral flux $\tilde{Z}$ at the mesopause altitude is shown in Fig. 4 for
the case of the 5 main greenhouse gases as well as CFC12.

Figure 5: Altitude dependence of frequency integrated flux $Z$ of (14). The
flux for three concentrations of CO2 are shown, the standard concentration of
Fig. 2, twice and half that value. The other greenhouse gases have the
standard concentrations of Fig. 2. The vertical dashed line is the flux
$\sigma_{\rm SB}T_{s}^{4}=394$ W m-2 for a transparent atmosphere with a
surface temperature $T_{s}=288.7$ K.

The frequency integrals of the flux (10) and the forcing (13) are

$\displaystyle Z$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\nu\,\tilde{Z},$ (14) $\displaystyle F$
$\displaystyle=$ $\displaystyle\int_{0}^{\infty}d\nu\,\tilde{F}=\sigma_{\rm
SB}T_{s}^{4}-Z,$ (15)

where $\sigma_{SB}$ is the Stefan Boltzmann constant. Fig. 5 shows the
frequency integrated flux versus altitude for three different concentrations
of CO2 in addition to the four other main greenhouse gases. A similar plot to
Fig. 5 where one would show the result of increasing the concentration of any
of the halogenated gases from its preindustrial level to its 2020
concentration results in an imperceptible change from the green curve.

Molecule |  $C^{\\{i\\}}$(2020) | $dC^{\\{i\\}}/dt$ | Altitude | Column Density $\hat{N}^{\\{i\\}}_{\rm sd}$ | Forcing (mW/m2)
---|---|---|---|---|---
$i$ | (ppt) | (ppt/yr) | Dependence | ($10^{14}$ mol/cm2) | $z_{\rm tp}$ | $z_{\rm mp}$
Chlorofluorocarbons |  |  |  |  |  |
CFC10 (CCl4) | 78 | -1.1 | No | 16.7 | 20.6 | 25.7
CFC11 | 226 | -1.9 | Yes | 44.5 | 84.6 | 104
|  |  | No | 48.2 | 91.3 | 118
CFC12 | 499 | -3.8 | Yes | 101 | 219 | 280
|  |  | No | 107 | 227 | 302
CFC113 | 69 | -0.64 | No | 14.7 | 32.3 | 40.7
CFC114 | 16 | -0.02 | No | 3.42 | 7.64 | 9.57
Hydrochlorofluorocarbons |  |  |  |  |  |
HCFC22 | 255 | 2.2 | No | 54.4 | 77.6 | 102
HCFC141b | 25 | 0.10 | No | 5.34 | 5.90 | 7.40
HCFC142b | 22 | -0.14 | No | 4.70 | 6.29 | 8.39
Hydrofluorocarbons |  |  |  |  |  |
HFC23 | 38 | 1.1 | No | 8.11 | 10.4 | 14.0
HFC125 | 32 | 3.1 | No | 6.83 | 10.9 | 14.6
HFC134a | 114 | 6.1 | No | 2.43 | 27.5 | 36.7
HFC143a | 25 | 1.5 | No | 5.34 | 6.07 | 8.14
Fully Fluorinated Species |  |  |  |  |  |
SF6 | 10 | 0.37 | Yes | 1.97 | 7.33 | 9.18
|  |  | No | 2.13 | 7.90 | 10.5
CF4 | 86 (34) | 0.78 | No | 18.4 | 6.55 | 8.75
|  |  |  | Halogenated Gas Total | 521 | 667
Main Greenhouse Gases |  |  |  |  |  |
CO2 | 4.13 (2.80) $\times 10^{8}$ | $2.5\times 10^{6}$ |  | $8.54\times 10^{7}$ | 2950 | 1670
CH4 | 1.9 (0.8) $\times 10^{6}$ | $8\times 10^{3}$ | Yes | $3.73\times 10^{5}$ | 834 | 761
N2O | 3.3 (2.7) $\times 10^{5}$ | $8\times 10^{2}$ | Yes | $6.56\times 10^{4}$ | 277 | 259
Table 2: Instantaneous Clear Sky Radiative Forcings at the altitude $z_{\rm
tp}=11$ km of the tropopause and $z_{\rm mp}=86$ km of the mesopause. The
total anthropogenic forcing was found relative to the preindustrial levels
which was zero unless listed in brackets. The observed rates of the
concentration changes of the halogenated gases is given in Fig. 1 while the
observed concentration changes for CH4 and N2O are discussed further in
reference [20]. The values for CFC11, CFC12 and SF6 lacking altitude
dependence were not considered when finding the halogenated gas total
radiative forcings.

The radiative forcings of the various gases resulting from the increase in the
respective concentration relative to their preindustrial levels are given in
Table 2 at the tropopause and mesopause altitudes. Table 2 as well as
subsequent tables do not list water vapor, the most significant greenhouse
gas. As shown in Table 2 of reference [10], the total radiative forcing of
water vapor is about twice that of CO2. Clouds further increase the importance
of H2O for radiative transfer in Earth’s atmosphere. All halogenated gases
except for CF4 were not present in the atmosphere prior to the industrial
revolution. The effect of the altitude dependence on forcing was examined for
CFC11, CFC12 and SF6. Taking a constant concentration at all altitudes
increased their forcings by as much as 8% at the tropopause and 14% at the
mesopause. It is reasonable to assume that the radiative forcings of the other
halogenated species are similarly affected.

The effect due to all 14 halogenated species present in the atmosphere gives a
forcing of 521 (667) mW/m2 at the tropopause (mesopause) which is $\leq$ 0.2%
of the total flux as shown in Fig. 5. Nearly 60% of these total forcings are
due to CFC11 and CFC12 that have the highest halogenated gas concentrations.
The rate of change of each halogenated species was used to find the rate of
change of the total halogenated gas radiative forcing to be 1.5 (2.2) mW/m2/yr
at the tropopause (mesopause). This means the total halogenated gas radiative
forcing will remain significantly less than that due to carbon dioxide and
methane for the foreseeable future.

##  4 Concentration Dependence of Forcing

The frequency integrated forcing, $F$, of (15) depends on the altitude $z$ and
on the column densities of the five main greenhouse gases as well as the 14
halogenated species given in Table 2.

$F=F(z,\hat{N}^{\\{1\\}},\ldots,\hat{N}^{\\{19\\}}).$ (17)

We assume the temperature $T$ and densities $N^{\\{i\\}}$ have the same
altitude profiles as shown in Fig. 2. An important special case of (17) is the
forcing, $F_{\rm sd}$, when each greenhouse gas $i$ is present at its 2020
column density $\hat{N}^{\\{i\\}}_{\rm sd}$ of Table 2,

$F_{\rm sd}(z)=F(z,\hat{N}^{\\{1\\}}_{\rm sd},\ldots,\hat{N}^{\\{19\\}}_{\rm
sd}).$ (18)

A second special case of (17) is the hypothetical, per molecule standard
forcing, $F^{\\{i\\}}_{\rm sd}$, when the atmosphere contains only molecules
of type $i$ at their standard column density,
$\hat{N}^{\\{i\\}}=\hat{N}^{\\{i\\}}_{\rm sd}$, and the concentrations of the
other greenhouse vanish, $\hat{N}^{\\{j\\}}=0$ if $j\neq i$,

$F^{\\{i\\}}_{\rm sd}(z)=F(z,0,\ldots,0,\hat{N}^{\\{i\\}}_{\rm
sd},0,\ldots,0).$ (19)

We define the forcing power per added molecule as

$P^{\\{i\\}}(z,\hat{N}^{\\{1\\}},\ldots,\hat{N}^{\\{n\\}})=\frac{\partial
F}{\partial\hat{N}^{\\{i\\}}}.$ (20)

The densities of greenhouse gases $j$ with $j\neq i$ are held constant in the
partial derivative of (20). If the units of $F$ are taken to be W m-2 and the
units of $\hat{N}^{\\{i\\}}$ are taken to be molecules m-2, then the units of
$P^{\\{i\\}}$ will be W molecule-1.

We define a finite forcing increment for the $i$th type of greenhouse molecule
as

$\Delta F^{\\{i\\}}(z,f)=F(z,\hat{N}^{\\{1\\}}_{\rm
sd},\ldots,\hat{N}^{\\{i-1\\}}_{\rm sd},f\hat{N}^{\\{i\\}}_{\rm
sd},\hat{N}^{\\{i+1\\}}_{\rm sd},\ldots,\hat{N}^{\\{n\\}}_{\rm sd})-F_{\rm
sd}.$ (21)

Molecule | $P^{\\{i\\}}_{\rm ot}(z)$ | $P^{\\{i\\}}_{\rm sd}(z,0)$ | $P^{\\{i\\}}_{\rm sd}(z,1)$
---|---|---|---
$i$ | $z_{\rm tp}$ | $z_{\rm mp}$ | $z_{\rm tp}$ | $z_{\rm mp}$ | $z_{\rm tp}$ | $z_{\rm mp}$
Chlorofluorocarbons |  |  |  |  |  |
CFC10 (CCl4) | 1.62 | 2.10 | 1.25 | 1.55 | 1.24 | 1.54
CFC11 | 2.09 | 2.64 | 2.03 | 2.41 | 1.87 | 2.31
CFC12 | 2.36 | 3.12 | 2.29 | 2.83 | 2.12 | 2.71
CFC113 | 2.58 | 3.50 | 2.22 | 2.82 | 2.15 | 2.74
CFC114 | 2.77 | 3.81 | 2.31 | 2.84 | 2.25 | 2.84
Hydrochlorofluorocarbons |  |  |  |  |  |
HCFC22 | 1.69 | 2.32 | 1.42 | 1.88 | 1.42 | 1.86
HCFC141b | 1.52 | 2.02 | 1.11 | 1.82 | 1.10 | 1.38
HCFC142b | 1.69 | 2.33 | 1.39 | 1.81 | 1.33 | 1.78
Hydrofluorocarbons |  |  |  |  |  |
HFC23 | 1.60 | 2.25 | 1.29 | 1.74 | 1.28 | 1.73
HFC125 | 2.25 | 3.13 | 1.62 | 2.14 | 1.60 | 2.13
HFC134a | 1.66 | 2.34 | 1.13 | 1.51 | 1.13 | 1.51
HFC143a | 1.70 | 2.37 | 1.14 | 1.53 | 1.14 | 1.52
Fully Fluorinated Species |  |  |  |  |  |
SF6 | 4.28 | 5.68 | 3.75 | 4.68 | 3.71 | 4.67
CF4 | 2.00 | 2.85 | 0.63 | 0.81 | 0.62 | 0.80
Main Greenhouse Gases |  |  |  |  |  |
CO2 [16] | 0.276 | 0.349 | 0.22 | 0.25 | $9.0\times 10^{-5}$ | $4.9\times 10^{-5}$
CH4 [16] | 0.052 | 0.072 | 0.021 | 0.027 | $2.8\times 10^{-3}$ | $2.6\times 10^{-3}$
N2O [16] | 0.170 | 0.227 | 0.073 | 0.091 | $2.1\times 10^{-2}$ | $2.0\times 10^{-2}$
Table 3: Forcing powers (20) per additional molecule in units of $10^{-21}$ W
at the altitude $z_{\rm tp}=11$ km of the tropopause and $z_{\rm mp}=86$ km of
the mesopause. Altitude dependences were taken into account as given in Fig. 2
and Fig. 3 while all other species had altitude independent concentrations.
The surface temperature was $T_{s}=288.7$ K, and the altitude profile of
temperature is given in Fig. 2. $P^{\\{i\\}}_{\rm ot}(z)$ of (24) is for the
optically-thin limit. $P^{\\{i\\}}_{\rm sd}(z,0)$ from (22) is for an
atmosphere that previously had no molecules of type $i$ (so
$\hat{N}^{\\{i\\}}=0$) but all other greenhouse molecules had standard
concentrations. $P^{\\{i\\}}_{\rm sd}(z,1)$ from (22) is for a single molecule
of type $i$ added to an atmosphere that previously had standard densities for
all 5 main greenhouse gases.

Differentiating (21) with respect to $f$ we find

$\frac{\partial\Delta F^{\\{i\\}}}{\partial f}(z,f)=\hat{N}^{\\{i\\}}_{\rm
sd}P^{\\{i\\}}_{\rm sd}(z,f),$ (22)

where $P^{\\{i\\}}_{\rm sd}(z,f)$ is the forcing power per additional molecule
of type $i$ when these molecules have the column density
$\hat{N}^{\\{i\\}}=f\hat{N}^{\\{i\\}}_{\rm sd}$ and all other types of
greenhouse molecules have their standard column densities.

We now consider the optically thin limit, where the concentrations of
greenhouse gases are sufficiently low that the optical depths $\tau$ of (3)
will be small, $\tau\ll 1$, for all frequencies $\nu$ and at all altitudes
$z$. The frequency integral of the spectral forcing (15) at altitude $z$ can
then be written as

$\displaystyle F_{\rm ot}(z)$ $\displaystyle=$
$\displaystyle\sum_{i}\hat{N}^{\\{i\\}}P^{\\{i\\}}_{\rm ot}(z)$ (23)

where the forcing power per molecule of type $i$ is

$P^{\\{i\\}}_{\rm
ot}(z)=\frac{1}{2}\int_{0}^{z}dz^{\prime}\frac{N^{\\{i\\}^{\prime}}}{\hat{N}^{\\{i\\}}}\left[\Pi^{\\{i\\}}(T^{\prime},T_{s})-\Pi^{\\{i\\}}(T^{\prime},T^{\prime})\right]+\frac{1}{2}\int_{z}^{\infty}dz^{\prime}\frac{N^{\\{i\\}^{\prime}}}{\hat{N}^{\\{i\\}}}\Pi^{\\{i\\}}(T^{\prime},T^{\prime}).$
(24)

Here $N^{\\{i\\}^{\prime}}=N^{\\{i\\}}(z^{\prime})$, and
$T^{\prime}=T(z^{\prime})$. The mean power absorbed by a greenhouse gas
molecule of temperature $T$ from thermal equilibrium radiation of temperature
$T^{\prime}$ is

$\Pi^{\\{i\\}}(T,T^{\prime})=4\pi\sum_{ul}S_{ul}^{\\{i\\}}(T)\tilde{B}(\nu_{ul},T^{\prime}).$
(25)

For $z^{\prime}>z$ we see from (24) that the $dz^{\prime}\,N^{\prime}$
molecules in the altitude interval $z^{\prime}$ to $z^{\prime}+dz^{\prime}$
each emit the power $\Pi(T^{\prime},T^{\prime})$, of which half goes to outer
space and half goes down through the reference plane at altitude $z$,
diminishing the net flux through the reference plane by
$\Pi(T^{\prime},T^{\prime})/2$. Molecules above the reference plane can only
cause positive forcing unlike molecules below the reference plane which can
cause either positive or negative forcing.

An important check of our work was to use (24) to calculate the forcing power
per molecule for the optically thin limit. The results displayed in column 2
of Table 3 agreed with that found by dividing the radiative forcing in the
limit of zero gas concentration, by the column density.

The forcing powers per molecule given in Table 3 are all on the order of
$10^{-21}$ W for the halogenated gases. The values are largest for the
optically thin case and slightly lower for $P^{\\{i\\}}_{\rm sd}(z,0)$ and
$P^{\\{i\\}}_{\rm sd}(z,1)$ because of saturation. The effect of satuation is
much larger for the 5 main greenhouse gas values as explained in reference
[16] because their number densities are up to $10^{7}$ times greater than
those of the halogenated gases.

The forcing power can be used to estimate the contribution of each gas to
global warming as follows. First, one assumes that the temperature increase
$\Delta T^{\\{i\\}}$ due to a change in the concentration of gas $i$ is
proportional to its radiative forcing change,

$\Delta T^{\\{i\\}}=c\ \Delta F^{\\{i\\}}$ (26)

Molecule | $\Delta F^{\\{i\\}}/{\Delta F_{\rm tot}}$
---|---
| (%)
Chlorofluorocarbons |
CFC10 (CCl4) | -0.51
CFC11 | -1.32
CFC12 | -2.99
CFC113 | -0.51
CFC114 | -0.02
Hydrochlorofluorocarbons |
HCFC22 | 1.16
HCFC141b | 0.04
HCFC142b | -0.07
Hydrofluorocarbons |
HFC23 | 0.52
HFC125 | 1.84
HFC134a | 2.56
HFC143a | 0.63
Fully Fluorinated Species |
SF6 | 0.51
CF4 | 0.18
Main Greenhouse Gases |
CO2 | 83.4
CH4 | 8.3
N2O | 6.2
Total | 100.0
Table 4: Relative Contribution of Greenhouse Gases to Total Radiative Forcing
Change. $\Delta F_{\rm{Tot}}$ represents the total radiative forcing change
caused by changes in the concentrations of all greenhouse gas molecules. The
relative forcing was found using the observed 2020 concentration change rate
of each gas $i$ as discussed in the text. A negative relative forcing change
occurs if the halogenated gas atmospheric concentration is decreasing as a
result of the Montreal Protocol.

where $c$ is a constant of proportionality whose estimation is difficult
because it includes feedback effects. The radiative forcing change at the
tropopause due to a concentration change occurring over time $\Delta t$ is
given by

$\Delta F^{\\{i\\}}={P^{\\{i\\}}_{\rm sd}(z_{\rm
tr},1)}{{dC^{\\{i\\}}}\over{dt}}\Delta t.$ (27)

Here, we have evaluated the forcing power per molecule at the tropopause which
is given in Table 3. The total radiative forcing change due to changes in all
$i$ greenhouse gases over time $\Delta t$ is given by

$\Delta F_{\rm Tot}=\sum_{i}\Delta F^{\\{i\\}}.$ (28)

Table 4 gives the relative contributions of the various gases to the total
radiative forcing change. This was found evaluating (27) and (28) using the
2020 observed concentration change rates given in Table 2. The halogenated gas
total is responsible for 2% of the total radiative forcing change.

The total temperature change due to changes in all $i$ greenhouse gas
concentrations over time $\Delta t$ is given by

$\Delta T_{\rm Tot}=\sum_{i}\Delta T^{\\{i\\}}.$ (29)

Equations (26), (28) and (29) imply

${{\Delta T^{\\{i\\}}}\over{\Delta T_{\rm Tot}}}={{\Delta
F^{\\{i\\}}}\over{\Delta F_{\rm Tot}}}.$ (30)

Surface temperatures have increased by about 1.1 K since 1850 according to the
IPCC [4]. More recent satellite measurements since 1979, show a slightly
higher warming rate of about 0.1 C per decade [19]. Since the observed
temperature increase is comparable to previous temperature increases when
there were negligible changes of greenhouse gas concentrations, there is no
way to know how much is due to currently increasing greenhouse gases. But if
one makes the extreme assumption that changes in greenhouse gas concentrations
have driven all of this observed warming, one finds the changes in the
halogenated gas concentrations are responsible for a warming of 0.002 C per
decade.

##  5 Global Warming Potential

A comparison of the greenhouse warming contribution of a gas $i$ relative to
carbon dioxide is provided by the global warming potential which is defined as
[20]

$\hbox{GWP}^{\\{i\\}}(T)=\frac{\langle\hbox{RF}^{\\{i\\}}\rangle}{\langle\hbox{RF}^{\\{{\rm
CO}_{2}\\}}\rangle}.$ (31)

The time-integrated radiative forcing at the tropopause per unit mass over an
observation time $T$ is

$\langle\hbox{RF}^{\\{i\\}}\rangle=\frac{P^{\\{i\\}}\theta^{\\{i\\}}(T)}{m^{\\{i\\}}}.$
(32)

The units of $\langle\hbox{RF}^{\\{i\\}}\rangle$ are joules per square meter
per unit mass. It is the amount of “forcing heat” acquired by a square-meter
column of the atmosphere during time $T$ after the pulse emission of one unit
mass of greenhouse gas of type $i$. The forcing time $\theta^{\\{i\\}}$ is

$\theta^{\\{i\\}}(T)=\int_{0}^{T}dt^{\prime}f^{\\{i\\}}(t^{\prime}).$ (33)

$f^{\\{i\\}}(t^{\prime})$ is the fraction of the excess greenhouse gas
molecules remaining in the atmosphere at a time $t^{\prime}$ after the “pulse
emission.” Excess greenhouse molecules can remain in the atmosphere long after
the greenhouse molecules of the original pulse have exchanged with the land
and oceans because they are replaced with equivalent molecules that continue
to provide radiative forcing. An example is the relatively short time needed
for the exchange of 14CO2 molecules from atmospheric nuclear weapon tests with
CO2 molecules of the ocean. It takes a much longer time for excess CO2
concentrations to decay away. The excess atmospheric concentration of
molecules, caused by the emitted pulse of the species $i$, can be removed from
the atmosphere by various mechanisms. For example, CO2 is absorbed by the
biosphere and oceans. CH4 is oxidized by OH radicals and other atmospheric
gases. N2O and the various halogenated gases are photodissociated by solar
ultraviolet radiation in the stratosphere, etc.

It is convenient to describe the excess fraction $f^{\\{i\\}}(t)=f(t)$ with a
sum of $n+1$ exponentially decaying components, of amplitudes $a_{j}$ and time
constants $\tau_{j}$ [21, 22]

$f(T)=\sum_{j=0}^{n}a_{j}e^{-T/\tau_{j}},\quad\hbox{where}\quad\sum_{j=0}^{n}a_{j}=1.$
(34)

Substituting (34) into (33) we find that the forcing time,
$\theta^{\\{i\\}}=\theta$, is

$\theta(T)=\sum_{j=0}^{n}a_{j}\tau_{j}\left(1-e^{-T/\tau_{j}}\right).$ (35)

For the limiting case of infinitely long time constants, $\tau_{j}\to\infty$,
one should make the replacement $\tau_{j}\left(1-e^{-T/\tau_{j}}\right)\to T$
in (35). From inspection of (35) one can conclude that the forcing time
$\theta^{\\{i\\}}$ is always less than or equal to the observation time,
$\theta^{\\{i\\}}(T)\leq T$. Using (35) in (32) we find that the global
warming potential is

$\hbox{GWP}^{\\{i\\}}(T)=\left(\frac{P^{\\{i\\}}}{P^{\\{{\rm
CO}_{2}\\}}}\right)\left(\frac{m^{\\{{\rm
CO}_{2}\\}}}{m^{\\{i\\}}}\right)\left(\frac{\theta^{\\{i\\}}(T)}{\theta^{\\{{\rm
CO}_{2}\\}}(T)}\right).$ (36)

Molecule | Mass | Lifetime[23] | This Work | IPCC 2021 [23]
---|---|---|---|---
$i$ |  $m^{\\{i\\}}$(amu) |  $\tau_{i}$(yr) | GWP(0) | GWP(20) | GWP(100) | GWP(20) | GWP(100)
Chlorofluorocarbons |  |  |  |  |  |  |
CFC10 (CCl4) | 151.1 | 32 | 4296 | 4179 | 2494 | 3810 | 2200
CFC11 | 137.4 | 52 | 6654 | 7225 | 5605 | 7430 | 5560
CFC12 | 120.9 | 102 | 8573 | 10182 | 10367 | 11400 | 11200
CFC113 | 187.4 | 93 | 5609 | 6601 | 6520 | 6860 | 6520
CFC114 | 170.9 | 89 | 6437 | 7540 | 7336 | 8260 | 9430
Hydrochlorofluorocarbons |  |  |  |  |  |  |
HCFC22 | 86.5 | 11.9 | 8026 | 5083 | 1812 | 5690 | 1960
HCFC141b | 116.9 | 9.4 | 4600 | 2491 | 820 | 2710 | 860
HCFC142b | 100.5 | 18 | 6470 | 5109 | 2201 | 5510 | 2300
Hydrofluorocarbons |  |  |  |  |  |  |
HFC23 | 70 | 228 | 8940 | 11195 | 13732 | 12400 | 14600
HFC125 | 120 | 30 | 6519 | 6223 | 3578 | 6740 | 3740
HFC134a | 102 | 14 | 5416 | 3770 | 1438 | 4140 | 1530
HFC143a | 84 | 51 | 6635 | 7179 | 5517 | 7840 | 5810
Fully Fluorinated Species |  |  |  |  |  |  |
SF6 | 146 | 3200 | 12423 | 16198 | 23207 | 18300 | 25200
CF4 | 88 | 50000 | 3444 | 4504 | 6529 | 5300 | 7380
Main Greenhouse Gases |  |  |  |  |  |  |
CO2 | 44 | See Caption | 1 | 1 | 1 | 1 | 1
CH4 [20] | 16 | 11.8 | 85.5 | 53.9 | 19.2 | 81.2 | 27.9
N2O [20] | 44 | 109 | 233 | 279 | 290 | 273 | 273
Table 5: Comparison of Global Warming Potentials calculated for this work
using equation (36) and the tropopause powers given in Table 3 with IPCC 2021
values. There are multiple pathways for CO2 to be removed from the atmosphere
as is discussed in the text.

Unlike the forcing powers $P^{\\{i\\}}$ per greenhouse molecule of type $i$,
which can be accurately calculated, the time dependence of the fractions
$f^{\\{i\\}}(T)$ is not well known. So, the global warming potentials GWP{i}
are of limited quantitative value and do not contribute in a clear way to
estimating the warming caused by a change in a greenhouse gas concentration.
To add further uncertainty, “indirect effects” are sometimes included in the
GWP{i}, for example, the effects of the CO2 and H2O molecules which result
from the oxidation of CH4. For this paper only the direct forcing of the
greenhouse gases was considered.

For CO2, Harvey [21] uses a five-exponential form of (34) with the parameters

$\left[\begin{array}[]{c}a_{0}\\\ a_{1}\\\ a_{2}\\\ a_{3}\\\
a_{4}\end{array}\right]=\left[\begin{array}[]{c}.131\\\ .201\\\ .321\\\
.249\\\
.098\end{array}\right]\quad\hbox{and}\quad\left[\begin{array}[]{c}\tau_{0}\\\
\tau_{1}\\\ \tau_{2}\\\ \tau_{3}\\\
\tau_{4}\end{array}\right]=\left[\begin{array}[]{c}\infty\\\ 362.9\\\ 73.6\\\
17.3\\\ 1.9\end{array}\right]\hbox{ yr}.$ (37)

Alternate parameterizations of $f^{\\{{{\rm CO}_{2}}\\}}$, for example, those
of Joos et al. [22], give GWPs that do not differ more than the uncertainties
of the estimates.

Table 5 shows the global warming potentials calculated with (36) using the
lifetimes of the various greenhouse gases given in the most recent IPCC report
[23] and equation (37). Our GWP values are consistent with those of the 2021
IPCC report. The somewhat larger values of GWP${}^{\\{\rm CH_{4}\\}}$ may be
due to indirect effects included in the IPCC calculations.

##  6 Summary

This work found the instantaneous clear sky radiative forcings of 14
halogenated molecules that have previously been found to exert the largest
greenhouse warming. CFC11 and CFC12 comprise about 60% of the total forcing
generated by the 14 halogenated molecules and their concentrations have been
decreasing in recent years at 2% and 4% per year, respectively. The
concentrations of 6 of the 14 halogenated gases studied are decreasing as a
result of the Montreal Protocol.

A number of factors affect the accuracy of the results. First, the
calculations used the absorption cross sections given by the Hitran database
which have an uncertainty of several percent. Second, the Hitran data does not
show a temperature dependence of the cross sections. This contrasts with our
previous work [16] that examined the effect of temperature on the frequency
integrated cross section which equals the sum of the line strengths.
Negligible temperature effect was found for the 5 main greenhouse gases as
well as CF4 but the absorption cross section for SF6 increased by a factor of
2 as temperature decreased from 300 to 200 K. Molecules made of many atoms,
like SF6, have many low-frequency, infrared-inactive modes which can be
significantly excited at higher temperatures [16]. Then the molecule spends
more time in nonradiating modes and emits less average power at the
frequencies of infrared active modes than a simpler molecule would. The
temperature increase causes a corresponding decrease in the absorption cross
section of the infrared active mode. The absorption cross sections of other
complicated halogenated moecules can be expected to have analogous temperature
dependences. Finally, the altitude dependence of the halogenated gas
concentrations is poorly known. For the cases of CFC11, CFC12 and SF6, the
radiative forcing increases by up to 14% if the altitude concentration
dependence is ignored.

Our calculations were done for the case of a clear sky. In reality, about 70%
of the Earth’s surface is covered by clouds [24]. A cloud layer will
particularly perturb the so called spectral window shown between the CO2 and
O3 absorptions in Fig. 4, where the halogenated gas absorptions most affect
the radiative forcing. It is therefore reasonable that incorporating clouds,
will reduce the halogenated gas radiative forcing. Some studies have found
clouds reduce the halogenated gas forcing by 25% [7].

This study computed the global warming potentials which are comparable to
those found by the IPCC2021 as shown in Table 5. It also evaluated the forcing
power per molecule defined as the derivative of the radiative forcing of a gas
with respect to its column density. The advantage of the per-molecule forcing
power (20) is that unlike the GWP it enables one to quickly estimate the
relative warming caused by a concentration change of a greenhouse gas. The
forcing power strongly depends on the gas concentration. At very low
concentrations, all gases have an optically thin power on the order of
$10^{-21}$ W per molecule. Table 3 shows the forcing power per molecule of the
halogenated gases decreases slightly at higher gas concentrations. However,
for the 5 main greenhouse gases which have much higher gas concentrations, the
forcing power per molecule decreases by orders of magnitude. The relative
radiative forcing change caused by a change in the concentration of each gas
is given in Table 4. Most of the total radiative forcing change is caused by
increases in CO2 (83%), CH4 (8%) and N2O (6%). The radiative forcing change
due to the 14 halogenated molecules is only 2% of the total. Assuming the
temperature change due to a changing gas concentration is proportional to its
radiative forcing change, the halogenated gases are responsible for a surface
warming of 0.002 C per decade compared to the total observed warming of about
0.1 C per decade.

## Acknowledgements

The Canadian Natural Science and Engineering Research Council provided
financial support of one of us.

## References

  * [1] M. J. Molina and F. S. Rowland, Stratospheric sink for chlorofluoromethanes: chlorine atom-catalysed destruction of ozone, Nature, 249 810-812 (1974).
  * [2] P. J. Crutzen, Estimates of possible future ozone reductions from continued use of fluoro-chloromethanes (CF2 Cl2, CFCl3), Geophys. Res. Lett. 1, No. 5, 205-208 (1974).
  * [3] About Montreal Protocol, https://www.unep.org/ozonaction/who-we-are/about-montreal-protocol, Accessed 31 May 2023.
  * [4] V. Masson-Delmotte et al, Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, UK, 2391 pp. doi: 10.1017/9781009157896.
  * [5] A. K. Jain, B. P. Briegleb, K. Mischwaner and D. J. Wuebbles, Radiative forcings and global warming potentials of 39 greenhouse gases, J. Geophys. Res. 105, No. D16, 20773-20790 (2000).
  * [6] K. Sihra, M. D. Hurley, K. P. Shine and T. J. Wallington, Updated radiative forcing estimates of 65 halocarbons and nonmethane hydrocarbons, J. Geophys. Res. 106, No. D17, 20493-20505 (2001).
  * [7] S. A. Buehler, M. Brath, O. Lemke, O. Hodnebrog, R. Pincus, P. Eriksson, I. Gordon and R. Larsson, A New Halocarbon Absorption Model Based on HITRAN Cross-Section Data and New Estimates of Halocarbon Instantaneous Clear-Sky Radiative Forcing, J. Advances in Modeling Earth Systems, 10.1029/2022MS003239 (2022).
  * [8] I. E. Gordon et al., The HITRAN2016 Molecular Spectroscopic Database, JQSRT 203, 3-69 (2017) https://hitran.org/ Accessed April 1, 2023.
  * [9] The U.S. Standard Atmosphere, NASA Report TM-X-74335 (1976).
  * [10] W. A. van Wijngaarden and W. Happer, Dependence of Earth’s Thermal Radiation on Five Most Abundant Greenhouse Gases, Atmos. & Oceanic Phys. arXiv: 2006.03098 (2020).
  * [11] G. P. Anderson, S. A. Clough, F. X. Kneizys, J. H. Chetwynd and E. P. Shettle, AFGL Atmospheric Constituent Profiles (0-120 km), AFGL-TR-86-0110 Air Force Geophysics Laboratory, Hanscom Air Force Base, Massachusetts, (1986).
  * [12] Global Greenhouse Gas Reference Network, NOAA Earth System Research Laboratory, https://gml.noaa.gov/aftp/data/hats/solvents/CCl4/combined/HATS-global-CCl4.txt, Accessed May 10, 2023.
  * [13] M. K. Vollmer et al, Atmospheric histories and emissions of chlorofluorocarbons CFC-13 (CClF3), $\Sigma$CFC-114 (C2Cl2F4) and CFC-115 (C2ClF5), Atmos. Chem. Phys. 18, 979-1002, (2018).
  * [14] K. M. Stanley, D. Say, J. Mühle, C. M. Harth, P. B. Krummel, D. Young, S. J. O’Doherty, P. K. Salamey, P. G. Simmonds, R. F. Weiss, R. G. Prinn, P. J. Fraser & M. Rigby, Increase in global emissions of HFC-23 despite near-total expected reductions, Nature Comm. 11, 397 (2020), doi.org/10.138/s41467-019-13899-4.
  * [15] A. T. Brown, C. M. Volk, M. R. Schoebert, C. D. Bonne & P. F. Bernath, Stratospheric lifetimes of CFC-12, CCl4, CH4, CH3Cl and N2O from measurements made by the Atmospheric Chemistry Experiment- Fourier Transform Spectrometer (ACE-FTS), Atmos. Chem. Phys. 13, 6921-6950 (2013).
  * [16] W. A. van Wijngaarden and W. Happer, Relative Potency of Greenhouse Molecules, Atmos. & Oceanic Phys. arXiv: 2103.16465 (2021)
  * [17] K. Schwarzschild, Über Diffusion und Absorption in der Sonnenatmosphäre, Nachr. K. Gesell. Wiss. Math.-Phys. Klasse 195, 41 (1906).
  * [18] J. J. Harrison, New and improved infrared absorption cross sections for trichlorofluoromethane (CFC-11), Atmos. Meas. Tech. 11, 5827-5836 (2018).
  * [19] R. Spencer, Global Warming, https://www.drroyspencer.com/latest-global-temperatures/, Accessed Oct. 2022.
  * [20] K. de Lange, J. D. Ferguson, W. Happer and W. A. van Wijngaarden, Nitrous Oxide and Climate, Atmos. & Oceanic Phys. arXiv: 2207.093978 (2022).
  * [21] L. D. Harvey, A guide to global warming potentials (GWPs), Energy Policy, p. 23, January 1993.
  * [22] E. Joos et al., Carbon dioxide and climate impulse response functions for the computation of greenhouse gas metrics: a multi-model analysis, Atmos. Chem. Phys. 13, 2293-2825, (2013).
  * [23] C. Smith, Z. Nicholls, K. Armour, W. Collins, P. Forster, M. Meinshausen, M. Palmer and M. Watanabe, IPCC 2021 Chapter 7: The Earth’s energy budget, climate feedbacks, and climate sensitivity - Supplementary Material, (2021).
  * [24] International Satellite Cloud Climatology Project (ISCCP) https://isccp.giss.nasa.gov/ Accessed 13 June 2023.