# QCD axion bubbles from the hidden SU(N) gauge symmetry breaking

Hai-Jun Li<EMAIL_ADDRESS>Center for Advanced Quantum Studies,
Department of Physics, Beijing Normal University, Beijing 100875, China

###### Abstract

The QCD axion bubbles can be formed due to an extra Peccei-Quinn (PQ) symmetry
breaking in the early Universe. In this paper, we investigate the QCD axion
bubbles formation from the PQ symmetry broken by hidden $SU(N)_{H}$ gauge
interactions after inflation, which leads to the multiple vacua. The axion
acquires a light mass and then settles down into different vacua. The QCD
axion bubbles are formed when the conventional QCD axion arises during the QCD
phase transition. In our scenario, the QCD axions that start to oscillate at
the large values $\sim 2\pi/3$ can lead to the high density axion bubbles with
$N=2$. The cosmological implications of the QCD axion bubbles are also
discussed, such as the primordial black holes (PBHs) and the axion
miniclusters. We find that the PBH mass is lager than $\sim\mathcal{O}(5\times
10^{5})M_{\odot}$ for the axion scale $f_{a}\sim\mathcal{O}(10^{16})\,\rm
GeV$.

††preprint: BNU-23-028

##  I Introduction

The strong CP violation in quantum chromodynamics (QCD) is a long-standing
problem and can be solved by the Peccei-Quinn (PQ) mechanism with a
spontaneously broken $U(1)$ PQ symmetry Peccei and Quinn (1977a, b). The
mechanism predicts the light pseudo Nambu-Goldstone (NG) boson axion Weinberg
(1978); Wilczek (1978), also called the QCD axion, which acquires a tiny mass
from the QCD non-perturbative effects ’t Hooft (1976a, b). When the potential
of QCD axion is generated by the QCD instanton, axion is stable at the minimum
value of CP conservation, which solves the strong CP problem. The QCD axion is
the potential cold dark matter (DM) candidate if non-thermally produced in the
early Universe through misalignment mechanism Preskill  _et al._ (1983);
Abbott and Sikivie (1983); Dine and Fischler (1983); Preskill  _et al._
(1983). The axion is massless at high temperatures, as the cosmic temperature
decreases, it acquires a non-zero mass at the QCD phase transition and starts
to oscillate when its mass becomes comparable to the Hubble parameter, which
explains the observed DM abundance. See $\rm e.g.$ Refs. Di Luzio  _et al._
(2020); Chadha-Day  _et al._ (2022); Adams  _et al._ (2022) for recent
reviews.

The PQ symmetry is supposed to be broken before or during inflation. In this
case, the QCD axion is massless during inflation and acquires the quantum
fluctuations. To suppress the isocurvature perturbation, a feasible method is
considering an extra PQ symmetry breaking Dine and Anisimov (2005); Jeong and
Takahashi (2013); Takahashi and Yamada (2015); Harigaya  _et al._ (2015). The
PQ symmetry is an approximate global symmetry and can be strongly broken in
the early Universe with the multiple vacua Kallosh  _et al._ (1995); Banks and
Seiberg (2011); Witten (2018); Harlow and Ooguri (2019); Jeong  _et al._
(2022). The QCD axion bubbles can be formed due to this extra PQ symmetry
breaking after inflation Kitajima and Takahashi (2020). In this case, axion
acquires a light mass and oscillates when this mass is greater than the Hubble
parameter. Note that the extra PQ symmetry is temporarily broken, the axion
potential will disappear before the QCD phase transition. Therefore, the final
QCD axion abundance can also be calculated through misalignment mechanism.
However, the initial misalignment angle is determined by the extra PQ symmetry
breaking, which can be split into different values. At the QCD scale, the QCD
axions start to oscillate with these initial angles, in which the high axion
density region can be formed with the large initial angle values. This high
axion density region is called the “axion bubble”. The concept of QCD axion
bubbles was first proposed in Ref. Kitajima and Takahashi (2020), which is
similar to the baryon bubbles in the inhomogeneous Affleck-Dine baryogenesis
Dolgov  _et al._ (2009); Hasegawa and Kawasaki (2019). They also studied the
formations of primordial black holes (PBHs) and axion miniclusters from the
QCD axion bubbles.

In this paper, we focus our attention on the QCD axion bubbles and investigate
the axion bubbles formation due to the extra PQ symmetry dynamically broken by
the hidden $SU(N)_{H}$ gauge interactions in the early Universe. After
inflation, the axion acquires a light mass due to this PQ symmetry breaking
with a large axion scale, and then settles down into the nearest potential
minimum. The multiple vacua split the initial misalignment angle into
different values. In this case, if the axion initial value is smaller than a
critical value, the axion is stabilized near the origin. On the contrary, the
axion will be stabilized at the another minimum if the axion initial value is
larger than the critical value. During the QCD phase transition, the QCD
axions start to oscillate near the origin and the another large value. The
former case can account for the cold DM abundance, and the latter one can form
the large energy density QCD axion bubbles. In our scenario, the QCD axions
that start to oscillate at the large values $\sim 2\pi/3$ can lead to the high
density axion bubbles with $N=2$. We also discuss the cosmological
implications of the QCD axion bubbles, focusing on the PBHs formation. The
PBHs are formed when the axions dominate the radiation in the bubbles, which
predicts the minimum PBH mass $\sim\mathcal{O}(5\times 10^{5})M_{\odot}$ for
the axion scale $f_{a}\sim\mathcal{O}(10^{16})\,\rm GeV$. On the other hand,
if the bubbles enter the horizontal plane before the axions dominate the
radiation, they will eventually form the axion miniclusters.

This paper is organized as follows. In Sec. II, we briefly review the QCD
axion and the misalignment mechanism. In Sec. III, we discuss the QCD axion
bubbles formation due to the extra PQ symmetry dynamically broken by the
hidden $SU(N)_{H}$ gauge interactions in the early Universe. The cosmological
implications are discussed in Sec. IV, including the axion DM abundance, PBHs,
and axion miniclusters. The discussion and conclusion are given in Sec. V.

##  II QCD axion and misalignment mechanism

Here we briefly review the QCD axion and the misalignment mechanism. The QCD
axion is a pseudo NG boson with a spontaneously broken $U(1)$ PQ symmetry. It
couples to gluons with the following effective Lagrangian

$\displaystyle\mathcal{L}_{agg}=-\frac{\alpha_{s}}{8\pi}\frac{a}{f_{a}}G^{a\,\mu\nu}\tilde{G}_{\mu\nu}^{a}\,,$
(1)

where $\alpha_{s}$ is the strong fine structure constant, $a$ is the axion
field, $f_{a}$ is the axion decay constant, $\theta=a/f_{a}$ is the axion
angle, $G^{a\,\mu\nu}$ and $\tilde{G}_{\mu\nu}^{a}$ are the gluon field
strength tensor and its dual tensor, respectively. The resulting effective
potential of QCD axion is given by

$\displaystyle V_{\rm
QCD}(a)=m_{a}^{2}(T)f_{a}^{2}\left[1-\cos\left(\frac{a}{f_{a}}\right)\right]\,,$
(2)

where $m_{a}(T)$ is the temperature-dependent axion mass for $T\gtrsim T_{\rm
QCD}$ ($\sim 150\,\rm MeV$) Borsanyi  _et al._ (2016)

$\displaystyle m_{a}(T)\simeq m_{a,0}\left(\frac{T}{T_{\rm
QCD}}\right)^{-4.08}\,,$ (3)

with the zero-temperature axion mass Grilli di Cortona  _et al._ (2016)

$\displaystyle m_{a,0}\simeq 5.70(7)\,{\mu\rm
eV}\left(\frac{f_{a}}{10^{12}\,{\rm GeV}}\right)^{-1}\,.$ (4)

As the cosmic temperature decreases, the QCD axion starts to oscillate when
its mass $m_{a}(T)$ becomes comparable to the Hubble parameter $H(T)$,
$3H(T_{a})=m_{a}(T_{a})$, with the oscillation temperature

$\displaystyle\begin{aligned} T_{a}&\simeq 9.59\times 10^{-1}\,{\rm
GeV}\left(\frac{g_{*}(T_{a})}{61.75}\right)^{-0.082}\\\
&\times\left(\frac{f_{a}}{10^{12}\,\rm GeV}\right)^{-0.16}\,,\end{aligned}$
(5)

where $g_{*}(T)$ is the number of effective degrees of freedom. The axion
number density at $T_{a}$ is given by

$\displaystyle
n_{a}(T_{a})=\frac{1}{2}m_{a}(T_{a})f_{a}^{2}\left\langle\theta_{i}^{2}f(\theta_{i})\right\rangle\chi\,,$
(6)

where $\theta_{i}$ is the initial misalignment angle, $\chi\simeq 1.44$ is a
numerical factor Turner (1986), and $f(\theta_{i})$ is the anharmonic factor

$\displaystyle
f(\theta_{i})\simeq\left[\ln\left(\frac{e}{1-\theta_{i}^{2}/\pi^{2}}\right)\right]^{1.16}\,,$
(7)

which is taken as $f(\theta_{i})\simeq 1$ for $|\theta_{i}|\ll\pi$ Lyth
(1992); Visinelli and Gondolo (2009). The present axion energy density
$\rho_{a}(T_{0})=m_{a,0}n_{a}(T_{0})$ is

$\displaystyle\rho_{a}(T_{0})=\frac{m_{a,0}m_{a}(T_{a})s(T_{0})}{2s(T_{a})}f_{a}^{2}\left\langle\theta_{i}^{2}f(\theta_{i})\right\rangle\chi\,,$
(8)

where $s(T)=2\pi^{2}g_{*}(T)T^{3}/45$ is the entropy density, and $T_{0}$ is
the present CMB temperature. Then we have the current QCD axion abundance
$\Omega_{a}h^{2}=\rho_{a}(T_{0})/\rho_{c}h^{2}$ as

$\displaystyle\begin{aligned} \Omega_{a}h^{2}&\simeq 1.43\times
10^{-1}\left(\frac{g_{*}(T_{0})}{3.94}\right)\left(\frac{g_{*}(T_{a})}{61.75}\right)^{-0.42}\\\
&\times\left(\frac{f_{a}}{10^{12}\,\rm
GeV}\right)^{1.16}\left\langle\theta_{i}^{2}f(\theta_{i})\right\rangle\,,\end{aligned}$
(9)

where $\rho_{c}=3H_{0}^{2}M_{\rm Pl}^{2}$ is the critical energy density,
$M_{\rm Pl}$ is the reduced Planck mass, and $h\simeq 0.68$ is the reduced
Hubble constant. In order to explain the observed cold DM abundance,
$\Omega_{\rm DM}h^{2}\simeq 0.12$ Aghanim  _et al._ (2020), we derive the
initial misalignment angle

$\displaystyle\begin{aligned} \theta_{i}&\simeq
0.87\left(\frac{g_{*}(T_{0})}{3.94}\right)^{1/2}\left(\frac{g_{*}(T_{a})}{61.75}\right)^{0.21}\\\
&\times\left(\frac{f_{a}}{10^{12}\,\rm GeV}\right)^{-0.58}\,.\end{aligned}$
(10)

Note that the scope of above discussion is $f_{a}\lesssim 10^{17}\,\rm GeV$.
We can see that the abundance of QCD axion will be many orders of magnitude
larger than that of the observed DM with high scales
$f_{a}\sim\mathcal{O}(10^{16}-10^{17})\,\rm GeV$, unless the initial
misalignment angle is less than $\sim\mathcal{O}(10^{-2})$.

##  III QCD axion bubbles

In this section, we discuss the QCD axion bubbles formation. The QCD axion
bubbles can be formed due to an extra PQ symmetry breaking after inflation,
and are formed at the QCD phase transition. The evolution can be summarized as
the following three parts.

  * •

During/After inflation.
During inflation, the QCD axion is massless and acquires the quantum
fluctuations. An extra PQ symmetry breaking after inflation is required. Then
the axion acquires a light mass due to this extra PQ symmetry breaking, which
can lead to the multiple vacua. The axion starts to oscillate when its mass is
comparable to the Hubble parameter, and then settles down into different
minimum depending on its initial position. Note that this extra symmetry is
temporarily broken, and the resulting axion potential will be disappeared
before the QCD phase transition.

  * •

During/After the QCD phase transition.
The QCD axion bubbles are formed when the QCD axion starts to oscillate at the
QCD scale. During the QCD phase transition, the conventional axion potential
$V_{\rm QCD}(\phi)$ arises. Therefore, one of the multiple vacua minimum
should be around at the minimum of $V_{\rm QCD}(\phi)$ to ensure the correct
cold DM abundance. In this case, if the axion initial value is smaller than a
critical value, the axion is stabilized near the origin. Otherwise, the axion
will be stabilized at the another minimum with a large initial value, leading
to the high density QCD axion bubbles.

  * •

The late time.
The QCD axion bubbles can lead to some further interesting phenomena, such as
the formations of PBHs and axion miniclusters, which will be briefly discussed
in Sec. IV.

###  III.1 PQ symmetry broken by hidden $SU(N)_{H}$

There are many scenarios for the extra PQ symmetry breaking in the early
Universe, such as the Witten effect of monopoles in hidden sectors Nomura  _et
al._ (2016); Kawasaki  _et al._ (2016, 2018), a larger scale of the
spontaneous PQ symmetry breaking with the higher dimensional term Chiba  _et
al._ (2004); Takahashi and Yamaguchi (2004); Higaki  _et al._ (2014); Co  _et
al._ (2020) and the hidden non-Abelian gauge interactions Takahashi and Yamada
(2015), and a stronger QCD with the large Higgs field expectation value Choi
_et al._ (1997); Banks and Dine (1997); Jeong and Takahashi (2013). In Ref.
Kitajima and Takahashi (2020), they considered the Witten effect as an example
of the extra PQ symmetry breaking to form the QCD axion bubbles. For
comparison, we consider the extra PQ symmetry is broken by the hidden non-
Abelian gauge interactions whose scales are temporarily enhanced by the
dynamics of flat directions, which was considered in Ref. Takahashi and Yamada
(2015) to suppress the isocurvature perturbations of the QCD axion. In this
case, the PQ symmetry is dynamically broken by the hidden $SU(N)_{H}$ gauge
interactions. Then the axion becomes so heavy during inflation and its
isocurvature perturbations can be significantly suppressed. However, in our
scenario the PQ symmetry is assumed to be broken after inflation with the
large axion scale, resulting in a light axion mass and the multiple vacua.

In Ref. Takahashi and Yamada (2015), they considered a supersymmetric (SUSY)
axion model, but this can be applied to the non-SUSY case. We begin with the
superpotential

$\displaystyle W=y_{ij}\psi
Q_{H}^{i}\bar{Q}_{H}^{j}+{y}_{kl}^{\prime}\phi{Q}_{H}^{{}^{\prime}k}{\bar{Q}}_{H}^{{}^{\prime}l}+W_{\cancel{\rm
PQ}}+W_{\phi}(\phi)\,,~{}$ (11)

where contains $N_{F}$ ($N^{\prime}_{F}$) flavors of the hidden quarks $Q_{H}$
($Q^{\prime}_{H}$) and anti-quarks $\bar{Q}_{H}$ ($\bar{Q}^{\prime}_{H}$) in
the representation of $SU(N)_{H}$ with (without) the PQ charges, respectively,
$\psi$ is the PQ breaking scalar field, $\phi$ is the singlet scalar field,
$W_{\cancel{\rm PQ}}$ is the PQ breaking term, and $y$ is the Yukawa coupling
with the flavor indices $i$, $j$, $k$, and $l$. Assuming that $\phi$ is a flat
direction, then the superpotential of $\phi$ is given by

$\displaystyle W_{\phi}(\phi)=\lambda\frac{\phi^{4}}{4M_{\rm Pl}}\,,$ (12)

where $\lambda\sim\mathcal{O}(1)$. The potential of $\phi$ during inflation is

$\displaystyle\begin{aligned}
V(\phi)&=\left(m_{\phi}^{2}-c_{H}H^{2}\right)|\phi|^{2}\\\
&-\left(a_{H}H\lambda\frac{\phi^{4}}{4M_{\rm Pl}}+{\rm
c.c.}\right)+\lambda^{2}\frac{|\phi|^{6}}{M_{\rm Pl}^{2}}\,,\end{aligned}$
(13)

where $c_{H}=a_{H}=1$, $m_{\phi}$ is the soft SUSY breaking mass, the three
terms on the right-hand side represent the sum of low-energy and Hubble-
induced soft mass terms, the Hubble-induced $A$-terms, and the contribution of
the non-renormalizable term, respectively. The flat direction $\phi$ acquires
a mass of order the Hubble parameter during inflation, and starts to oscillate
when its soft mass becomes comparable to the Hubble parameter. Then it settles
down into the potential minimum

$\displaystyle\langle|\phi|\rangle(t)\simeq\left(\frac{H(t)M_{\rm
Pl}}{\sqrt{3}\lambda}\right)^{1/2}\,.$ (14)

Note that the phase direction of $\phi$ also has a mass, which comes from the
Hubble-induced $A$-terms. The $F$-component of $\phi$ is given by

$\displaystyle\langle F_{\phi}\rangle\simeq\left(\frac{H^{3}(t)M_{\rm
Pl}}{3^{3/2}\lambda}\right)^{1/2}\,,$ (15)

which hints that the $SU(N)_{H}$ gaugino acquires a soft mass due to the gauge
mediated SUSY breaking

$\displaystyle
m_{\lambda}\simeq\frac{N^{\prime}_{F}g_{H}^{2}}{16\pi^{2}}\frac{F_{\phi}}{\langle\phi\rangle}\simeq\frac{N^{\prime}_{F}g_{H}^{2}}{16\sqrt{3}\pi^{2}}H(t)\,,$
(16)

where $g_{H}$ is the gauge coupling.

(a) $N=2$. (b) $N=3$.

Figure 1: Illustration of the QCD axion bubbles formation with the axion
potentials. The left and right panels correspond to the cases $N=2$ and $3$,
respectively. The red and blue lines represent the potentials $V_{\rm
QCD}(\phi)$ and $V_{\cancel{\rm PQ}}(\phi)$ with a small $\theta_{i}$,
respectively. The amplitude of $V_{\cancel{\rm PQ}}(\phi)$ and the value of
$\theta_{i}$ are exaggerated for illustrative purposes.

The Lagrangian of the $SU(N)_{H}$ gauge field is Takahashi and Yamada (2015)

$\displaystyle\mathcal{L}=\frac{1}{2}\int{\rm
d}^{2}\theta\left(\frac{1-2m_{\lambda}\theta^{2}}{2g_{H}^{2}}-\frac{i\theta_{H}}{16\pi^{2}}\right)W^{\alpha}W_{\alpha}+{\rm
h.c.}\,,~{}~{}~{}$ (17)

where $\theta_{H}$ is the vacuum angle. The axion comes from the gauge kinetic
function with the notation

$\displaystyle\theta_{H}\to\theta_{H}-3N_{F}\frac{\phi}{v_{a}}\,.$ (18)

Using $f_{a}=v_{a}/N_{\rm DW}=v_{a}/(NN_{F})$, we have

$\displaystyle\theta_{H}\to\theta_{H}-\frac{3}{N}\frac{\phi}{f_{a}}\,,$ (19)

where $N_{\rm DW}$ is the domain wall number. The gaugino condensation induces
the effective superpotential

$\displaystyle W_{\rm eff}(\phi)=N\tilde{\Lambda}_{H}^{3}(\phi)\,,$ (20)

where $\Lambda_{H}(\phi)$ is the condensation scale. Then the axion field
acquires the effective potential

$\displaystyle\begin{aligned} V_{\cancel{\rm PQ}}(\phi)&=-\int{\rm
d}^{2}\theta W_{\rm eff}(\phi)+{\rm h.c.}\\\
&=\frac{32\pi^{2}}{g_{H}^{2}}m_{\lambda}\Lambda_{H}^{3}(\phi)\cos\left(\frac{\theta_{H}}{N}-\frac{3}{N^{2}}\frac{\phi}{f_{a}}\right)+\cdot\cdot\cdot\,,~{}~{}\end{aligned}$
(21)

and the axion effective mass is given by Takahashi and Yamada (2015)

$\displaystyle
m_{a}^{2}(\phi)=\frac{\hat{c}H(t)\Lambda_{H}^{3}(\phi)}{f_{a}^{2}}\,,$ (22)

where $\hat{c}=6\sqrt{3}N^{\prime}_{F}/N^{4}$.

It should be emphasized that in our scenario the PQ symmetry is assumed to be
broken after inflation with the large axion scale $f_{a}$. The axion acquires
a light effective mass and starts to oscillate when its mass becomes
comparable to the Hubble parameter. Then the axion will be stabilized at the
potential minimum with

$\displaystyle\phi_{\rm min}^{n}=\left(\theta_{H}+2\pi
nN\right)\frac{N}{3}f_{a}\,,$ (23)

where $n$ is the integer, which represents the multiple vacua.

###  III.2 Axion bubbles formation

Here we show the QCD axion bubbles formation due to above PQ symmetry breaking
in the early Universe. During inflation, the QCD axion is massless and
acquires the quantum fluctuations with a Gaussian distribution. After
inflation, the axion acquires a light axion mass due to the PQ symmetry broken
by the hidden $SU(N)_{H}$ gauge interactions and settles down into the nearest
potential minimum $\phi_{\rm min}^{n}$. During the QCD phase transition, the
conventional axion potential $V_{\rm QCD}(\phi)$ arises. Therefore, one of the
minimum of $\phi_{\rm min}^{n}$ should be around at the minimum of $V_{\rm
QCD}(\phi)$ to ensure the correct DM abundance. In this paper, we consider the
large axion scale $f_{a}\sim\mathcal{O}(10^{16})\,\rm GeV$ and assume that it
is unchanged. In order to explain the DM abundance, we derive the small
initial misalignment angle

$\displaystyle\begin{aligned} \theta_{i}&\simeq 4.29\times
10^{-3}\left(\frac{g_{*}(T_{0})}{3.94}\right)^{1/2}\left(\frac{g_{*}(T_{a})}{61.75}\right)^{0.21}\\\
&\times\left(\frac{f_{a}}{10^{16}\,\rm GeV}\right)^{-0.58}\,.\end{aligned}$
(24)

Note that this initial angle is corresponding to the QCD axion.

We first consider the case with $N=2$, see Fig. 1 (a). The red and blue lines
represent the potentials $V_{\rm QCD}(\phi)$ and $V_{\cancel{\rm PQ}}(\phi)$,
respectively. After inflation, the axion is stabilized at the potential
minimum of $V_{\cancel{\rm PQ}}(\phi)$ with

$\displaystyle\phi_{\rm min}^{0}\simeq 0\,,\quad\phi_{\rm
min}^{1}\simeq\frac{8\pi}{3}f_{a}\,,\quad\cdot\cdot\cdot\,,$ (25)

which correspond to the effective initial misalignment angle $\theta_{i,n}$
with

$\displaystyle\theta_{i,0}=0-\theta_{i}\,,\quad\theta_{i,1}=\frac{2\pi}{3}-\theta_{i}\,,\quad\cdot\cdot\cdot\,.$
(26)

When the QCD axion $V_{\rm QCD}(\phi)$ arises during the QCD phase transition,
the state $\phi_{\rm min}^{0}$ with the effective initial angle $\theta_{i,0}$
explains the DM. On the other hand, if the initial value is larger than
$\phi_{\rm crit}$, the axion will settle down into the minimum $\phi_{\rm
min}^{1}$ with the initial angle $\theta_{i,1}$. In this case, since the large
initial misalignment angle $\sim 2\pi/3$, the local axion density at the
minimum $\phi_{\rm min}^{1}$ becomes much higher than that at $\phi_{\rm
min}^{0}$. This high axion density region is called the “axion bubble”.

When the axion dominates the radiation in the bubbles, $\rm i.e.$, the axion
energy density is equal to the radiation energy density, we define the cosmic
background temperature, $T_{B}$. The local axion energy density in the bubbles
at $T_{B}$ is given by

$\displaystyle\rho_{a}(T_{B})=\frac{m_{a,0}m_{a}(T_{a})s(T_{B})}{2s(T_{a})}f_{a}^{2}\left\langle\theta_{i,1}^{2}f(\theta_{i,1})\right\rangle\chi\,.$
(27)

The radiation energy density is defined by

$\displaystyle\rho_{R}(T_{B})=\frac{\pi^{2}}{30}g_{*}(T_{B})T_{B}^{4}\,.$ (28)

Considering $\rho_{a}(T_{B})=\rho_{R}(T_{B})$, and the effective initial angle
$\theta_{i,1}=2\pi/3-\theta_{i}$, we derive the temperature

$\displaystyle T_{B}\simeq 0.23\,{\rm
MeV}\left(\frac{g_{*}(T_{a})}{61.75}\right)^{-0.42}\left(\frac{f_{a}}{10^{16}\,\rm
GeV}\right)^{1.16}\,.~{}~{}$ (29)

Note that here the initial angle is taken as $\theta_{i}\simeq 4.29\times
10^{-3}$ for $f_{a}=10^{16}\,\rm GeV$ and we do not make other
simplifications.

Additionally, we also consider the case with $N=3$, which is shown in Fig. 1
(b). In this case, the axion will be stabilized at the potential minimum with

$\displaystyle\phi_{\rm min}^{0}\simeq 0\,,\quad\phi_{\rm min}^{1}\simeq 6\pi
f_{a}\,,\quad\cdot\cdot\cdot\,,$ (30)

corresponding to the effective initial angle $\theta_{i,n}$ with

$\displaystyle\theta_{i,0}=0-\theta_{i}\,,\quad\theta_{i,1}=0-\theta_{i}\,,\quad\cdot\cdot\cdot\,.$
(31)

Since the effective initial angle $\theta_{i,0}=\theta_{i,1}\simeq 0$, $\rm
i.e.$, the local axion density at $\phi_{\rm min}^{0}$ is same with that at
$\phi_{\rm min}^{1}$, we find no QCD axion bubbles is generated in this case.
Therefore, in the following we just discuss the axion bubbles in the case with
$N=2$.

Then we discuss the QCD axion bubbles abundance, which is related to the
inflationary fluctuations. The probability density function obeys the Fokker-
Planck equation Starobinsky (1982); Linde (1982); Starobinsky and Yokoyama
(1994). Under the initial condition $P(N,\phi)=\delta(\phi-\phi_{i})$, and
assuming the Hubble parameter during inflation is constant, its solution can
be given by the Gaussian distribution Kitajima and Takahashi (2020)

$\displaystyle
P(N,\phi)=\frac{1}{\sqrt{2\pi}\sigma(N)}\exp{\left(-\frac{\left(\phi-\phi_{i}\right)^{2}}{2\sigma^{2}(N)}\right)}\,,$
(32)

with the variance

$\displaystyle\sigma(N)=\frac{H_{\rm inf}}{2\pi}\sqrt{N}\,.$ (33)

Note that here $N$ is the e-folding number, and $H_{\rm
inf}\sim\mathcal{O}(10^{13})\,\rm GeV$ is the Hubble parameter during
inflation. Then the volume fraction of the QCD axion bubbles is given by
Kitajima and Takahashi (2020)

$\displaystyle\begin{aligned} \frac{{\rm d}\beta}{{\rm
d}N}&\simeq\frac{\partial}{\partial N}\int_{\phi_{\rm crit}}^{\infty}{\rm
d}\phi P(N,\phi)\\\ &=\frac{1}{2}\frac{\partial}{\partial N}{\rm
erfc}\left(\frac{\phi_{\rm crit}-\phi_{i}}{\sqrt{2}\sigma(N)}\right)\\\
&=\frac{\phi_{\rm crit}-\phi_{i}}{2}P(N,\phi_{\rm crit})\,,\end{aligned}$ (34)

where ${\rm erfc}(x)$ is the complementary error function. With ${\rm d}N={\rm
d}\ln k$, where $k$ is the wave number, then we obtain the volume fraction

$\displaystyle\frac{{\rm d}\beta}{{\rm d}\ln k}\simeq\frac{\phi_{\rm
crit}-\phi_{i}}{2}P(\ln(k/k_{*})+N_{*},\phi_{\rm crit})\,,$ (35)

where $k_{*}=0.002\,\rm Mpc^{-1}$ and $N_{*}=0$.

##  IV Cosmological implications

In this section, we discuss the cosmological implications of the QCD axion
bubbles, including the axion DM abundance, PBHs, and axion miniclusters.

###  IV.1 Axion dark matter abundance

The QCD axion abundance depends on the initial misalignment angle $\theta_{i}$
and the decay constant $f_{a}$. In our scenario, the current axion abundance
can also be described by the misalignment mechanism with the effective initial
angle $\theta_{i,0}$. For the large scale $f_{a}\sim\mathcal{O}(10^{16})\,\rm
GeV$, we have the current QCD axion abundance

$\displaystyle\begin{aligned} \Omega_{a}h^{2}&\simeq 6.52\times
10^{3}\left(\frac{g_{*}(T_{0})}{3.94}\right)\left(\frac{g_{*}(T_{a})}{61.75}\right)^{-0.42}\\\
&\times\left(\frac{f_{a}}{10^{16}\,\rm
GeV}\right)^{1.16}\left\langle\theta_{i,0}^{2}f(\theta_{i,0})\right\rangle\,.\end{aligned}$
(36)

In order to explain the observed cold DM abundance, the effective initial
angle is given by Eq. 24 with $|\theta_{i,0}|\sim\mathcal{O}(10^{-3})$.

###  IV.2 Primordial black holes & axion miniclusters

The PBHs can be generated due to large density perturbations in the early
Universe, they are also the DM candidates Bird  _et al._ (2016); Carr and
Kuhnel (2020); Green and Kavanagh (2021); Carr and Kuhnel (2022); Li (2022).
The initial PBH mass at the formation time $t_{f}$ is given by $M_{\rm
PBH}=4\pi\gamma\rho_{R}/(3H_{f}^{3})$ Carr  _et al._ (2010), where
$\gamma\simeq 0.2$ is the gravitational collapse factor Carr (1975),
$\rho_{R}$ is the radiation energy density, and $H_{f}$ is the Hubble
parameter at $t_{f}$. Then we derive the PBH mass at the formation time

$\displaystyle\begin{aligned} \frac{M_{\rm PBH}}{M_{\odot}}&\simeq
0.03\times\left(\frac{\gamma}{0.2}\right)\left(\frac{g_{*}(T_{f})}{10.75}\right)^{-1/2}\\\
&\times\left(\frac{T_{f}}{1\,\rm GeV}\right)^{-2}\,,\end{aligned}$ (37)

where $T_{f}$ is the temperature at $t_{f}$ and $M_{\odot}$ is the solar mass.
We also have the relation Di and Gong (2018)

$\displaystyle\begin{aligned} \frac{M_{\rm PBH}}{M_{\odot}}&\simeq
3.68\times\left(\frac{\gamma}{0.2}\right)\left(\frac{g_{*}(T_{f})}{10.75}\right)^{-1/6}\\\
&\times\left(\frac{k}{10^{6}\rm\,Mpc^{-1}}\right)^{-2}\,,\end{aligned}$ (38)

where $k$ is the wave number.

Figure 2: The PBH fractional abundance $f_{\rm PBH}$ as a function of the PBH
mass (in the solar mass $M_{\odot}$). The red, orange, and blue lines
correspond to the results with $f_{a}=5\times 10^{16}\,\rm GeV$, $1\times
10^{16}\,\rm GeV$, and $5\times 10^{15}\,\rm GeV$, respectively. The cut-off
represent the corresponding minimum PBH mass $M_{\rm PBH}^{\rm min}$.

In our scenario, the PBHs will be formed when the axions dominate the
radiation in the bubbles, $\rm i.e.$,

$\displaystyle T_{f}\lesssim T_{B}\,.$ (39)

The condition for PBHs formation is that when the bubbles enter the horizontal
plane, the local energy density inside the bubbles is significantly greater
than the background radiation density. Substituting Eq. 29 in Eq. 37, then we
have the minimum PBH mass

$\displaystyle\begin{aligned} \frac{M_{\rm PBH}^{\rm min}}{M_{\odot}}&\simeq
5.65\times
10^{5}\left(\frac{\gamma}{0.2}\right)\left(\frac{g_{*}(T_{f})}{10.75}\right)^{-1/2}\\\
&\times\left(\frac{g_{*}(T_{a})}{61.75}\right)^{0.84}\left(\frac{f_{a}}{10^{16}\,\rm
GeV}\right)^{-2.33}\,.\end{aligned}$ (40)

The PBH mass at this region may account for the LIGO-Virgo gravitational waves
(GWs) events Abbott  _et al._ (2016a, b, 2017a, 2017b) and the seeds of
supermassive black holes (SMBHs) Richstone  _et al._ (1998); Ferrarese and
Merritt (2000); Kormendy and Ho (2013). The PBH energy density at $T_{B}$ is
given by

$\displaystyle\rho_{\rm PBH}(T_{B})=\frac{3}{4}T_{B}s(T_{B})\frac{{\rm
d}\beta}{{\rm d}\ln k}\,,$ (41)

where $s(T_{B})$ is the entropy density. Finally, we derive the fractional
abundance of PBH at present

$\displaystyle f_{\rm PBH}=\frac{\Omega_{\rm PBH}}{\Omega_{\rm
DM}}=\frac{3}{4}T_{B}s(T_{0})\frac{{\rm d}\beta}{{\rm d}\ln
k}\frac{1}{\rho_{c}}\frac{1}{\Omega_{\rm DM}}\,,$ (42)

with the total cold DM abundance $\Omega_{\rm DM}\simeq 0.268$ Aghanim  _et
al._ (2020). With Eq. 35, we show the PBH fractional abundance as a function
of the PBH mass $M_{\rm PBH}$ in Fig. 2. The three typical values of
$f_{a}=5\times 10^{16}\,\rm GeV$ (red), $1\times 10^{16}\,\rm GeV$ (orange),
and $5\times 10^{15}\,\rm GeV$ (blue) are selected for comparisons. The cut-
off represent the minimum PBH mass $M_{\rm PBH}^{\rm min}$ (Eq. 40). Here we
take $\phi_{\rm crit}-\phi_{i}=4.5H_{\rm inf}$ as a typical value. The
fractional abundance would be higher for the large scale
$f_{a}\sim\mathcal{O}(10^{17})\,\rm GeV$.

Additionally, the another interesting phenomenon is the formation of axion
miniclusters from the QCD axion bubbles. The axion miniclusters Hogan and Rees
(1988); Kolb and Tkachev (1993); Fairbairn  _et al._ (2018); Ellis  _et al._
(2022); Dandoy  _et al._ (2022) are gravitationally bound clumps of the axion
DM. Their mass and size depend on the Hubble volume when the QCD axion starts
to oscillate. In our scenario, the axion miniclusters can be formed when the
bubbles enter the horizontal plane before the axions dominate the radiation in
the bubbles, $\rm i.e.$,

$\displaystyle T_{f}>T_{B}\,,$ (43)

where $T_{f}$ is the temperature at the miniclusters formation time. The
miniclusters mass is related to the QCD axion oscillation temperature $T_{a}$
(Eq. 5). Since we are more concerned about the phenomena of PBHs formation,
the miniclusters will not be further discussed in this context, and more
details about the axion miniclusters from the QCD axion bubbles can be found
in Ref. Kitajima and Takahashi (2020).

##  V Discussion and conclusion

In this paper, we have presented the QCD axion bubbles formation from the
dynamically broken PQ symmetry in the early Universe. First we introduce the
QCD axion and the misalignment mechanism. Then we consider that the PQ
symmetry is broken by the hidden $SU(N)_{H}$ gauge interactions, which is
assumed to be broken after inflation with the multiple vacua. The axion
acquires a light effective mass and then settles down into different vacua.
Note that this symmetry is temporarily broken, and the resulting axion
potential will be disappeared before the QCD axion arises. During the QCD
phase transition, the QCD axion starts to oscillate near the origin and the
another large value. In our scenario, we consider the two cases with $N=2$ and
3. We find that the QCD axions start to oscillate at the large values $\sim
2\pi/3$ can lead to the high density QCD axion bubbles with $N=2$. However,
since the effective initial angle $\theta_{i,0}=\theta_{i,1}\simeq 0$, we find
no axion bubbles is generated with $N=3$.

We also discuss the cosmological implications of the QCD axion bubbles. Since
the large axion scale $f_{a}\sim\mathcal{O}(10^{16})\,\rm GeV$, the initial
misalignment angle should be small $\sim\mathcal{O}(10^{-3})$ to explain the
cold DM abundance. The interesting phenomenon in our scenario is the PBHs
formation. The PBHs will be formed when the axions dominate the radiation in
the bubbles. We find that the PBH mass is lager than $\sim\mathcal{O}(5\times
10^{5})M_{\odot}$ for the scale $f_{a}\sim\mathcal{O}(10^{16})\,\rm GeV$. The
PBH mass at this region may explain the LIGO-Virgo GWs events and the seeds of
SMBHs. In addition, if the bubbles enter the horizontal plane before the
axions dominate the radiation, the bubbles will eventually cluster into the
axion miniclusters.

## Acknowledgments

The author would like to thank Wei Chao, Naoya Kitajima, Shota Nakagawa, and
Fuminobu Takahashi for helpful discussions and valuable comments. This work
was supported by the National Natural Science Foundation (NNSF) of China
(Grants No. 11775025 and No. 12175027).

## References

  * Peccei and Quinn (1977a) R.D. Peccei and Helen R. Quinn, “Constraints Imposed by CP Conservation in the Presence of Instantons,” Phys. Rev. D 16, 1791–1797 (1977a).
  * Peccei and Quinn (1977b) R.D. Peccei and Helen R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38, 1440–1443 (1977b).
  * Weinberg (1978) Steven Weinberg, “A New Light Boson?” Phys. Rev. Lett. 40, 223–226 (1978).
  * Wilczek (1978) Frank Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40, 279–282 (1978).
  * ’t Hooft (1976a) Gerard ’t Hooft, “Symmetry Breaking Through Bell-Jackiw Anomalies,” Phys. Rev. Lett. 37, 8–11 (1976a).
  * ’t Hooft (1976b) Gerard ’t Hooft, “Computation of the Quantum Effects Due to a Four-Dimensional Pseudoparticle,” Phys. Rev. D 14, 3432–3450 (1976b), [Erratum: Phys.Rev.D 18, 2199 (1978)].
  * Preskill  _et al._ (1983) John Preskill, Mark B. Wise, and Frank Wilczek, “Cosmology of the Invisible Axion,” Phys. Lett. B 120, 127–132 (1983).
  * Abbott and Sikivie (1983) L.F. Abbott and P. Sikivie, “A Cosmological Bound on the Invisible Axion,” Phys. Lett. B 120, 133–136 (1983).
  * Dine and Fischler (1983) Michael Dine and Willy Fischler, “The Not So Harmless Axion,” Phys. Lett. B 120, 137–141 (1983).
  * Di Luzio  _et al._ (2020) Luca Di Luzio, Maurizio Giannotti, Enrico Nardi, and Luca Visinelli, “The landscape of QCD axion models,” Phys. Rept. 870, 1–117 (2020), arXiv:2003.01100 [hep-ph] .
  * Chadha-Day  _et al._ (2022) Francesca Chadha-Day, John Ellis, and David J. E. Marsh, “Axion dark matter: What is it and why now?” Sci. Adv. 8, abj3618 (2022), arXiv:2105.01406 [hep-ph] .
  * Adams  _et al._ (2022) C. B. Adams _et al._ , “Axion Dark Matter,” in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.14923 [hep-ex] .
  * Dine and Anisimov (2005) Michael Dine and Alexey Anisimov, “Is there a Peccei-Quinn phase transition?” JCAP 07, 009 (2005), arXiv:hep-ph/0405256 .
  * Jeong and Takahashi (2013) Kwang Sik Jeong and Fuminobu Takahashi, “Suppressing Isocurvature Perturbations of QCD Axion Dark Matter,” Phys. Lett. B 727, 448–451 (2013), arXiv:1304.8131 [hep-ph] .
  * Takahashi and Yamada (2015) Fuminobu Takahashi and Masaki Yamada, “Strongly broken Peccei-Quinn symmetry in the early Universe,” JCAP 10, 010 (2015), arXiv:1507.06387 [hep-ph] .
  * Harigaya  _et al._ (2015) Keisuke Harigaya, Masahiro Ibe, Kai Schmitz, and Tsutomu T. Yanagida, “Peccei-Quinn Symmetry from Dynamical Supersymmetry Breaking,” Phys. Rev. D 92, 075003 (2015), arXiv:1505.07388 [hep-ph] .
  * Kallosh  _et al._ (1995) Renata Kallosh, Andrei D. Linde, Dmitri A. Linde, and Leonard Susskind, “Gravity and global symmetries,” Phys. Rev. D 52, 912–935 (1995), arXiv:hep-th/9502069 .
  * Banks and Seiberg (2011) Tom Banks and Nathan Seiberg, “Symmetries and Strings in Field Theory and Gravity,” Phys. Rev. D 83, 084019 (2011), arXiv:1011.5120 [hep-th] .
  * Witten (2018) Edward Witten, “Symmetry and Emergence,” Nature Phys. 14, 116–119 (2018), arXiv:1710.01791 [hep-th] .
  * Harlow and Ooguri (2019) Daniel Harlow and Hirosi Ooguri, “Constraints on Symmetries from Holography,” Phys. Rev. Lett. 122, 191601 (2019), arXiv:1810.05337 [hep-th] .
  * Jeong  _et al._ (2022) Kwang Sik Jeong, Kohei Matsukawa, Shota Nakagawa, and Fuminobu Takahashi, “Cosmological effects of Peccei-Quinn symmetry breaking on QCD axion dark matter,” JCAP 03, 026 (2022), arXiv:2201.00681 [hep-ph] .
  * Kitajima and Takahashi (2020) Naoya Kitajima and Fuminobu Takahashi, “Primordial Black Holes from QCD Axion Bubbles,” JCAP 11, 060 (2020), arXiv:2006.13137 [hep-ph] .
  * Dolgov  _et al._ (2009) A. D. Dolgov, M. Kawasaki, and N. Kevlishvili, “Inhomogeneous baryogenesis, cosmic antimatter, and dark matter,” Nucl. Phys. B 807, 229–250 (2009), arXiv:0806.2986 [hep-ph] .
  * Hasegawa and Kawasaki (2019) Fuminori Hasegawa and Masahiro Kawasaki, “Primordial Black Holes from Affleck-Dine Mechanism,” JCAP 01, 027 (2019), arXiv:1807.00463 [astro-ph.CO] .
  * Borsanyi  _et al._ (2016) Sz. Borsanyi _et al._ , “Calculation of the axion mass based on high-temperature lattice quantum chromodynamics,” Nature 539, 69–71 (2016), arXiv:1606.07494 [hep-lat] .
  * Grilli di Cortona  _et al._ (2016) Giovanni Grilli di Cortona, Edward Hardy, Javier Pardo Vega, and Giovanni Villadoro, “The QCD axion, precisely,” JHEP 01, 034 (2016), arXiv:1511.02867 [hep-ph] .
  * Turner (1986) Michael S. Turner, “Cosmic and Local Mass Density of Invisible Axions,” Phys. Rev. D 33, 889–896 (1986).
  * Lyth (1992) D. H. Lyth, “Axions and inflation: Sitting in the vacuum,” Phys. Rev. D 45, 3394–3404 (1992).
  * Visinelli and Gondolo (2009) Luca Visinelli and Paolo Gondolo, “Dark Matter Axions Revisited,” Phys. Rev. D 80, 035024 (2009), arXiv:0903.4377 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Nomura  _et al._ (2016) Yasunori Nomura, Surjeet Rajendran, and Fabio Sanches, “Axion Isocurvature and Magnetic Monopoles,” Phys. Rev. Lett. 116, 141803 (2016), arXiv:1511.06347 [hep-ph] .
  * Kawasaki  _et al._ (2016) Masahiro Kawasaki, Fuminobu Takahashi, and Masaki Yamada, “Suppressing the QCD Axion Abundance by Hidden Monopoles,” Phys. Lett. B 753, 677–681 (2016), arXiv:1511.05030 [hep-ph] .
  * Kawasaki  _et al._ (2018) Masahiro Kawasaki, Fuminobu Takahashi, and Masaki Yamada, “Adiabatic suppression of the axion abundance and isocurvature due to coupling to hidden monopoles,” JHEP 01, 053 (2018), arXiv:1708.06047 [hep-ph] .
  * Chiba  _et al._ (2004) Takeshi Chiba, Fuminobu Takahashi, and Masahide Yamaguchi, “Baryogenesis in a flat direction with neither baryon nor lepton charge,” Phys. Rev. Lett. 92, 011301 (2004), [Erratum: Phys.Rev.Lett. 114, 209901 (2015)], arXiv:hep-ph/0304102 .
  * Takahashi and Yamaguchi (2004) Fuminobu Takahashi and Masahide Yamaguchi, “Spontaneous baryogenesis in flat directions,” Phys. Rev. D 69, 083506 (2004), arXiv:hep-ph/0308173 .
  * Higaki  _et al._ (2014) Tetsutaro Higaki, Kwang Sik Jeong, and Fuminobu Takahashi, “Solving the Tension between High-Scale Inflation and Axion Isocurvature Perturbations,” Phys. Lett. B 734, 21–26 (2014), arXiv:1403.4186 [hep-ph] .
  * Co  _et al._ (2020) Raymond T. Co, Lawrence J. Hall, and Keisuke Harigaya, “Axion Kinetic Misalignment Mechanism,” Phys. Rev. Lett. 124, 251802 (2020), arXiv:1910.14152 [hep-ph] .
  * Choi  _et al._ (1997) Kiwoon Choi, Hang Bae Kim, and Jihn E. Kim, “Axion cosmology with a stronger QCD in the early universe,” Nucl. Phys. B 490, 349–364 (1997), arXiv:hep-ph/9606372 .
  * Banks and Dine (1997) Tom Banks and Michael Dine, “The Cosmology of string theoretic axions,” Nucl. Phys. B 505, 445–460 (1997), arXiv:hep-th/9608197 .
  * Starobinsky (1982) Alexei A. Starobinsky, “Dynamics of Phase Transition in the New Inflationary Universe Scenario and Generation of Perturbations,” Phys. Lett. B 117, 175–178 (1982).
  * Linde (1982) Andrei D. Linde, “Scalar Field Fluctuations in Expanding Universe and the New Inflationary Universe Scenario,” Phys. Lett. B 116, 335–339 (1982).
  * Starobinsky and Yokoyama (1994) Alexei A. Starobinsky and Junichi Yokoyama, “Equilibrium state of a selfinteracting scalar field in the De Sitter background,” Phys. Rev. D 50, 6357–6368 (1994), arXiv:astro-ph/9407016 .
  * Bird  _et al._ (2016) Simeon Bird, Ilias Cholis, Julian B. Muñoz, Yacine Ali-Haïmoud, Marc Kamionkowski, Ely D. Kovetz, Alvise Raccanelli, and Adam G. Riess, “Did LIGO detect dark matter?” Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
  * Carr and Kuhnel (2020) Bernard Carr and Florian Kuhnel, “Primordial Black Holes as Dark Matter: Recent Developments,” Ann. Rev. Nucl. Part. Sci. 70, 355–394 (2020), arXiv:2006.02838 [astro-ph.CO] .
  * Green and Kavanagh (2021) Anne M. Green and Bradley J. Kavanagh, “Primordial Black Holes as a dark matter candidate,” J. Phys. G 48, 043001 (2021), arXiv:2007.10722 [astro-ph.CO] .
  * Carr and Kuhnel (2022) Bernard Carr and Florian Kuhnel, “Primordial black holes as dark matter candidates,” SciPost Phys. Lect. Notes 48, 1 (2022), arXiv:2110.02821 [astro-ph.CO] .
  * Li (2022) Hai-Jun Li, “Primordial black holes induced stochastic axion-photon oscillations in primordial magnetic field,” JCAP 11, 045 (2022), arXiv:2208.04605 [astro-ph.CO] .
  * Carr  _et al._ (2010) B. J. Carr, Kazunori Kohri, Yuuiti Sendouda, and Jun’ichi Yokoyama, “New cosmological constraints on primordial black holes,” Phys. Rev. D 81, 104019 (2010), arXiv:0912.5297 [astro-ph.CO] .
  * Carr (1975) Bernard J. Carr, “The Primordial black hole mass spectrum,” Astrophys. J. 201, 1–19 (1975).
  * Di and Gong (2018) Haoran Di and Yungui Gong, “Primordial black holes and second order gravitational waves from ultra-slow-roll inflation,” JCAP 07, 007 (2018), arXiv:1707.09578 [astro-ph.CO] .
  * Abbott  _et al._ (2016a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence,” Phys. Rev. Lett. 116, 241103 (2016b), arXiv:1606.04855 [gr-qc] .
  * Abbott  _et al._ (2017a) Benjamin P. Abbott _et al._ (LIGO Scientific, VIRGO), “GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2,” Phys. Rev. Lett. 118, 221101 (2017a), [Erratum: Phys.Rev.Lett. 121, 129901 (2018)], arXiv:1706.01812 [gr-qc] .
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “GW170814: A Three-Detector Observation of Gravitational Waves from a Binary Black Hole Coalescence,” Phys. Rev. Lett. 119, 141101 (2017b), arXiv:1709.09660 [gr-qc] .
  * Richstone  _et al._ (1998) D. Richstone _et al._ , “Supermassive black holes and the evolution of galaxies,” Nature 395, A14–A19 (1998), arXiv:astro-ph/9810378 .
  * Ferrarese and Merritt (2000) Laura Ferrarese and David Merritt, “A Fundamental relation between supermassive black holes and their host galaxies,” Astrophys. J. Lett. 539, L9 (2000), arXiv:astro-ph/0006053 .
  * Kormendy and Ho (2013) John Kormendy and Luis C. Ho, “Coevolution (Or Not) of Supermassive Black Holes and Host Galaxies,” Ann. Rev. Astron. Astrophys. 51, 511–653 (2013), arXiv:1304.7762 [astro-ph.CO] .
  * Hogan and Rees (1988) C. J. Hogan and M. J. Rees, “AXION MINICLUSTERS,” Phys. Lett. B 205, 228–230 (1988).
  * Kolb and Tkachev (1993) Edward W. Kolb and Igor I. Tkachev, “Axion miniclusters and Bose stars,” Phys. Rev. Lett. 71, 3051–3054 (1993), arXiv:hep-ph/9303313 .
  * Fairbairn  _et al._ (2018) Malcolm Fairbairn, David J. E. Marsh, Jérémie Quevillon, and Simon Rozier, “Structure formation and microlensing with axion miniclusters,” Phys. Rev. D 97, 083502 (2018), arXiv:1707.03310 [astro-ph.CO] .
  * Ellis  _et al._ (2022) David Ellis, David J. E. Marsh, Benedikt Eggemeier, Jens Niemeyer, Javier Redondo, and Klaus Dolag, “Structure of axion miniclusters,” Phys. Rev. D 106, 103514 (2022), arXiv:2204.13187 [hep-ph] .
  * Dandoy  _et al._ (2022) Virgile Dandoy, Thomas Schwetz, and Elisa Todarello, “A self-consistent wave description of axion miniclusters and their survival in the galaxy,” JCAP 09, 081 (2022), arXiv:2206.04619 [astro-ph.CO] .


# Accessible Computation of Tight Symbolic Bounds on Causal Effects using an
Intuitive Graphical Interface

by Gustav Jonzon, Michael C Sachs, and Erin E Gabriel

###### Abstract

Strong untestable assumptions are almost universal in causal point estimation.
In particular settings, bounds can be derived to narrow the possible range of
a causal effect. Symbolic bounds apply to all settings that can be depicted
using the same directed acyclic graph (DAG) and for the same effect of
interest. Although the core of the methodology for deriving symbolic bounds
has been previously developed, the means of implementation and computation
have been lacking. Our R-package causaloptim (Sachs et al., 2022b) aims to
solve this usability problem by implementing the method of (Sachs et al.,
2022a) and providing the user with a graphical interface through shiny that
allows for input in a way that most researchers with an interest in causal
inference will be familiar; a DAG (via a point-and-click experience) and
specifying a causal effect of interest using familiar counterfactual notation.

## Introduction and Background

A common goal in many different areas of scientific research is to determine
causal relationships between one or more exposure variables and an outcome.
Prior to any computation or inference, we must clearly state all assumptions
made, i.e., all subject matter knowledge available, regarding the causal
relationships between the involved variables as well as any additional
variables, called confounders, that may not be measured but influence at least
two other variables of interest. These assumptions are usually encoded in a
causal directed acyclic graph (DAG), with directed edges encoding direct
causal influences, which conveniently depicts all relevant information and has
become a familiar tool in applied research (Greenland et al., 1999). Such a
DAG not only clearly states the assumptions made by the researcher, but also
comes with a sound methodology for causal inference, in the form of
identification results as well as derivation of causal estimators (Pearl,
2009).

Unfortunately, point identification of a desired causal effect typically
requires an assumption of no unmeasured confounders, in some form. When there
are unmeasured confounders, it is sometimes still possible to derive bounds on
the effect, i.e., a range of possible values for the causal effect in terms of
the observed data distribution. Symbolic bounds are algebraic expressions for
the bounds on the causal effect written in terms of probabilities that can be
estimated using observed data. Alexander Balke and Judea Pearl first used
linear programming to derive tight symbolic bounds in a simple binary
instrumental variable setting (Balke and Pearl, 1997). Balke wrote a program
in C++ to take a linear programming problem as text file input, perform
variable reduction, conversion of equality constraints into inequality
constraints, and perform the vertex enumeration algorithm of (Mattheiss,
1973). This program has been used by researchers in the field of causal
inference (Balke and Pearl, 1997; Cai et al., 2008; Sjölander, 2009; Sjölander
et al., 2014) but it is not particularly accessible because of the technical
challenge of translating the DAG plus causal query into the constrained
optimization problem and to determine whether it is linear. Moreover, the
program is not optimized and hence does not scale well to more complex
problems. Since they only cover a simple instrumental variable setting, it has
also not been clear to what extent their techniques extend to more general
settings, nor how to apply them to more complex queries. Thus, applications of
this approach have been limited to a small number of settings and few attempts
to generalize the method to more widely applicable settings have been made.

Recent developments have expanded the applicability by generalizing the
techniques and the causal DAGs and effects to which they apply (Sachs et al.,
2022a). These new methods have been applied in novel observational and
experimental settings (Gabriel et al., 2020, 2021, 2022). Moreover, through
the R package causaloptim (Sachs et al., 2022b), these computations are now
accessible. With causaloptim, the user needs only to give input in a way they
would usually express their causal assumptions and state their target causal
estimand; through a DAG and counterfactual expression. Providing DAGs through
textual input is an awkward experience for most users, as DAGs are generally
communicated pictorially. Our package causaloptim provides a user-friendly
graphical interface through a web browser, where the user can draw their DAG
in a way that is familiar to them. The methodology that underpins causaloptim
is not universal however; some restrictions on the DAG and query are imposed.
These are validated and communicated to the user through the graphical
interface, which guides the user through providing the DAG and query, adding
any extra conditions beyond those encoded in the DAG, computing, interpreting
and exporting the bounds for various further analyses.

There exist few other R-packages related to causal bounds and none to our
knowledge for computation of symbolic bounds. bpbounds (Palmer et al., 2018)
provides a text-based interface to compute numeric bounds for the original
single instrumental variable example of Balke and Pearl and extends this by
being able to compute bounds given different types of data input including a
ternary rather than binary instrument. There is also a standalone program
written in Java by the TETRAD Project (https://github.com/cmu-phil/tetrad)
that includes a graphical user interface and has a wrapper for R. Its focus,
however, is on causal discovery in a given sample data set, and although it
can also compute bounds, it can do so only numerically for the given data set.

In this paper we describe our R package causaloptim, first focusing on the
graphical and programmatic user interfaces in the next 2 sections. Then we
highlight some of our interesting functions and data structures that may be
useful in other contexts. We provide a summary of the theoretical background
and methods, while referring to the companion paper (Sachs et al., 2022a) for
the details. We illustrate the use of the package with some numeric examples
and close with a discussion and summary.

## Graphical User Interface

In the following, we will work through the binary instrumental variable
example, where we have 3 observed binary variables $X$, $Y$, $Z$, and we want
to determine the average causal effect of $X$ on $Y$ given by total causal
risk difference, in the presence of unmeasured confounding by $U_{R}$ and an
instrumental variable $Z$. Our causal DAG is given by $Z\to X\to Y$ and
$X\leftarrow U_{R}\to Y$ and our causal query is $P(Y(X=1)=1)-P(Y(X=0)=1)$,
where we use $Y(X=x)$ to denote the potential outcome for $Y$ if $X$ were
intervened upon to have value $x$.

causaloptim includes a graphical user interface using shiny (Chang et al.,
2021). The interface is launched in the user’s default web browser by calling
specify_graph(). Once the shiny app is launched, the user is presented with an
interactive display as shown in Figure 1, in which they can draw their causal
DAG. This display is divided into a left side $\mathcal{L}$ and right side
$\mathcal{R}$ to classify the vertices according to the class of DAGs that the
method covers. In particular, the existence of unmeasured confounders is
assumed within each of these sides, but not between them, and any causal
influence between the two sides must originate in $\mathcal{L}$. Thus, for the
example, we would want to put the instrumental variable on the left side, but
the exposure and outcome on the right side. In the web version of this article
an interactive version of this interface is shown at the end of this section.

Figure 1: The Shiny web interface at lauch

### Specifying the setting by drawing a causal diagram and adding attributes

(a) Adding and naming variables

(b) Adding directed edges

Figure 2: Constructing the DAG

The DAG is drawn using a point-and-click device (e.g., a mouse) to add
vertices representing variables (by Shift-click) and name them (using any
valid variable name in R), and to draw edges representing direct causal
influences (Shift+drag) between them. The vertices may also be moved around,
renamed and deleted (as can the edges) as also described in an instruction
text preceding the DAG interface. As shown in Figure 2, for the example we add
a vertex $Z$ on the left side, and vertices $X$ and $Y$ on the right side.
Then the $Z\to X$ and $X\to Y$ edges are added by Shift+clicking on a parent
vertex and dragging to the child vertex. There is no need to add the
unmeasured confounder variable $U_{R}$ as it is assumed and added
automatically.

Importantly, the nodes may be selected and assigned additional information. In
$\mathcal{R}$ a variable may be assigned as unobserved (click+‘u’). All
observed variables are assumed categorical and their cardinality (i.e., number
of levels) may be set (click+‘c’ brings up a prompt for this this number;
alternatively a short-cut click+‘any digit’ is provided), with the default
being binary. Although the causal query (i.e., the causal effect of interest)
is entered subsequently, the DAG interface provides a convenient short-cut; a
node $X$ may be assigned as an exposure (click+‘e’) and another $Y$ as outcome
(click+‘y’), whereupon the default query is the total causal risk difference
$P(Y(X=1)=1)-P(Y(X=0)=1)$. Finally, an edge may be assigned as representing an
assumed monotonic influence (click+‘m’). The nodes and edges change appearance
according to their assigned characteristics (Figure 3) and violations to the
restrictions characterizing the class of DAGs are detected and communicated to
the user.

(a) Setting the number of categories

(b) Confirmation message

(c) Setting exposure and outcome

Figure 3: Setting attributes

Once the DAG has been drawn, the user may click the button “Analyze the
graph”, upon which the DAG is interpreted and converted into an annotated
igraph-object (Csardi and Nepusz, 2006) as described in the implementation
details below, and the results are displayed in graphical form to the user
(Figure 4). The addition of $U_{R}$, the common unmeasured cause of $X$ and
$Y$, is added and displayed in this static plot.

(a) Graphical summary of the DAG with added confounding

(b) Computing the bounds

Figure 4: The causal DAG and bounds

### Specifying the causal query

Next, the user is asked to specify the causal query, i.e., causal effect of
interest. If no outcome variable has been assigned in the DAG then the input
field for the causal query is left blank and a query needs to be specified. In
our example, since we have assigned an exposure and outcome using the DAG
interface, the total causal risk difference $P(Y(X=1)=1)-P(Y(X=0)=1)$ is
suggested.

### Specifying optional additional constraints

Finally, the user is given the option to provide any additional constraints
besides those imposed by the DAG. This may be considered an optional advanced
feature where, e.g., monotonicity of a certain direct influence of $Z$ on $X$
may be assumed by entering $X(Z=1)\geq X(Z=0)$, with any such extra
constraints separated by line breaks. If this feature is used, the input is
followed by clicking the button “Parse”, which identifies and fixes them.

### Computing the symbolic tight bounds on the query under the given
constraints

As the final step, the button “Compute the bounds” is clicked, whereupon the
constraints and objective are compiled into an optimization problem which is
then solved for tight causal bounds on the query symbolically in terms of
observational quantities (conditional probabilities of the observed variables
in the DAG) and the expressions are displayed alongside information on how the
parameters are to be interpreted in terms of the given variable names (Figure
4). During computation, a progress indicator is shown, and the user should be
aware that complex and/or high-dimensional problems may take significant time.
The interface also provides a feature to subsequently convert the bounds to
LaTeX-code using standard probabilistic notation for publication purposes.

Once done, clicking “Exit and return objects to R” stops the shiny app and
returns all information about the DAG, query and computed bounds to the
R-session. This information is bundled in a list containing the graph, query,
parameters and their interpretation, the symbolic tight bounds as expressions
as well as implementations as R-functions and further log information about
the formulation and optimization procedures.

## Programmatic user interface

Interaction may also be done entirely programmatically as we illustrate with
the same binary instrumental variable example. First we create the igraph
object using the graph_from_literal function. Once the basic graph is created,
the necessary vertex and edge attributes are added. The risk difference is
defined as a character object. The analyze_graph function is the workhorse of
causaloptim; it translates the causal graph, constraints, and causal effect of
interest into a linear programming problem. This linear programming object,
stored in obj in the code below, gets passed to optimize_effect_2 which
performs vertex enumeration to obtain the bounds as symbolic expressions in
terms of observable probabilities.

graph <\- igraph::graph_from_literal(Z -+ X, X -+ Y, Ul -+ Z, Ur -+ X, Ur -+
Y)V(graph)$leftside <\- c(1, 0, 0, 1, 0)V(graph)$latent <\- c(0, 0, 0, 1,
1)V(graph)$nvals <\- c(2, 2, 2, 2, 2)E(graph)$rlconnect <\- c(0, 0, 0, 0,
0)E(graph)$edge.monotone <\- c(0, 0, 0, 0, 0)riskdiff <\- "p{Y(X = 1) = 1} -
p{Y(X = 0) = 1}"obj <\- analyze_graph(graph, constraints = NULL, effectt =
riskdiff)bounds <\- optimize_effect_2(obj)bounds

#> lower bound =#> MAX {#> p00_0 - p00_1 - p10_1 - p01_1,#> p00_0 - p00_1 -
p10_0 - p10_1 - p01_0,#> p00_0 - p00_1 + p10_0 - 2p10_1 - 2p01_1,#> -p10_1 -
p01_1,#> -p10_0 - p01_0,#> -p00_0 + p00_1 - 2p10_0 + p10_1 - 2p01_0,#> -p00_0
+ p00_1 - p10_0 - p10_1 - p01_1,#> -p00_0 + p00_1 - p10_0 - p01_0#> }#>
\----------------------------------------#> upper bound =#> MIN {#> 1 - p10_1
- p01_0,#> 1 + p00_0 + p10_0 - 2p10_1 - p01_1,#> 2 - p00_1 - p10_0 - p10_1 -
2p01_0,#> 1 - p10_1 - p01_1,#> 1 - p10_0 - p01_0,#> 1 + p00_1 - 2p10_0 + p10_1
- p01_0,#> 2 - p00_0 - p10_0 - p10_1 - 2p01_1,#> 1 - p10_0 - p01_1#> }

The resulting bounds object contains character strings representing the bounds
and logs containing detailed information from the vertex enumeration
algorithm. The bounds are printed to the console but more features are
available to facilitate their use. The interpret_bounds function takes the
bounds and parameter names as input and returns an R function implementing
vectorized forms of the symbolic expressions for the bounds.

bounds_function <\- interpret_bounds(bounds$bounds,
obj$parameters)str(bounds_function)

#> function (p00_0 = NULL, p00_1 = NULL, p10_0 = NULL, p10_1 = NULL, p01_0 =
NULL,#> p01_1 = NULL, p11_0 = NULL, p11_1 = NULL)

The results can also be used for numerical simulation using simulate_bounds.
This function randomly generates counterfactuals and probability distributions
that satisfy the constraints implied by the DAG and optional constraints. It
then computes and returns the bounds as well as the true causal effect.

If one wants to bound a different effect using the same causal graph, the
update_effect function can be used to save some computation time. It takes the
object returned by analyze_graph and the new effect string then returns an
object of class linearcausalproblem that can be optimized: obj2 <\-
update_effect(obj, "p{Y(X = 1) = 1}").

Finally, LaTeX-code may also be generated using the function latex_bounds as
in latex_bounds(bounds$bounds, obj$parameters) yielding

Lower bound
$\displaystyle=\mbox{max}\left.\begin{cases}P(X=0,Y=0|Z=0)-P(X=0,Y=0|Z=1)-P(X=1,Y=0|Z=1)-P(X=0,Y=1|Z=1),\\\
P(X=0,Y=0|Z=0)-P(X=0,Y=0|Z=1)-P(X=1,Y=0|Z=0)-P(X=1,Y=0|Z=1)-P(X=0,Y=1|Z=0),\\\
P(X=0,Y=0|Z=0)-P(X=0,Y=0|Z=1)+P(X=1,Y=0|Z=0)-2P(X=1,Y=0|Z=1)-2P(X=0,Y=1|Z=1),\\\
-P(X=1,Y=0|Z=1)-P(X=0,Y=1|Z=1),\\\ -P(X=1,Y=0|Z=0)-P(X=0,Y=1|Z=0),\\\
-P(X=0,Y=0|Z=0)+P(X=0,Y=0|Z=1)-2P(X=1,Y=0|Z=0)+P(X=1,Y=0|Z=1)-2P(X=0,Y=1|Z=0),\\\
-P(X=0,Y=0|Z=0)+P(X=0,Y=0|Z=1)-P(X=1,Y=0|Z=0)-P(X=1,Y=0|Z=1)-P(X=0,Y=1|Z=1),\\\
-P(X=0,Y=0|Z=0)+P(X=0,Y=0|Z=1)-P(X=1,Y=0|Z=0)-P(X=0,Y=1|Z=0)\end{cases}\right\\}$
Upper bound
$\displaystyle=\mbox{min}\left.\begin{cases}1-P(X=1,Y=0|Z=1)-P(X=0,Y=1|Z=0),\\\
1+P(X=0,Y=0|Z=0)+P(X=1,Y=0|Z=0)-2P(X=1,Y=0|Z=1)-P(X=0,Y=1|Z=1),\\\
2-P(X=0,Y=0|Z=1)-P(X=1,Y=0|Z=0)-P(X=1,Y=0|Z=1)-2P(X=0,Y=1|Z=0),\\\
1-P(X=1,Y=0|Z=1)-P(X=0,Y=1|Z=1),\\\ 1-P(X=1,Y=0|Z=0)-P(X=0,Y=1|Z=0),\\\
1+P(X=0,Y=0|Z=1)-2P(X=1,Y=0|Z=0)+P(X=1,Y=0|Z=1)-P(X=0,Y=1|Z=0),\\\
2-P(X=0,Y=0|Z=0)-P(X=1,Y=0|Z=0)-P(X=1,Y=0|Z=1)-2P(X=0,Y=1|Z=1),\\\
1-P(X=1,Y=0|Z=0)-P(X=0,Y=1|Z=1)\end{cases}\right\\}.$

## Implementation and Program Overview

An overview of the main functions and their relations is depicted as a flow
chart in Figure 5. All functions may be called individually by the user at the
R-console and all input, output and interaction available through the shiny
app has corresponding availability at the R-console as well.

(Sachs et al., 2022a) define the following class of problems for which the
query in general is not identifiable, but for which a methodology to derive
symbolic tight bounds on the query is provided. The causal DAG consists of a
finite set
$\mathcal{W}=\\{W_{1},\dots,W_{n}\\}=\mathcal{W}_{\mathcal{L}}\cup\mathcal{W}_{\mathcal{R}}$
of categorical variables with
$\mathcal{W}_{\mathcal{L}}\cap\mathcal{W}_{\mathcal{R}}=\varnothing$, no edges
going from $\mathcal{W}_{\mathcal{R}}$ to $\mathcal{W}_{\mathcal{L}}$ and no
external common parent between $\mathcal{W}_{\mathcal{L}}$ and
$\mathcal{W}_{\mathcal{R}}$, but _importantly_ external common parents
$\mathcal{U}_{\mathcal{L}}$ and $\mathcal{U}_{\mathcal{R}}$ of variables
within $\mathcal{W}_{\mathcal{R}}$ and $\mathcal{W}_{\mathcal{R}}$ may not be
ruled out. Nothing is assumed about any characteristics of these confounding
variables $\mathcal{U}_{\mathcal{L}}$ and $\mathcal{U}_{\mathcal{R}}$.

The causal query may be any linear combination of joint probabilities of
factual and counterfactual outcomes expressed in terms of the variables in
$\mathcal{W}$ and may always be expressed as a sum of probabilities of
response function variables of the DAG. It is subject to the restriction that
each outcome variable is in $\mathcal{W}_{\mathcal{R}}$ and if
$\mathcal{W}_{\mathcal{L}}\neq\varnothing$ it is also subject to a few
regularity conditions as detailed in (Sachs et al., 2022a). Tight bounds on
the query may then be derived symbolically in terms of conditional
probabilities of the observable variables $\mathcal{W}$.

Algorithms 1 and 2 in (Sachs et al., 2022a) construct the constraint space and
causal query in terms of the joint probabilities of the response function
variables and in causaloptim are implemented in the functions create_R_matrix
and create_effect_vector respectively. Both are called as sub-procedures of
the function analyze_graph to translate the causal problem to that of
optimizing an objective function over a constraint space. The implementation
of Algorithm 1 involves constructing the response functions themselves as
actual R-functions. Evaluating these correspond to evaluating the structural
equations of the causal DAG.

The conditions on the DAG suffice to ensure that the causal query will depend
only on the response functions corresponding to the variables in
$\mathcal{W}_{\mathcal{R}}$ and that the exhaustive set of constraints of
their probabilities are linear in a subset of conditional probabilities of
observable variables (Proposition 2 in (Sachs et al., 2022a)), and the
conditions on the query in turn ensure that it may be expressed as a linear
combination of joint probabilities of the response functions of the variables
in $\mathcal{W}_{\mathcal{R}}$ (Proposition 3 in (Sachs et al., 2022a)).

Once this formulation of the causal problem as a linear program has been set
up, a vertex enumeration method is employed to compute the extrema
symbolically in terms of conditional probabilities of the observable
variables.

The main and interesting functions will be described in some detail below. We
begin however with an overview of how they are tied together by the shiny app.

Figure 5: Function Overview Flow Chart

#### specify_graph

The graphical interface is launched by specify_graph(), or preferably results
<\- specify_graph(). Once the shiny app is stopped, the input, output and
other useful information is returned by the function, so we recommend
assigning it to a variable so they are saved in the R-session and may easily
be further analyzed and processed. All further function calls will take place
automatically as the user interacts with the web interface. Thus, from a basic
user perspective, specify_graph is the main function. The core functionality
however is implemented in the functions analyze_graph, optimize_effect_2 and
their subroutines.

The JavaScript that handles the communication between the shiny server and the
input as the user draws a DAG through the web interface uses on the project
directed-graph-creator, an interactive tool for creating directed graphs,
created using d3.js and hosted at https://github.com/cjrd/directed-graph-
creator, which has been modified for the purpose of causal diagrams. The
modification binds the user inputs as they interact with the graph to shiny so
that the directed graph and its attributes set by the user are reactively
converted into an igraph-object for further processing. Since directed graphs
are common in many computational and statistical problems, this shiny
interface may also be valuable to many other R-package authors and maintainers
who may wish to provide their users with an accessible and intuitive way to
interact with their software.

The server listens to a reactive function that, as the user draws the DAG,
collects information about the current edges, collects and annotates vertices,
adds left- and right-side confounding, and returns an annotated igraph-object,
comprising information about the connectivity along with some additional
attributes; for each variable, its name, cardinality, latency-indicator and
side-indicator, and for each edge, a monotonicity-indicator and (to detect and
communicate violations on direction) a right-to-left-indicator. The server
meanwhile also monitors the DAG for any violation of the restriction that each
edge between $\mathcal{L}$ and $\mathcal{R}$ must go _from_ $\mathcal{L}$ _to_
$\mathcal{R}$, and if detected directly communicates this to the user through
a text message in the shiny app.

#### analyze_graph

The function analyze_graph takes a DAG (in the form of an igraph object),
optional constraints, and a string representing the causal effect of interest
and proceeds to construct and return a linear optimization problem (a
linearcausalproblem-object) from these inputs.

First, some basic data structures are created to keep track of the observed
variables, their possible values, the latent variables, and whether they are
in $\mathcal{L}$ or $\mathcal{R}$. Once these basic data-structures have been
created, the first task of the algorithm is to create the response function
variables (for each variable, observed or not, except $U_{\mathcal{L}}$ and
$U_{\mathcal{R}}$). Probabilities of these will be the entities $\mathbf{q}$
in which the objective function (representing the target causal effect) is
expressed and will constitute the points in the space it is optimized over,
where this space itself is constrained by the the relationships between them
and observed conditional probabilities $\mathbf{p}$.

#### create_response_function

The function create_response_function returns a list respvars that has a named
entry for each observed variable, containing its response function variable
and response function. If $X$ is an observed variable with $n$ response
functions, then they are enumerated by $\\{0,\dots,n-1\\}$. Its entry
respvars$X contains the response function variable $R_{X}$ of $X$, and is a
list with two entries. The first, respvars$X$index, is a vector containing all
the possible values of $R_{X}$, i.e., the integers $(0,\dots,n-1)$. The
second, respvars$X$values is itself a list with $n$ entries; each containing
the particular response function of $X$ corresponding to its index. Each such
response function is an actual R-function and may be evaluated by passing it
any possible values of the parents of $X$ as arguments.

Next, the response function variables are used in the creation of a matrix of
unobserved probabilities. Specifically the joint probabilities
$P(\mathbf{R}_{\mathcal{R}}=\mathbf{r}_{\mathcal{R}})$ for each possible
value-combination $\mathbf{r}_{\mathcal{R}}$ of the response function
variables $\mathbf{R}_{\mathcal{R}}$ of the right-side-variables
$\mathbf{W}_{\mathcal{R}}$. In (Sachs et al., 2022a), the possible value-
combinations $\mathbf{r}_{\mathcal{R}}$ are enumerated by
$\gamma\in\\{1,\dots,\aleph_{\mathcal{R}}\\}$ with corresponding probabilities
$q_{\gamma}:=P(\mathbf{R}_{\mathcal{R}}=\mathbf{r}_{\gamma})$ being components
of the vector $\mathbf{q}\in[0,1]^{\aleph_{\mathcal{R}}}$.

#### create_R_matrix

The constraints that the DAG and observed conditional probabilities
$\mathbf{p}$ (in p.vals) impose on the unobserved probabilities $\mathbf{q}$
(represented by variables) are linear. Specifically, there exists a matrix
whose entries are the coefficients relating p.vals to variables. This matrix
is called $P$ in (Sachs et al., 2022a), where its existence is guaranteed by
Proposition 2 and its construction is detailed in Algorithm 1, which is
implemented in the function create_R_matrix. This function returns back a list
with two entries; a vector of strings representing the linear constraints on
the unobserved $\mathbf{q}\in[0,1]^{\aleph_{\mathcal{R}}}$ imposed by and in
terms of the observed $\mathbf{p}\in[0,1]^{B}$ and the numeric matrix
$R\in\\{0,1\\}^{(B+1)\times\aleph_{\mathcal{R}}}$ of coefficients
corresponding to these constraints as well as the probabilistic ones and given
by $R=\begin{pmatrix}\mathbf{1}\\\ P\end{pmatrix}$ where
$P\in\\{0,1\\}^{B\times\aleph_{\mathcal{R}}}:\mathbf{p}=P\mathbf{q}$, so
$R\mathbf{q}=\begin{pmatrix}1\\\ \mathbf{p}\end{pmatrix}$.

This determines the constraint space as a compact convex polytope in
$\mathbf{q}$-space, i.e., in $\mathbb{R}^{\aleph_{\mathcal{R}}}$. To create
the matrix, we define a recursive function gee_r that takes two arguments; a
positive integer i being the index $i\in\\{1,\dots,n\\}$ of a variable
$W_{i}\in\mathcal{W}$ (i.e. the $i^{th}$ component of $\mathbf{W}$ or,
equivalently, the ith entry of obsvars) and a vector r being a value
$\mathbf{r}\in\nu(\mathbf{R})$ in the set $\nu(\mathbf{R})$ of all possible
value-vectors of the joint response function variable $\mathbf{R}$. This
recursive function is called for each variable in obsvars and for each
possible value of the response function variable vector. The base case is
reached if the variable has no parents, in which case the list corresponding
to the response function variable $R_{W_{i}}$ of $W_{i}$ is extracted from
respvars. From this list, the entry whose index matches the ith index of r
(i.e. the one corresponding to the response function variable value
$r_{i}=$r[i]) is extracted and finally its value, i.e., the corresponding
response function itself, is extracted and is evaluated on an empty list of
arguments, since it is a constant function and determined only by the value
$r_{i}$.

The recursive case is encountered when parents is non-empty. If so, then for
each parent in parents, its index in obsvars is determined and gee_r is
recursively called with the same vector r as first argument but now with this
particular index (i.e. that of the current parent) as second argument. The
numeric values returned by these recursive calls are then sequentially stored
in a vector lookin, whose entries are named by those in parents. Just as in
the base case, the response function corresponding to the particular value
$r_{i}$ of the response function variable $R_{W_{i}}$ (i.e. the response
function of the variable obsvars[i] that has the index r[i]) is extracted from
respvars and is now evaluated with arguments given by the list lookin. Note
that gee_r(r, i) corresponds to the value $w_{i}=g^{*}_{W_{i}}(\mathbf{r})$ in
(Sachs et al., 2022a).

Then the values that match the observed probabilities are recorded, the
corresponding entries in the current row of the matrix R are set to 1 and a
string representing the corresponding equation is constructed and added to the
vector of constraints.

#### parse_effect

Now that the constraint space has been determined, the objective function
representing the causal query needs to be specified as a linear function of
the components of $\mathbf{q}$, i.e., variables. First the causal query that
has been provided by the user as a text-string in effectt is passed to the
function parse_effect, which identifies its components including nested
counterfactuals and creates a data structure representing the causal query.
This structure includes nested lists which represent all interventional paths
to each outcome variable in the query.

Once the nested list effect is returned back to analyze_graph, it checks that
the requirements (see Proposition 3 in (Sachs et al., 2022a)) on the query are
fulfilled before creating the linear objective function. Despite these
regularity conditions, a large set of possible queries may be entered using
standard counterfactual notation, using syntax described in the accompanying
instruction text along with examples such as
$P(Y(M(X=0),X=1)=1)-P(Y(M(X=0),X=0)=1)$; the natural direct effect (Pearl,
2001) of a binary exposure $X$ at level $M=0$ on a binary outcome $Y$ _not_
going through the mediator $M$, in the presence of unmeasured confounding
between $M$ and $Y$ (Sjölander, 2009).

#### create_effect_vector

Now that the required characteristics of the query have been established, the
corresponding objective function will be constructed by the function
create_effect_vector which returns a list var.eff of string-vectors; one for
each term in the query. Each such vector contains the names (strings in
variables) of the response function variables of the right-side (i.e. the
components of $\mathbf{q}$) whose sum corresponds the that particular term.
The function create_effect_vector implements Algorithm 2 of (Sachs et al.,
2022a) with the additional feature that if the user has entered a query that
is incomplete in the sense that there are omitted mediating variables on paths
from base/intervention variables to the outcome variable, then this is
interpreted as the user intending the effects of the base/intervention
variables to be propagated through the mediators, so that they are set to
their “natural” values under this intervention. These mediators are detected
and their values are set accordingly.

We define a recursive function gee_rA that takes three arguments; a positive
integer i (the index $i$ of a variable $W_{i}\in\mathcal{W}=$obsvars), a
vector r (a value $\mathbf{r}\in\nu(\mathbf{R})$ in the set $\nu(\mathbf{R})$
of all possible value-vectors of the joint response function variable
$\mathbf{R}$) and a string path that represents an interventional path and is
of the form “X -> … -> Y” if not NULL. The base case is reached either if path
is non-NULL and corresponds to a path to the intervention set or if parents is
empty. In the former case, the corresponding numeric intervention-value is
returned, and in the latter case, the value of the corresponding response
function called on the empty list of arguments is returned just as in the base
case of gee_r. The recursive case is encountered when path is NULL and parents
is non-empty. This recursion proceeds just as in gee_r, but now rather with a
recursive call to gee_rA, whose third argument is now path = paste(gu, "->",
path) where the string in gu is the name of the parent variable in parents
whose index i in obsvars is the second argument of this recursive call. This
construction traces the full path taken from the outcome of interest to the
variable being intervened upon. Note that gee_rA(r, i, path) corresponds to
the value $w_{i}=h^{A_{i}}_{W_{i}}(\mathbf{r},W_{i})$ in (Sachs et al.,
2022a). A matrix is now created just as in the observational case, but this
time using gee_rA instead of gee_r .

#### optimize_effect_2

Once the constraints on $\mathbf{q}$ as well as the effect of interest in
terms of $\mathbf{q}$ have been established, it remains only to optimize this
expression over the constraint space. Here, $\mathbf{c}$ denotes the constant
gradient vector of the linear objective function and $P$ denotes the
coefficient matrix of the linear restrictions on $\mathbf{q}$ in terms of
$\mathbf{q}$ imposed by the causal DAG. By adding the probabilistic
constraints on $\mathbf{q}$ we have arrived at e.g. the following linear
program giving a tight lower bound on the average causal effect
$\theta_{\mathbf{q}}=P\\{Y(X=1)=1\\}-P\\{Y(X=0)=1\\}$ in the simple
instrumental variable problem of the introductory section:

$\displaystyle\min_{\mathbf{q}}\theta_{\mathbf{q}}$
$\displaystyle=\min\\{\mathbf{c}^{\top}\mathbf{q}\mid\mathbf{q}\in\mathbb{R}^{16},\mathbf{q}\geq\mathbf{0}_{16\times
1},\mathbf{1}_{1\times 16}\mathbf{q}=1,P\mathbf{q}=\mathbf{q}\\}$
$\displaystyle=\max\\{\begin{pmatrix}1&\mathbf{q}^{\top}\end{pmatrix}\mathbf{y}\mid\mathbf{y}\in\mathbb{R}^{9},\mathbf{y}\geq\mathbf{0}_{9\times
1},\begin{pmatrix}\mathbf{1}_{16\times
1}&P^{\top}\end{pmatrix}\leq\mathbf{c}\\}$
$\displaystyle=\max\\{\begin{pmatrix}1&\mathbf{q}^{\top}\end{pmatrix}\bar{\mathbf{y}}\mid\bar{\mathbf{y}}\text{
is a vertex of
}\\{\mathbf{y}\in\mathbb{R}^{9}\mid\mathbf{y}\geq\mathbf{0}_{9\times
1},R^{\top}\leq\mathbf{c}\\}\\}$

Since we allow the user to provide additional linear inequality constraints
(e.g. it may be quite reasonable to assume the proportion of “defiers” in the
study population of our example to be quite low), the actual primal and dual
linear programs may look slightly more complicated, but this small example
still captures the essentials.
In general, given the matrix of linear constraints on the observable
probabilities implied by the DAG and an optional user-provided matrix
inequality, we construct the coefficient matrix and right hand side vector of
the dual polytope.

The optimization via vertex enumeration step in causaloptim is implemented in
the function optimize_effect_2 which uses the double description method for
vertex enumeration, as implemented in the rcdd package (Geyer et al., 2021).
This step of vertex enumeration has previously been the major computational
bottleneck. The approach is now based on cddlib
(https://people.inf.ethz.ch/fukudak/cdd_home/), which has an implementation of
the Double Description Method (dd). Any convex polytope can be dually
described as either an intersection of half-planes (which is the form we get
our dual constraint space in) or as a minimal set of vertices of which it is
the convex hull (which is the form we want it in) and the dd algorithm
efficiently converts between these two descriptions. cddlib also uses exact
rational arithmetic, so there is no need to worry about any numerical
instability issues. The vertices of the dual polytope are obtained and stored
as rows of a matrix with hrep <\- rcdd::makeH(a1, b1); vrep <\-
rcdd::scdd(hrep); vertices <\- vrep$output[vrep$output[, 1] == 0 &
vrep$output[, 2] == 1, -c(1, 2), drop=FALSE].

The rest is simply a matter of plugging them into the dual objective function,
evaluating the expression and presenting the results. The first part of this
is done by apply(vertices, 1, function(y) evaluate_objective(c1_num, p, y))
(here (c1_num,p)$=(\begin{pmatrix}b_{\ell}^{\top}&1\end{pmatrix},p)$ separates
the dual objective gradient into its numeric and symbolic parts).

causaloptim also contains a precursor to to optimize_effect_2, called
optimize_effect. This legacy function uses the original optimization procedure
written in C++ by Alexander Balke and involves linear program formulation
followed by the vertex enumeration algorithm of (Mattheiss, 1973). This has
worked well for very simple settings but has struggled severely with even
remotely complex ones and thus been insufficient for the ambitions of
causaloptim. The efficiency gains of optimize_effect_2 over the legacy code
have reduced the computation time for several setting from hours to
milliseconds.

## Numeric Examples

### A Mediation Analysis

In (Sjölander, 2009), the author derives bounds on natural direct effects in
the presence of confounded intermediate variables and applies them to data
from the Lipid Research Clinics Coronary Primary Prevention Trial (Freedman et
al., 1992), where subjects were randomized to cholestyramine treatment and
presence of coronary heart disease events as well as levels of cholesterol
were recorded after a 1-year follow-up period. We let $X$ be a binary
treatment indicator, with $X=0$ indicating actual cholestyramine treatment and
$X=1$ indicating placebo. We further let $Y$ be an indicator of the occurrence
of coronary heart disease events within follow-up, with $Y=0$ indicating
event-free follow-up and $Y=1$ indicating an event. We finally let $M$ be a
dichotomized (cut-off at $280\ mg/dl$) cholesterol level indicator, with $M=0$
indicating levels $<280\ mg/dl$ and $M=1$ indicating levels $\geq 280\ mg/dl$.
The causal assumptions are summarized in the DAG shown in Figure 6, where
$U_{l}$ and $U_{r}$ are unmeasured and the latter confounds the effect of $M$
on $Y$.

Figure 6: Causal DAG for mediation example

b <\- igraph::graph_from_literal(X -+ Y, X -+ M, M -+ Y, Ul -+ X, Ur -+ Y, Ur
-+ M)V(b)$leftside <\- c(1, 0, 0, 1, 0)V(b)$latent <\- c(0, 0, 0, 1,
1)V(b)$nvals <\- c(2, 2, 2, 2, 2)E(b)$rlconnect <\- c(0, 0, 0, 0, 0,
0)E(b)$edge.monotone <\- c(0, 0, 0, 0, 0, 0)

Using the data from Table IV of (Sjölander, 2009), we compute the observed
conditional probabilities.

# parameters of the form pab_c, which represents# the probability P(Y = a, M = b | X = c)p00_0 <\- 1426/1888 # P(Y=0,M=0|X=0)p10_0 <\- 97/1888 # P(Y=1,M=0|X=0)p01_0 <\- 332/1888 # P(Y=0,M=1|X=0)p11_0 <\- 33/1888 # P(Y=1,M=1|X=0)p00_1 <\- 1081/1918 # P(Y=0,M=0|X=1)p10_1 <\- 86/1918 # P(Y=1,M=0|X=1)p01_1 <\- 669/1918 # P(Y=0,M=1|X=1)p11_1 <\- 82/1918 # P(Y=1,M=1|X=1)

We proceed to compute bounds on the controlled direct effect
$CDE(0)=P(Y(M=0,X=1)=1)-P(Y(M=0,X=0)=1)$ of $X$ on $Y$ not passing through $M$
at level $M=0$, the controlled direct effect
$CDE(1)=P(Y(M=1,X=1)=1)-P(Y(M=1,X=0)=1)$ at level $M=1$, the natural direct
effect $NDE(0)=P(Y(M(X=0),X=1)=1)-P(Y(M(X=0),X=0)=1)$ of $X$ on $Y$ at level
$X=0$ and the natural direct effect
$NDE(1)=P(Y(M(X=1),X=1)=1)-P(Y(M(X=1),X=0)=1)$ at level $X=1$.

CDE0_query <\- "p{Y(M = 0, X = 1) = 1} - p{Y(M = 0, X = 0) = 1}"CDE0_obj <\-
analyze_graph(b, constraints = NULL, effectt = CDE0_query)CDE0_bounds <\-
optimize_effect_2(CDE0_obj)CDE0_boundsfunction <\- interpret_bounds(bounds =
CDE0_bounds$bounds, parameters = CDE0_obj$parameters)CDE0_numericbounds <\-
CDE0_boundsfunction(p00_0 = p00_0, p00_1 = p00_1, p10_0 = p10_0, p10_1 =
p10_1, p01_0 = p01_0, p01_1 = p01_1, p11_0 = p11_0, p11_1 = p11_1)CDE1_query
<\- "p{Y(M = 1, X = 1) = 1} - p{Y(M = 1, X = 0) = 1}"CDE1_obj <\-
update_effect(CDE0_obj, effectt = CDE1_query)CDE1_bounds <\-
optimize_effect_2(CDE1_obj)CDE1_boundsfunction <\- interpret_bounds(bounds =
CDE1_bounds$bounds, parameters = CDE1_obj$parameters)CDE1_numericbounds <\-
CDE1_boundsfunction(p00_0 = p00_0, p00_1 = p00_1, p10_0 = p10_0, p10_1 =
p10_1, p01_0 = p01_0, p01_1 = p01_1, p11_0 = p11_0, p11_1 = p11_1)NDE0_query
<\- "p{Y(M(X = 0), X = 1) = 1} - p{Y(M(X = 0), X = 0) = 1}"NDE0_obj <\-
update_effect(CDE0_obj, effectt = NDE0_query)NDE0_bounds <\-
optimize_effect_2(NDE0_obj)NDE0_boundsfunction <\- interpret_bounds(bounds =
NDE0_bounds$bounds, parameters = NDE0_obj$parameters)NDE0_numericbounds <\-
NDE0_boundsfunction(p00_0 = p00_0, p00_1 = p00_1, p10_0 = p10_0, p10_1 =
p10_1, p01_0 = p01_0, p01_1 = p01_1, p11_0 = p11_0, p11_1 = p11_1)NDE1_query
<\- "p{Y(M(X = 1), X = 1) = 1} - p{Y(M(X = 1), X = 0) = 1}"NDE1_obj <\-
update_effect(CDE0_obj, effectt = NDE1_query)NDE1_bounds <\-
optimize_effect_2(NDE1_obj)NDE1_boundsfunction <\- interpret_bounds(bounds =
NDE1_bounds$bounds, parameters = NDE1_obj$parameters)NDE1_numericbounds <\-
NDE1_boundsfunction(p00_0 = p00_0, p00_1 = p00_1, p10_0 = p10_0, p10_1 =
p10_1, p01_0 = p01_0, p01_1 = p01_1, p11_0 = p11_0, p11_1 = p11_1)

We obtain the same symbolic bounds as (Sjölander, 2009) and the resulting
numeric bounds are given in Table 1 which of course agree with those of Table
V in (Sjölander, 2009).

Table 1: Bounds on the controlled and natural direct effects. | lower | upper
---|---|---
CDE(0) | -0.20 | 0.39
CDE(1) | -0.78 | 0.63
NDE(0) | -0.07 | 0.56
NDE(1) | -0.55 | 0.09

### A Mendelian Randomization Study of the Effect of Homocysteine on
Cardiovascular Disease

Mendelian randomization (Davey Smith and Ebrahim, 2003) assumes certain
genotypes may serve as suitable instrumental variables for investigating the
causal effect of an associated phenotype on some disease outcome.

In (Palmer, 2011), the authors investigate the effect of homocysteine on
cardiovascular disease using the 677CT polymorphism (rs1801133) in the
Methylenetetrahydrofolate Reductase gene as an instrument. They use
observational data from (Meleady et al., 2003) in which the outcome is binary,
the treatment has been made binary by a suitably chosen cut-off at $15\mu
mol/L$, and the instrument is ternary (this polymorphism can take three
possible genotype values).

With $X$ denoting the treatment, $Y$ the outcome and $Z$ the instrument, the
conditional probabilities are given as follows.

params <\- list(p00_0 = 0.83, p00_1 = 0.88, p00_2 = 0.72, p10_0 = 0.11, p10_1
= 0.05, p10_2 = 0.20, p01_0 = 0.05, p01_1 = 0.06, p01_2 = 0.05, p11_0 = 0.01,
p11_1 = 0.01, p11_2 = 0.03)

The computation using causaloptim is done using the following code.

# Input causal DAGb <\- graph_from_literal(Z -+ X, Ul -+ Z, X -+ Y, Ur -+ X,
Ur -+ Y)V(b)$leftside <\- c(1, 0, 1, 0, 0)V(b)$latent <\- c(0, 0, 1, 0,
1)V(b)$nvals <\- c(3, 2, 2, 2, 2)E(b)$rlconnect <\- c(0, 0, 0, 0,
0)E(b)$edge.monotone <\- c(0, 0, 0, 0, 0)# Construct causal problemobj <\-
analyze_graph(b, constraints = NULL, effectt = "p{Y(X = 1) = 1} - p{Y(X = 0) =
1}")# Compute bounds on querybounds <\- optimize_effect_2(obj)# Construct
bounds as function of parametersboundsfunction <\- interpret_bounds(bounds =
bounds$bounds, parameters = obj$parameters)# Insert observed conditional
probabilitiesnumericbounds <\- do.call(boundsfunction,
as.list(params))round(numericbounds, 2)

#> lower upper#> 1 -0.09 0.74

Our computed bounds agree with those computed using bpbounds as well as those
estimated using Theorem 2 of (Richardson and Robins, 2014), who independently
derived expressions for tight bounds that are applicable to this setting.

## Summary and Discussion

The methods and algorithms described in (Sachs et al., 2022a) to compute
symbolic expressions for bounds on non-identifiable causal effects are
implemented in the package causaloptim. Our aim was to provide a user-friendly
interface to these methods with a graphical interface to draw DAGs,
specification of causal effects using standard notation for potential
outcomes, and an efficient implementation of vertex enumeration to reduce
computation times. These methods are applicable to a wide variety of causal
inference problems which appear in biomedical research, economics, social
sciences and more. Aside from the graphical interface, programming with the
package is encouraged to promote reproducibility and advanced use. Our package
includes automated unit tests and also tests for correctness by comparing the
symbolic bounds derived using our program to independently derived bounds in
particular settings.

Our implementation uses a novel approach to draw DAGs using JavaScript in a
web browser that can then be passed to R using shiny. This graphical approach
can be adapted and used in other settings where graphs need to be specified
and computed on, such as other causal inference settings, networks, and multi-
state models. Other algorithms and data structures that could be more broadly
useful include the representation of structural equations as R functions,
recursive evaluation of response functions, and parsing of string equations
for causal effects and constraints.

## References

  * Balke and Pearl (1997) A. Balke and J. Pearl.  Bounds on treatment effects from studies with imperfect compliance.  _Journal of the American Statistical Association_ , 92(439):1171–1176, 1997.
  * Cai et al. (2008) Z. Cai, M. Kuroki, J. Pearl, and J. Tian.  Bounds on direct effects in the presence of confounded intermediate variables.  _Biometrics_ , 64(3):695–701, 2008.
  * Chang et al. (2021) W. Chang, J. Cheng, J. Allaire, C. Sievert, B. Schloerke, Y. Xie, J. Allen, J. McPherson, A. Dipert, and B. Borges.  _shiny: Web Application Framework for R_ , 2021.  URL https://CRAN.R-project.org/package=shiny.  R package version 1.7.1.
  * Csardi and Nepusz (2006) G. Csardi and T. Nepusz.  The igraph software package for complex network research.  _InterJournal_ , Complex Systems:1695, 2006.  URL https://igraph.org.
  * Davey Smith and Ebrahim (2003) G. Davey Smith and S. Ebrahim.  ‘Mendelian randomization’: can genetic epidemiology contribute to understanding environmental determinants of disease?*.  _International Journal of Epidemiology_ , 32(1):1–22, 02 2003.  ISSN 0300-5771.  doi: 10.1093/ije/dyg070.  URL https://doi.org/10.1093/ije/dyg070.
  * Freedman et al. (1992) L. S. Freedman, B. I. Graubard, and A. Schatzkin.  Statistical validation of intermediate endpoints for chronic diseases.  _Statistics in Medicine_ , 11(2):167–178, 1992\.  doi: https://doi.org/10.1002/sim.4780110204.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.4780110204.
  * Gabriel et al. (2020) E. E. Gabriel, M. C. Sachs, and A. Sjölander.  Causal bounds for outcome-dependent sampling in observational studies.  _Journal of the American Statistical Association_ , pages 1–12, 2020\.
  * Gabriel et al. (2021) E. E. Gabriel, A. Sjölander, and M. C. Sachs.  Nonparametric bounds for causal effects in imperfect randomized experiments.  _Journal of the American Statistical Association_ , pages 1–9, 2021\.
  * Gabriel et al. (2022) E. E. Gabriel, M. C. Sachs, and A. Sjölander.  Sharp nonparametric bounds for decomposition effects with two binary mediators.  _Journal of the American Statistical Association_ , pages 1–8, 2022\.
  * Geyer et al. (2021) C. J. Geyer, G. D. Meeden, and incorporates code from cddlib written by Komei Fukuda.  _rcdd: Computational Geometry_ , 2021.  URL https://CRAN.R-project.org/package=rcdd.  R package version 1.5.
  * Greenland et al. (1999) S. Greenland, J. Pearl, and J. M. Robins.  Causal diagrams for epidemiologic research.  _Epidemiology_ , pages 37–48, 1999.
  * Mattheiss (1973) T. H. Mattheiss.  An algorithm for determining irrelevant constraints and all vertices in systems of linear inequalities.  _Operations Research_ , 21(1):247–260, 1973.
  * Meleady et al. (2003) R. Meleady, P. M. Ueland, H. Blom, A. S. Whitehead, H. Refsum, L. E. Daly, S. E. Vollset, C. Donohue, B. Giesendorf, I. M. Graham, A. Ulvik, Y. Zhang, and A.-L. Bjorke Monsen.  Thermolabile methylenetetrahydrofolate reductase, homocysteine, and cardiovascular disease risk: the European Concerted Action Project.  _The American journal of clinical nutrition_ , 77(1):63–70, Jan. 2003.  ISSN 0002-9165.  doi: 10.1093/ajcn/77.1.63.  Place: United States.
  * Palmer et al. (2018) T. Palmer, R. Ramsahai, V. Didelez, and N. Sheehan.  _bpbounds: R package implementing Balke-Pearl bounds for the average causal effect_ , 2018.  URL https://github.com/remlapmot/bpbounds.
  * Palmer (2011) T. M. Palmer.  Nonparametric bounds for the causal effect in a binary instrumental-variable model.  _Stata Journal_ , 11(3):345–367(23), 2011.  URL //<?echo(www)%****␣main.bbl␣Line␣125␣****?>.stata-journal.com/article.html?article=st0232.
  * Pearl (2001) J. Pearl.  Direct and indirect effects.  In _Proceedings of the Seventeenth Conference on Uncertainty and Artificial Intelligence, 2001_ , pages 411–420. Morgan Kaufman, 2001.
  * Pearl (2009) J. Pearl.  _Causality_.  Cambridge university press, 2009.
  * Richardson and Robins (2014) T. S. Richardson and J. M. Robins.  Ace bounds; sems with equilibrium conditions.  _Statistical Science_ , 29(3):363–366, 2014.
  * Sachs et al. (2022a) M. C. Sachs, G. Jonzon, A. Sjölander, and E. E. Gabriel.  A general method for deriving tight symbolic bounds on causal effects.  _Journal of Computational and Graphical Statistics_ , 0(0):1–10, 2022a.  doi: 10.1080/10618600.2022.2071905.  URL https://doi.org/10.1080/10618600.2022.2071905.
  * Sachs et al. (2022b) M. C. Sachs, A. Sjölander, and E. E. Gabriel.  _causaloptim: An Interface to Specify Causal Graphs and Compute Bounds on Causal Effects_ , 2022b.  URL https://github.com/sachsmc/causaloptim.  R package version 0.9.2.
  * Sjölander (2009) A. Sjölander.  Bounds on natural direct effects in the presence of confounded intermediate variables.  _Statistics in Medicine_ , 28(4):558–571, 2009\.
  * Sjölander et al. (2014) A. Sjölander, W. Lee, H. Källberg, and Y. Pawitan.  Bounds on causal interactions for binary outcomes.  _Biometrics_ , 70(3):500–505, 2014.
  * Sjölander (2009) A. Sjölander.  Bounds on natural direct effects in the presence of confounded intermediate variables.  _Statistics in Medicine_ , 28(4):558–571, 2009\.  doi: https://doi.org/10.1002/sim.3493.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.3493.

_Gustav Jonzon
Department of Medical Epidemiology and Biostatistics, Karolinska Institutet

https://ki.se/meb
<EMAIL_ADDRESS>
_Michael C Sachs
Department of Public Health, University of Copenhagen

https://biostat.ku.dk/
ORCiD: 0000-0002-1279-8676
<EMAIL_ADDRESS>
_Erin E Gabriel
Department of Public Health, University of Copenhagen

https://biostat.ku.dk/
ORCiD: 0000-0002-0504-8404
<EMAIL_ADDRESS>