#  How Viable Is a QCD Axion near 10 MeV?

Sudhakantha Girmohanta, Shota Nakagawa, Yuichiro Nakai and Junxuan Xu
Tsung-Dao Lee Institute, Shanghai Jiao Tong University,
No. 1 Lisuo Road, Pudong New Area, Shanghai, 201210, China
School of Physics and Astronomy, Shanghai Jiao Tong University,
800 Dongchuan Road, Shanghai, 200240, China

###### Abstract

There has been an attempt to revive the visible QCD axion at the 10 MeV scale
assuming that it exclusively couples to the first-generation quarks and the
electron. This variant of the QCD axion is claimed to remain
phenomenologically viable, partly due to a clever model construction that
induces tree-level pion-phobia and exploits uncertainties inherent in the
chiral perturbation theory. We confront this model with the cosmological
domain wall problem, the quality issue and constraints arising from the
electron electric dipole moment. It is also pointed out that the gluon loop-
generated axion-top coupling can provide a very large contribution to rare
$B$-meson decays, such that the present LHCb data for $B^{0}\to
K^{*0}e^{+}e^{-}$ rule out the model for the axion mass larger than 30 MeV.
There is a strong motivation for pushing the experimental analysis of $B\to
K^{(*)}e^{+}e^{-}$ to a lower $e^{+}e^{-}$ invariant mass window, which will
conclusively determine the fate of the model, as its contribution to this
branching ratio significantly exceeds the Standard Model prediction.

##  I Introduction

The Standard Model (SM) has achieved incredible success in accounting for the
experimental data observed to date, yet it leaves several theoretical issues
unresolved. The strong CP problem, namely why the CP-violating parameter
originating from the topological vacuum nature in QCD ($\bar{\theta}$) is
smaller than $\sim 10^{-10}$, remains one of the outstanding puzzles.111 It
was noted however that, setting $\bar{\theta}=0$ at some UV scale, the
Cabibbo–Kobayashi–Maskawa (CKM) phase only generates $\bar{\theta}\sim
10^{-16}$ Ellis and Gaillard (1979). As the viability of the massless up quark
has been ruled out from the lattice simulations Workman  _et al._ (2022);
Alexandrou  _et al._ (2020), the most popular dynamical solution is due to
Peccei and Quinn (PQ) where a global $U(1)_{\rm PQ}$ symmetry is spontaneously
broken, giving rise to a (pseudo-)Nambu-Goldstone boson $a$, the axion Peccei
and Quinn (1977); Weinberg (1978); Wilczek (1978). Once the axion settles down
into the minimum of the potential induced by the strong interaction, the
strong CP problem is dynamically solved.

The original visible QCD axion, with a decay constant $f_{a}$ at the
electroweak scale, was soon ruled out as a result of constraints from rare
meson decays, and subsequently, attention was shifted to the invisible axion
with a much larger $f_{a}$ (and smaller mass) Kim (1979); Dine  _et al._
(1981). Nevertheless, it has been claimed in Ref. Alves and Weiner (2018) that
the exclusion for the axion at the 10 MeV scale is not conclusive, and a
particular pion-phobic variant, where the axion only couples to the first-
generation quarks and the electron by model construction, can still be viable.
This possibility is also motivated by the still indecisive claim of an
$e^{+}e^{-}$ excess in the 8Be nuclear de-excitation spectrum Krasznahorkay
_et al._ (2016). Progress has been made toward achieving a realistic UV
complete model of this QCD axion variant in Ref. Liu  _et al._ (2021), where
the PQ symmetry breaking is intertwined with the electroweak symmetry breaking
such that the observed CKM mixing matrix can be reproduced.

In the present paper, we discuss several new theoretical and experimental
consequences of the aforementioned visible QCD axion. As the decay constant
$f_{a}$ is around GeV for the axion at the 10 MeV scale, the PQ symmetry
breaking is likely to take place in the Universe after inflation. Further,
since this variant necessarily couples to the first-generation quarks, the
domain wall number is larger than 1. Hence, the PQ phase transition generates
a string-domain wall network, which results in a catastrophic deviation from
the standard cosmology unless it decays rapidly. The introduction of a
sufficiently large bias term in the form of operators that explicitly break
the PQ symmetry enables the network to decay. However, the strength of these
operators must be suppressed enough not to spoil the quality of the PQ
symmetry to achieve $\bar{\theta}\lesssim 10^{-10}$. We will study the
interplay of these two conditions for the axion at the 10 MeV scale.

The $U(1)_{\rm PQ}$ charge of the up and down quarks must be different to
realize the required pion-phobia Krauss and Nash (1988). Therefore, in a
renormalizable UV completion, one generically needs multiple additional Higgs
doublets, inevitably introducing new CP-violating phases in the scalar
potential. Working with the benchmark model in Ref. Liu  _et al._ (2021), we
will present those CP-violating phases. Then, electric dipole moments (EDMs)
are generated as a result of the new complex phases. We will find that the
electron EDM generated at one/two-loop in the model gives the most stringent
constraint on the new physics mass scale.

Although the tree-level interactions of the visible axion to the higher
generation quarks are zero by the model construction, they are generated via
gluon loops. While these loop generated interactions to the $b$ and
$c$-quarks, constrained by the decays of $J/\Psi$ and $\Upsilon$, have been
considered in the literature Alves and Weiner (2018); Alves (2021), the
similar coupling to the top quark was neglected, because we cannot produce the
corresponding topped bound state, and the top quark decays even before it can
hadronize. Here, however, we will argue that the axion-top coupling generated
through a gluon loop gives a dominant new physics contribution to decays like
$B\to Ka$. In particular, as the present axion variant promptly decays to a
$e^{+}e^{-}$ pair to be consistent with the beam dump experiments, the rare
decay mode $B^{0}\to K^{*0}e^{+}e^{-}$ can be used to constrain the model. We
show that the contribution of the visible QCD axion to the branching ratio of
the decay is much larger than that of the SM. Then, the current LHCb data Aaij
_et al._ (2013), applicable to the invariant mass of the $e^{+}e^{-}$ pair in
the range of $30$-$1000$ MeV, rule out the model for the axion mass $m_{a}>30$
MeV. We emphasize that there is a strong motivation for pushing the analysis
for this mode into the lower invariant mass regime, which will conclusively
determine the fate of the present visible axion variant.

The rest of the paper is organized as follows. In Sec. II, we briefly give a
review of the QCD axion variant at the 10 MeV scale and its realistic UV
completion. We discuss the quality and domain wall problems in Sec. III. Sec.
IV deals with the induced electron EDM, while Sec. V discusses the constraints
from the loop-mediated $B$-meson decay. Finally, we conclude in Sec. VI.

##  II QCD axion near 10 MeV

Let us begin with a review of the QCD axion variant with mass of
$\mathcal{O}(10)\ {\rm MeV}$ discussed in Refs. Alves and Weiner (2018); Alves
(2021); Liu  _et al._ (2021). It had been widely believed that such a visible
axion model can be easily excluded by beam dump experiments Donnelly  _et al._
(1978) or measurements of rare heavy meson decays Wilczek (1977); Hall and
Wise (1981). However, the authors of Ref. Alves and Weiner (2018) claimed a
possible way out of those severe constraints on the model, which is composed
of three points: (i) the axion must couple to only the first-generation of the
SM fermions, (ii) the axion is required to be pion-phobic Krauss and Nash
(1988), and (iii) there exist large uncertainties in theory and experiment.
For point (i), only the axion couplings to the first-generation quarks ($u,d$)
are allowed to avoid exotic decays of the heavy quarkonia, e.g. $J/\psi$ and
$\Upsilon$ Albrecht  _et al._ (1986); Hsueh and Palestini (1992); Armstrong
_et al._ (1993); Bai  _et al._ (1995); Adams  _et al._ (2006). The axion must
couple to the electron as otherwise it would become too long-lived and would
have been already detected by beam dump experiments (the lifetime is bounded
as $\tau\lesssim 10^{-13}\sec$) Bergsma  _et al._ (1985); Bjorken  _et al._
(1988); Davier and Nguyen Ngoc (1989); Blumlein  _et al._ (1991). The axion
coupling to the muon is prohibited by the measurement of the muon magnetic
dipole moment unless the decay constant $f_{a}$ is larger than the electroweak
scale. Regarding point (ii), pion-phobia enables the axion-pion mixing to be
suppressed, leading to the tiny rate of an exotic decay of the pion. The
SINDRUM collaboration Eichler  _et al._ (1986) searches for the process
$\pi^{+}\rightarrow e^{+}\nu_{e}a(\rightarrow e^{+}e^{-})$ and provides a
constraint on the branching ratio, resulting in the upper bound on the axion-
pion mixing. To realize pion-phobia, one can take special PQ charge
assignments of $u$ and $d$,

$\displaystyle\frac{Q_{u}}{Q_{d}}=2,$ (1)

leading to the cancellation due to the mass ratio between $u$ and $d$. For
point (iii), the axion-$\eta$ meson mixing induces a $K^{+}\to\pi^{+}a$ decay,
however, the inherent uncertainties in the chiral perturbation theory have
been used to argue that the QCD axion variant remains viable Alves and Weiner
(2018).

For other requirements, we should mention the bound from the electron magnetic
moment, $(g-2)_{e}$. The axion-electron and -photon couplings induce loop
contributions to $(g-2)_{e}$ and thus receive a constraint from the current
measurements.222 There is a discrepancy of values of $(g-2)_{e}$ calculated by
using values of the fine-structure constant measured by two independent
experiments. While the cesium recoil experiment Parker  _et al._ (2018) leads
to a smaller value of $(g-2)_{e}$, the experiment using the rubidium atom
Morel  _et al._ (2020) leads to a larger value. The recent measurement of
$(g-2)_{e}$ Fan  _et al._ (2023) does not solve the tension. In addition, the
axion coupling to the electron is constrained by KLOE Anastasi  _et al._
(2015) and NA64 experiments Banerjee  _et al._ (2020); Depero  _et al._
(2020). Taking account of both bounds, we obtain the allowed range of the PQ
charge of the electron, $1/3\lesssim Q_{e}\lesssim 2$.

| $H$ | $H_{u}$ | $H_{d}$ | $H_{e}$ | $\Phi_{u}$ | $\Phi_{d}$ | $\Phi_{e}$
---|---|---|---|---|---|---|---
$SU(2)_{L}$ | $\bm{2}$ | $\bm{2}$ | $\bm{2}$ | $\bm{2}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$
$U(1)_{Y}$ | $1/2$ | $-1/2$ | $1/2$ | $1/2$ | $0$ | $0$ | $0$
$U(1)_{\rm PQ}$ | $0$ | $-Q_{u}$ | $-Q_{d}$ | $-Q_{e}$ | $-Q_{u}$ | $-Q_{d}$ | $-Q_{e}$
Table 1: The charge assignments of scalar fields.

We now consider a UV model to satisfy the requirements, following the
discussion of Ref. Liu  _et al._ (2021). The model introduces four doublet
Higgs fields $H,H_{u},H_{d},H_{e}$ and three SM-singlet scalar fields
$\Phi_{u},\Phi_{d},\Phi_{e}$. Here, $H$ corresponds to the usual SM Higgs
doublet. The charge assignments of these scalar fields are summarized in Tab.
1, where we define the PQ charges of the right-handed up quark $u_{R}$, down
quark $d_{R}$ and electron $e_{R}$ as $Q_{u},Q_{d}$ and $Q_{e}$, respectively.
In the following discussion, we set them as

$\displaystyle Q_{u}=2,\quad Q_{d}=1,\quad Q_{e}=\frac{1}{n},$ (2)

with $n=2,3$. This choice makes the axion pion-phobic and evades the bound on
$Q_{e}$ conservatively. The color anomaly is estimated as $(Q_{u}+Q_{d})/2$,
and the domain wall number is three. Yukawa interactions with the new Higgs
doublets and the SM one are respectively given by

$\displaystyle\mathcal{L}_{\rm PQ}^{Y}$ $\displaystyle=$
$\displaystyle-\sum_{i=1,2,3}\left(\bar{Q}^{i}Y_{u}^{i1}H_{u}u_{R}^{1}+\bar{Q}^{i}Y_{d}^{i1}H_{d}d_{R}^{1}\right.$
(3) $\displaystyle+$
$\displaystyle\left.\bar{L}^{i}Y_{e}^{i1}H_{e}e_{R}^{1}\right)+{\rm h.c.},$
$\displaystyle\mathcal{L}_{\rm SM}^{Y}$ $\displaystyle=$
$\displaystyle-\sum_{i=1,2,3}\sum_{j=2,3}\left(\bar{Q}^{i}Y_{u}^{ij}\tilde{H}u_{R}^{j}+\bar{Q}^{i}Y_{d}^{ij}Hd_{R}^{j}\right.$
(4) $\displaystyle+$
$\displaystyle\left.\bar{L}^{i}Y_{e}^{ij}He_{R}^{j}\right)+{\rm h.c.},$

where $Q,L$ respectively denote the left-handed quark and lepton, the indices
$i,j$ represent the generation and $\tilde{H}\equiv i\sigma_{2}H^{\ast}$. The
first Lagrangian involves the first generation right-handed fermions, while
the second one is relevant only to the second and third generation right-
handed fermions. With these Yukawa interactions, the CKM matrix can be derived
successfully Liu  _et al._ (2021).

The potential terms of the scalar fields to break the electroweak and PQ
symmetries include

$\displaystyle V_{\rm PQ}$ $\displaystyle=$
$\displaystyle\bigl{(}A_{1}HH_{u}\Phi_{u}^{*}+A_{2}HH_{d}^{\dagger}\Phi_{d}+A_{3}H_{e}H^{\dagger}\Phi_{e}^{*}$
(5) $\displaystyle+$ $\displaystyle
A_{4}\Phi_{u}^{*}\Phi_{d}^{2}+A_{5}\Phi_{d}^{*}\Phi_{e}^{n}\bigr{)}+\rm{h.c.},$
$\displaystyle V_{\rm dia}$ $\displaystyle=$
$\displaystyle\sum_{\Psi}-\mu_{\Psi}\Psi^{\dagger}\Psi+\lambda_{\Psi}(\Psi^{\dagger}\Psi)^{2}.$
(6)

Here, $\Psi=H,H_{u},H_{d},H_{e},\Phi_{u},\Phi_{d},\Phi_{e}$ and
$A_{1},A_{2},A_{3},A_{4}$ and $\mu_{\Psi}$ are the parameters with mass
dimension one and two, respectively. The mass dimension of the parameter
$A_{5}$ depends on the choice of $n$, and $\lambda_{\Psi}$ is a dimensionless
coupling. Although other PQ-symmetric operators can be written down, they are
not relevant for the PQ mechanism and omitted here. However, in Sec. IV, we
will show that such extra terms induce sizable contributions to EDMs. In terms
of seven neutral pseudoscalars,
$\vec{\Phi}^{I}=(h^{I},-h_{u}^{I},h_{d}^{I},h_{e}^{I},\phi_{u}^{I},\phi_{d}^{I},\phi_{e}^{I})^{T}$
where we take $-h_{u}^{I}$ for $\tilde{H}_{u}$ whose charge is the same as
that of $H_{d,e}$, the physical axion degree of freedom is written as

$\displaystyle a$ $\displaystyle\simeq$
$\displaystyle\frac{1}{\sqrt{\sum_{f}Q_{f}^{2}(v_{f}^{2}+v_{\Phi_{f}}^{2})}}\biggl{(}-\sum(-1)^{f}Q_{f}v_{f}^{2}/v,-Q_{u}v_{u},$
(7) $\displaystyle
Q_{d}v_{d},Q_{e}v_{e},Q_{u}v_{\Phi_{u}},Q_{d}v_{\Phi_{d}},Q_{e}v_{\Phi_{e}}\biggr{)}\,\vec{\Phi}^{I}$
$\displaystyle\equiv$ $\displaystyle\vec{G}_{\rm PQ}^{T}\vec{\Phi}^{I},$

where $v\,(\simeq 246\ {\rm
GeV}),v_{u},v_{d},v_{e},v_{\Phi_{u}},v_{\Phi_{d}},v_{\Phi_{e}}$ denote vacuum
expectation values (VEVs) of $H,H_{u},H_{d},H_{e},\Phi_{u},\Phi_{d}$ and
$\Phi_{e}$, respectively, and $(-1)^{f}=+1$ for $f=d,e$ and $(-1)^{f}=-1$ for
$f=u$. In estimating the prefactor, we have used $v_{f}\ll v$. The axion-
fermion coupling is then given by

$\displaystyle\mathcal{L}_{aff}$ $\displaystyle=$
$\displaystyle\sum_{f=u,d,e}\frac{m_{f}}{v_{a}}Q_{f}ia\bar{f}\gamma_{5}f$ (8)
$\displaystyle-$
$\displaystyle\frac{\sum_{f}(-1)^{f}Q_{f}v_{f}^{2}}{v^{2}}\sum_{F={\rm
2nd,3rd}}\frac{m_{F}}{v_{a}}ia\bar{F}\gamma_{5}F\,,$

with the PQ breaking scale,
$v_{a}\equiv\sqrt{\sum_{f}Q_{f}^{2}(v_{f}^{2}+v_{\Phi_{f}}^{2})}$. Note that
we have the axion couplings to the 2nd and 3rd generation fermions, but they
are suppressed by $(v_{f}/v)^{2}$. In Sec. V, we will point out extra
contributions to those couplings induced by gluon loops.

The five physical pseudoscalar bosons (other than the axion and the Nambu-
Goldstone mode absorbed by the $Z$ boson) need to be heavy enough to evade the
current experimental searches. In Ref. Liu  _et al._ (2021), the authors take
the magnitudes of the VEVs and the $A$-term coefficients as

$\displaystyle A_{4}\gg A_{k(=1,2,3)}\simeq 20\ {\rm GeV},$ $\displaystyle
v_{f(=u,d,e)}\simeq 20\ {\rm MeV},\quad v_{\Phi_{f(=u,d,e)}}\simeq 1\ {\rm
GeV}.$ (9)

The value of $v_{\Phi_{f}}$ sets the PQ breaking scale to $\mathcal{O}(1)\
{\rm GeV}$ ($m_{a}=\mathcal{O}(10)\ {\rm MeV}$). Then, the three pseudoscalars
are approximately given by linear combinations of $h_{f}^{I}$s, and their
masses are estimated as $A_{k}vv_{\Phi_{f}}/v_{f}\sim(300\ {\rm GeV})^{2}$.
Looking at the fourth term in Eq. (5), one massive mode is described as a
linear combination of $\Phi_{u}$ and $\Phi_{d}$, and its mass is given by
$A_{4}v_{\Phi_{f}}\gg(1\ {\rm GeV})^{2}$. We can take a large value of $A_{4}$
to make this mode sufficiently heavy. For $n=2$, $A_{5}$ is dimensionful, and
hence the last pseudoscalar mode can be heavy by taking a large value of
$A_{5}$. On the other hand, for $n=3$, $A_{5}$ is a dimensionless constant, so
that the last mode has a relatively small mass of
$A_{5}v_{\Phi_{f}}^{2}\sim(1\ {\rm GeV})^{2}$. This mode is dominantly
composed of $\Phi_{e}$ without other $A$-term interactions, and then mainly
couples to the electron. The coupling strength is estimated as
$|Y_{e}|v_{f}/v_{\Phi_{f}}\simeq m_{e}/v_{\Phi_{f}}\simeq 5\cdot 10^{-4}$, and
constrained by dark photon searches, such as BaBar Lees  _et al._ (2014) for
the mass range of $1-10\ {\rm GeV}$. Hence, the case of $n=3$ is marginally
consistent with the current bound. In addition to those pseudoscalars, there
are seven CP-even scalars in the model. First, we have the SM-like Higgs boson
with mass-squared of $\sim\lambda_{H}v^{2}$. Then, as in the case of the
pseudoscalar modes, three modes obtain masses of $\sim
A_{k}vv_{\Phi_{f}}/v_{f}$ and two modes have masses of $\sim
A_{4}v_{\Phi_{f}}$. For $n=2$, taking a large value of $A_{5}$, all the CP-
even scalars are sufficiently heavy. For $n=3$, there remains a light mode
whose mass is of $\mathcal{O}(1)\ {\rm GeV}$, and the similar experimental
constraint as that of the pseudoscalar is applied to this mode.333For $n=1$,
the model can be more minimal, composed of five Higgs fields. All the three
residual pseudoscalars (and CP-even scalars) are sufficiently heavy, because
their masses are given by $A_{k}vv_{\Phi_{f}}/v_{f}$ and $A_{4}v_{\Phi_{f}}$.
However, as it has been pointed out in Refs. Alves (2021); Liu  _et al._
(2021) that the reported $e^{+}e^{-}$ excess in the de-excitation of 8Be and
4He Krasznahorkay  _et al._ (2016, 2019) can be conservatively explained by
the current axion model for $n=2,3$, we mainly consider the case of $n=2,3$ in
the present paper. We assume the magnitudes of Eq. (9) in the following
discussion.

##  III Quality and Cosmology

The PQ mechanism can solve the strong CP problem dynamically under the
assumption that the PQ symmetry is well preserved. However, the PQ symmetry is
introduced as a global symmetry, and we naturally expect Planck-suppressed
operators breaking the PQ symmetry explicitly Dine and Seiberg (1986); Barr
and Seckel (1992); Kamionkowski and March-Russell (1992a, b); Holman  _et al._
(1992); Kallosh  _et al._ (1995); Carpenter  _et al._ (2009a, b). For
invisible axion models with large $f_{a}$, they destroy the PQ solution, which
is known as the axion quality problem. Here, we discuss whether the similar
problem exists or not in the visible axion model with
$f_{a}=\mathcal{O}(1)\,\rm GeV$. We also explore a possibility that those PQ
violating operators provide bias terms to solve the domain wall problem.

###  III.1 Explicit PQ breaking operators

For the visible QCD axion model presented in Sec. II, PQ-breaking higher-
dimensional operators that are the most dangerous for the PQ solution are

$\displaystyle\mathcal{L}_{\cancel{\rm
PQ}}=\sum_{f=u,d,e}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle\left(g_{1}\frac{(HH^{\dagger})^{2}\Phi_{f}}{M_{\rm
Pl}}+g_{2}\frac{HH^{\dagger}\Phi_{f}^{3}}{M_{\rm Pl}}\right.$ (10)
$\displaystyle+$ $\displaystyle g_{3}\left.\frac{\Phi_{f}^{5}}{M_{\rm
Pl}}+\cdots\right)+{\rm h.c.},$

where each $g_{i}=|g_{i}|e^{i\delta_{i}}$ $(i=1,2,3)$ is a complex coupling
constant with a phase shift $\delta_{i}$ and is assumed to be independent of
$f$ for simplicity. Here, we have omitted to write down combinations such as
$\Phi_{u}^{3}\Phi_{d}^{2}$, which will be taken into account in the final
results. The three terms in Eq. (10) respectively generate the following axion
potential terms:

$\displaystyle V^{(1)}_{\cancel{\rm PQ}}$ $\displaystyle=$
$\displaystyle-\frac{|g_{1}|}{2\sqrt{2}}\frac{v^{4}}{M_{\rm
Pl}}\sum_{f}v_{\Phi_{f}}\cos\left(Q_{f}\frac{a}{v_{a}}+\delta_{1}\right),$
(11) $\displaystyle V^{(2)}_{\cancel{\rm PQ}}$ $\displaystyle=$
$\displaystyle-\frac{|g_{2}|}{2\sqrt{2}}\frac{v^{2}}{M_{\rm
Pl}}\sum_{f}v_{\Phi_{f}}^{3}\cos\left(3Q_{f}\frac{a}{v_{a}}+\delta_{2}\right),$
(12) $\displaystyle V^{(3)}_{\cancel{\rm PQ}}$ $\displaystyle=$
$\displaystyle-\frac{|g_{3}|}{2\sqrt{2}}\sum_{f}\frac{v_{\Phi_{f}}^{5}}{M_{\rm
Pl}}\cos\left(5Q_{f}\frac{a}{v_{a}}+\delta_{3}\right).$ (13)

They shift the axion VEV from the CP conserving minimum, which is set by the
potential from QCD nonperturbative effects,
$V_{a}=m_{a}^{2}f_{a}^{2}(1-\cos(a/f_{a}))$. Here the axion mass $m_{a}$ is
given by

$\displaystyle
m_{a}=\frac{\sqrt{m_{u}m_{d}}}{m_{u}+m_{d}}\frac{m_{\pi}f_{\pi}}{f_{a}}\,,$
(14)

where the decay constant of the axion is defined as $f_{a}\equiv v_{a}/N_{\rm
DW}$ with domain wall number $N_{\rm DW}(=Q_{u}+Q_{d}=3)$, and $m_{\pi}\simeq
135\ {\rm MeV}$ and $f_{\pi}\simeq 92\ {\rm MeV}$ denote the pion mass and
decay constant, respectively. For the term $V_{\cancel{\rm PQ}}^{(1)}$, the
shift of the axion VEV is evaluated as

$\displaystyle\frac{|\langle
a\rangle|}{f_{a}}\simeq\frac{|g_{1}|}{2\sqrt{2}N_{\rm DW}}\frac{v^{4}}{M_{\rm
Pl}m_{a}^{2}f_{a}^{2}}\sum_{f}Q_{f}v_{\Phi_{f}}\sin\delta_{1}.$ (15)

Those corresponding to $V_{\cancel{\rm PQ}}^{(2)}$ and $V_{\cancel{\rm
PQ}}^{(3)}$ can be estimated similarly.

The current neutron EDM measurement provides an upper bound on $\bar{\theta}$
Baker  _et al._ (2006); Pendlebury  _et al._ (2015), such that $|\langle
a\rangle|/f_{a}\lesssim 10^{-10}$. This bound can be read out as upper bounds
on the coupling coefficients $|g_{i}|$. Considering each operator as the
dominant one, we obtain

$\displaystyle|g_{1}|$ $\displaystyle\lesssim$ $\displaystyle 8\times
10^{-6}\left(\frac{m_{a}}{10\ {\rm MeV}}\right),$ (16) $\displaystyle|g_{2}|$
$\displaystyle\lesssim$ $\displaystyle
0.03\left(\frac{\Delta^{(2)}}{10}\right)^{-1}\left(\frac{3}{N_{\rm
DW}}\right)^{2}\left(\frac{m_{a}}{10\ {\rm MeV}}\right)^{3},$ (17)
$\displaystyle|g_{3}|$ $\displaystyle\lesssim$ $\displaystyle
200\left(\frac{\Delta^{(3)}}{100}\right)^{-1}\left(\frac{3}{N_{\rm
DW}}\right)^{4}\left(\frac{m_{a}}{10\ {\rm MeV}}\right)^{5}.$ (18)

Here, we use $f_{a}\simeq\sqrt{\sum Q_{f}^{2}v_{\Phi_{f}}^{2}}/N_{\rm DW}$ and
take $n=2$ and $\delta_{i}=1$. Since we expect $\mathcal{O}(10)$ and
$\mathcal{O}(100)$ number of operators in Eq. (10) giving similar
contributions to Eq. (12) and Eq. (13), respectively, the factors of
$\Delta^{(2)}=\mathcal{O}(10)$ and $\Delta^{(3)}=\mathcal{O}(100)$ are
included for the bounds on $|g_{2}|$ and $|g_{3}|$. The results show that
fine-tuning of $\mathcal{O}(10^{-5})$ is required for the PQ quality in the
presence of the coupling $g_{1}$. The $f_{a}$ (or $m_{a}$) dependence of
bounds on $|g_{i}|$ is also summarized in Fig. 1. The results for $n=3$ are
not significantly different from those of $n=2$.

While, compared to the invisible QCD axion models, the PQ quality for the
visible axion at the 10 MeV scale is much better, the complete PQ quality
requires removal of higher-dimensional operators involving the SM Higgs. This
conclusion comes from the fact that the PQ breaking scale is below the
electroweak scale. A small axion decay constant is not sufficient to achieve
the high quality of the PQ symmetry. However, this issue might be addressed by
a further model-building.

Figure 1:  Left panel : upper bounds on the coupling constants $|g_{1,2}|$
from the neutron EDM measurement as functions of $f_{a}$ (lower axis) or
$m_{a}$ (upper axis) for the case of $n=2$. The red solid and blue dashed
lines represent the cases of $\delta_{1,2}=1,10^{-2},10^{-4}$ for Eq. (11) and
Eq. (12), respectively. Right panel : bounds on the coupling constant
$|g_{3}|$ of Eq. (13) from the neutron EDM measurement for $n=2$. The red
solid and dotted lines denote the bounds for $\delta_{3}=1,10^{-2},10^{-4}$.
In the blue shaded region, the domain wall survives after BBN. In both panels,
the orange shaded regions are excluded by beam dump experiments Riordan  _et
al._ (1987).

###  III.2 Domain wall problem

Spontaneous breaking of the PQ symmetry after inflation generates topological
defects which can give a non-negligible impact on the evolution of the
Universe. At the PQ phase transition, cosmic strings are formed through the
Kibble mechanism Kibble (1976); Kibble  _et al._ (1982), and when the axion
acquires a potential from nonperturbative effects of QCD, domain walls are
produced and the string-wall system appears. Its cosmological consequence
depends on the domain wall number $N_{\rm DW}$ of a model. If $N_{\rm DW}=1$,
the string-wall system is unstable and collapses with a short lifetime. On the
other hand, if $N_{\rm DW}>1$, it is stable and may quickly dominate the
Universe. This is known as the cosmological domain wall problem (see Refs.
Beyer and Sarkar (2023); Dine (2023) for recent reviews). One possible
solution to the domain wall problem is to introduce an explicit PQ breaking
operator (often called bias term), which lifts up the degenerate vacua and
makes the string-wall system unstable Vilenkin (1981); Gelmini  _et al._
(1989); Larsson  _et al._ (1997). However, for invisible axion models, the
solution is in tension with the quality of the PQ symmetry, and one has to
adjust the axion potential from the explicit breaking operator so that its
minimum is aligned with that of the potential from nonperturbative effects of
QCD Kawasaki  _et al._ (2015).

The situation is improved considerably for the visible axion model with
$f_{a}=\mathcal{O}(1)\,\rm GeV$. In order to solve the domain wall problem, it
is sufficient to find a condition that the domain wall collapses through a PQ-
breaking operator before the Big Bang Nucleosynthesis (BBN) so as not to
destroy the light elements. The PQ-breaking operator makes an energy
difference $\Delta V$ between the vacua, which induces a volume pressure on
the domain wall. The wall becomes unstable when the volume pressure is
comparable to the domain wall tension:

$\displaystyle\Delta V\sim\frac{\beta\sigma_{\rm wall}}{t_{\rm dec}}\,.$ (19)

Here, $\sigma_{\rm wall}\simeq 9.25m_{a}f_{a}^{2}$ is the surface density of
the domain wall Huang and Sikivie (1985); Hiramatsu  _et al._ (2013), $\beta$
is a constant which is $\mathcal{O}(1)$ by dimensional analysis, and $t_{\rm
dec}$ is the cosmic time when the domain wall decays. The potential difference
between the origin and its neighbor vacuum is approximately given by

$\displaystyle\Delta V\sim|V_{\cancel{\rm PQ}}(a=2\pi v_{a}/N_{\rm
DW})-V_{\cancel{\rm PQ}}(0)|\,.$ (20)

For $V_{\cancel{\rm PQ}}^{(3)}$ in Eq. (13), the constraint from the domain
wall lifetime is shown in the blue shaded region of the right panel of Fig. 1.
Although we take $n=2$ here, a similar result can be obtained for $n=3$. We
can see that the constraint is still satisfied without causing the axion
quality problem. For $V_{\cancel{\rm PQ}}^{(1,2)}$ in Eqs. (11), (12), the
domain wall decays even faster, and the constraint does not appear in the left
panel of Fig. 1 although there is a quality issue for those operators. In
conclusion, if we can somehow suppress higher-dimensional operators involving
the SM Higgs which lead to $V_{\cancel{\rm PQ}}^{(1,2)}$ but allow the third
term of Eq. (10) leading to Eq. (13) by model-building, the domain wall
problem is solved with the high-quality of the PQ symmetry.

##  IV Electron EDM

(a) One-loop diagram in the case of $n=2$

(b) One-loop diagram in the case of $n=3$

(c) Two-loop Barr-Zee diagram with a $W$-loop

Figure 2: Feynman diagrams for dominant contributions to the electron EDM. The
arrows in the diagram represent charge flow.

The model contains three Higgs doublet fields $H_{f}$ ($f=u,d,e$) and SM-
singlet scalars $\Phi_{f}$ as well as the SM Higgs field. According to the
charge assignment summarized in Tab. 1, many scalar potential terms other than
those presented in Eqs. (5), (6) are actually allowed by the SM gauge and PQ
global symmetries. In the case of $n=2$, all the possible (off-diagonal)
scalar potential terms are

$\displaystyle V^{(n=2)}_{\rm scalar}$
$\displaystyle=A_{1}HH_{u}\Phi_{u}^{*}+A_{2}HH_{d}^{\dagger}\Phi_{d}+A_{3}H_{e}H^{\dagger}\Phi_{e}^{*}$
$\displaystyle+A_{4}\Phi_{u}^{*}\Phi_{d}^{2}+A_{5}\Phi_{d}^{*}\Phi_{e}^{2}+A_{6}H_{d}H_{e}^{\dagger}\Phi_{e}^{*}$
$\displaystyle+B_{1}HH_{u}{\Phi_{d}^{*}}^{2}+B_{2}HH_{d}^{\dagger}\Phi_{u}\Phi_{d}^{*}+B_{3}HH_{d}^{\dagger}\Phi_{e}^{2}$
$\displaystyle+B_{4}H^{\dagger}H_{e}\Phi_{d}^{*}\Phi_{e}+B_{5}H_{u}H_{e}\Phi_{u}^{*}\Phi_{e}^{*}+B_{6}H_{u}H_{d}\Phi_{u}^{*}\Phi_{d}^{*}$
$\displaystyle+B_{7}H_{d}^{\dagger}H_{e}\Phi_{d}\Phi_{e}^{*}+B_{8}\Phi_{u}\Phi_{d}^{*}{\Phi_{e}^{*}}^{2}+B_{9}H_{e}^{2}H^{\dagger}H_{d}^{\dagger}$
$\displaystyle+\rm{h.c.}\,,$ (21)

where $A_{j}$ ($j=1,2,\cdots,6$) and $B_{k}$ ($k=1,2,\cdots,9$) are parameters
with mass dimension one and dimensionless coupling constants, respectively.
All the Higgs doublets and singlets get VEVs. Those parameters and VEVs are
generally complex. Field redefinition rotates away some of the phases but is
not sufficient to remove all, leading to observable phases.

The physical complex phases violate CP symmetry and lead to EDMs. Fig. 2(a)
describes a one-loop diagram generating the electron EDM. The diagram includes
the interactions, $A_{3}^{*}H_{e}^{\dagger}H\Phi_{e}$ and
$B_{4}^{*}HH_{e}^{\dagger}\Phi_{d}\Phi_{e}^{*}$. The contribution to the
electron EDM is then estimated as

$\displaystyle\frac{d_{e}}{e}\sim\frac{\theta}{16\pi^{2}}|A_{3}^{*}B_{4}^{*}v^{2}v_{\Phi_{d}}|\frac{m_{e}}{M^{6}}\,.$
(22)

Here, $\theta$ and $M$ denote the relevant physical phase and a representative
mass scale of particles in the loop, respectively. The current experimental
upper limit on the electron EDM is $d_{e}\leq 4.1\times 10^{-30}e\,\rm cm$
Roussy  _et al._ (2023). With our choice of Eq. (9), the constraint can be
satisfied for the mass scale,

$\displaystyle M\gtrsim 500\,\rm GeV\,,$ (23)

where we have assumed the phase $\theta$ and the dimensionless coupling
$B_{4}$ are of order one. As we have seen in Sec. II, all the (pseudo-)scalar
particles inside the loop can have masses of $\mathcal{O}(100)\,\rm GeV$ in
the case of $n=2$. Therefore, the condition (23) is marginally satisfied.

In the case of $n=3$, the possible (off-diagonal) scalar potential terms
allowed by the symmetries are

$\displaystyle V^{(n=3)}_{\rm scalar}$
$\displaystyle=A_{1}HH_{u}\Phi_{u}^{*}+A_{2}HH_{d}^{\dagger}\Phi_{d}+A_{3}H_{e}H^{\dagger}\Phi_{e}^{*}$
$\displaystyle+A_{4}\Phi_{u}^{*}\Phi_{d}^{2}+A^{\prime}_{5}\Phi_{d}^{*}\Phi_{e}^{3}+B_{1}HH_{u}{\Phi_{d}^{*}}^{2}$
$\displaystyle+B_{2}HH_{d}^{\dagger}\Phi_{u}\Phi_{d}^{*}+B_{5}H_{u}H_{e}\Phi_{u}^{*}\Phi_{e}^{*}+B_{6}H_{u}H_{d}\Phi_{u}^{*}\Phi_{d}^{*}$
$\displaystyle+B_{7}H_{d}^{\dagger}H_{e}\Phi_{d}\Phi_{e}^{*}+B_{10}H_{d}H_{e}^{\dagger}{\Phi_{e}^{*}}^{2}+\rm{h.c.}\,,$
(24)

where $A^{\prime}_{5}$ and $B_{10}$ denote dimensionless coupling constants.
The diagram shown in Fig. 2(b) contains the interactions,
$A_{3}^{*}H_{e}^{\dagger}H\Phi_{e}$ and
$B_{10}H_{d}H_{e}^{\dagger}{\Phi_{e}^{*}}^{2}$. The contribution to the
electron EDM is then estimated as

$\displaystyle\frac{d_{e}}{e}\sim\frac{\delta}{16\pi^{2}}|A_{3}^{*}B_{10}vv_{d}v_{\Phi_{e}}^{*}|\frac{m_{e}}{M^{6}}\,,$
(25)

with a relevant physical phase $\delta$. The experimental upper limit can be
satisfied for

$\displaystyle M\gtrsim 100\,\rm GeV\,.$ (26)

Compared to the case of $n=2$, the lower bound on the mass scale $M$ becomes
weaker. However, as described in Sec. II, the mass of $\Phi_{e}$ in the loop
is typically $\mathcal{O}(1)\,\rm GeV$ in the case of $n=3$. Hence, it is
likely that some fine-tuning in CP-violating phases is required to satisfy the
EDM constraint. A more precise numerical estimate for the amount of tuning is
left for a future study.

In addition to the one-loop contributions, a two-loop diagram with a fewer
number of small VEV insertions may lead to a sizable contribution to the
electron EDM. Fig. 2(c) describes a two-loop Bar-Zee-type diagram with inner
$W$-loop. The diagram involves the CP-violating interaction
$A_{3}H_{e}H^{\dagger}\Phi_{e}^{*}$, and its contribution to the electron EDM
is estimated as Barr and Zee (1990) (the general formula for a two-loop Bar-
Zee-type diagram has been presented in Ref. Nakai and Reece (2017) (see also
Ref. Cesarotti  _et al._ (2019)))

$\displaystyle\frac{d_{e}}{e}\sim\frac{\alpha_{e}\zeta}{(4\pi)^{3}}\sqrt{2}G_{F}$
$\displaystyle
m_{e}\left(3f\left(\frac{m_{W}^{2}}{m_{H}^{2}}\right)+5g\left(\frac{m_{W}^{2}}{m_{H}^{2}}\right)\right)$
$\displaystyle\times\left|\frac{A_{3}v_{\Phi_{e}}^{*}}{2\sqrt{2}M^{2}}\right|\left(\sin^{2}\beta\tan\beta\right)\,,$
(27)

where $\alpha_{e}$ and $G_{F}$ denote the fine-structure constant and Fermi
coupling constant, respectively, $\zeta$ is the relevant physical phase and

$\displaystyle f(z)$
$\displaystyle\equiv\frac{1}{2}z\int_{0}^{1}dx\frac{1-2x(1-x)}{x(1-x)-z}\ln\left(\frac{x(1-x)}{z}\right)\,,$
$\displaystyle g(z)$
$\displaystyle\equiv\frac{1}{2}z\int_{0}^{1}dx\frac{1}{x(1-x)-z}\ln\left(\frac{x(1-x)}{z}\right)\,.$

The parameter $\beta$ is defined as

$\displaystyle\tan\beta\equiv\frac{v}{v_{e}}\sim 10^{4}\,.$ (28)

The experimental constraint can be satisfied for

$\displaystyle M\gtrsim 8\,\rm TeV\,,$ (29)

which gives a stronger bound (independent of the value of $Q_{e}$), compared
to those of the one-loop contributions. Fine-tuning in CP-violating phases
and/or some mechanism to suppress the phases are necessary.

In addition to the electron EDM constraint, we are able to consider those of
the neutron EDM Baker  _et al._ (2006) and mercury EDM Graner  _et al._ (2016)
where the chromo-EDM of the down quark and the electron EDM give main
contributions. The resulting constraints are found to be weaker than the one
directly given by the electron EDM measurement.

##  V Loop induced B-decay

Figure 3: Feynman diagrams that generate the axion couplings to gluons and the
top quark, mediating the $B$-meson decay.

To evade stringent constraints from rare decay modes of heavy mesons such as
$B$, $\Upsilon$ mediated by the QCD axion variant under consideration, it has
been assumed that the axion only couples to the first-generation quarks $u$,
$d$, and to the electron at the tree level. However, from a UV completion
perspective, the PQ symmetry breaking, together with the electroweak symmetry
breaking, has to generate the observed CKM mixing matrix. Therefore,
inevitably, couplings of the axion to heavier generation quarks will be
generated via mixing with the SM Higgs. Progress towards understanding this
has been made in Ref. Liu  _et al._ (2021), where a hierarchy between the VEVs
of additional Higgses and the SM Higgs has been deemed to suppress the axion
couplings to the 2nd and 3rd generation quarks to $\sim{\cal{O}}(10^{-7})$
level.

Here we point out that a gluon loop will also induce the axion couplings to
the heavier quarks. Namely,

$\displaystyle{\cal L}_{aFF}$ $\displaystyle=Q_{F,\rm eff}^{\rm
PQ}\frac{m_{F}}{v_{a}}a\bar{F}i\gamma_{5}F\ ,$ $\displaystyle Q_{F,\rm
eff}^{\rm PQ}$ $\displaystyle\simeq
C_{F}\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\ ,$ (30)

where $F$ can be any 2nd or 3rd generation quark, and $C_{F}=4/3$ is the
quadratic Casimir, and $\alpha_{s}$ is the strong coupling constant. The two-
loop suppression effect in Eq. (30) comes from the fact that the axion-gluon
coupling inherits a loop factor, generated through the up and down quarks
flowing in the triangle anomaly diagrams. Furthermore, the anomalous dimension
of this coupling is $-6C_{F}(\alpha_{s}/4\pi)^{2}$, and it is not changed
significantly upon running. Therefore, the 1-loop gluon-induced coupling to
the bottom quark is $Q_{b,\rm eff}^{\rm PQ}\simeq 3\times 10^{-4}$, where we
have used the running coupling $\alpha_{s}$ relevant for the process
$\Upsilon\to\gamma a$, which requires $Q_{b,\rm eff}^{\rm PQ}\lesssim
0.8\times 10^{-2}$ Alves and Weiner (2018). Hence, the 1-loop gluon-generated
coupling to the $b$ quark is consistent with bounds from rare $\Upsilon$
decays.

Let us now consider the loop-mediated decays of the $B$-meson such as $B\to
Ka$. This is obtained by integrating out the $W$-loop as in Fig. 3. The
effective $b$-$s$-$a$ vertex is denoted as

${\cal L}_{bsa}=-ig_{bsa}\bar{s}_{L}b_{R}a+{\rm h.c.}\ ,$ (31)

where we utilized the chirality-flipping nature of the diagram to infer that
the coefficient of $\bar{s}_{R}b_{L}$ is suppressed as $m_{s}/m_{b}$. The
leading-log result is well-known Batell  _et al._ (2011),

$g_{bsa}=\frac{G_{F}m_{W}^{2}}{4\sqrt{2}\pi^{2}}\frac{m_{b}}{v_{a}}\sum_{u_{i}=u,c,t}Q_{u_{i},{\rm{eff}}}^{\rm
PQ}\frac{m_{u_{i}}^{2}}{m_{W}^{2}}V_{u_{i}s}^{*}V_{u_{i}b}\ln{\left(\frac{\Lambda_{\rm
UV}^{2}}{m_{u_{i}}^{2}}\right)}\ ,$ (32)

where we have included a summation over all the up-type quarks, with an
effective coupling to the higher generation understood as given in Eq. (30),
and $\Lambda_{\rm UV}$ is taken as the UV cut-off of the theory. In order to
make the leading log expansion sensible, the log factor is assumed to be close
to unity. Although the top-axion coupling is loop generated, it contributes
dominantly to $g_{bsa}$ as evident in Eq. (32). To be consistent with the beam
dump decay experiment, the axion should decay into an $e^{+}e^{-}$ pair with
lifetime $\tau_{a}\lesssim 10^{-13}$ s. Hence, we employ the experimental
measurement on the branching ratio of the rare decay mode $B^{0}\to
K^{*0}e^{+}e^{-}$ by LHCb Aaij  _et al._ (2013),444 Here the errors are
respectively statistical, systematic, and relating to the uncertainties on the
$B^{0}\to J/\Psi K^{*0}$ and $J/\Psi\to e^{+}e^{-}$ branching fractions, which
were used as normalization channels.

${\cal B}(B^{0}\to K^{*0}e^{+}e^{-})=\left(3.1^{+0.9+0.2}_{-0.8-0.3}\pm
0.2({\cal B})\right)\times 10^{-7}\ .$ (33)

Note that this measurement is only applicable for dilepton invariant mass in
the range of $30$-$1000$ MeV, which is slightly higher than the axion mass we
are interested in, however, it should be noted that the theoretical
expectation from the SM is not very sensitive to the dilepton invariant mass,
and can be predicted up to the order one logarithmic correction factor as
Jäger and Martin Camalich (2013)

${\cal B}^{\rm SM}(B^{0}\to K^{*0}e^{+}e^{-})\simeq 2.43^{+0.66}_{-0.47}\times
10^{-7}\ .$ (34)

Given the axion coupling in Eq. (31), the BSM contribution can be written as

$\displaystyle{\cal B}^{\rm BSM}$ $\displaystyle(B^{0}\to K^{*0}e^{+}e^{-})$
$\displaystyle=\frac{|g_{bsa}|^{2}}{16\pi
m_{B}\Gamma_{B}}\lambda_{B,K^{*}a}^{1/2}|\langle
K^{*}|\bar{s}_{L}b_{R}|B\rangle|^{2}\ ,$ (35)

where

$\displaystyle\lambda_{x,yz}\equiv\left[1-\left(\frac{m_{y}}{m_{x}}-\frac{m_{z}}{m_{x}}\right)^{2}\right]\left[1-\left(\frac{m_{y}}{m_{x}}+\frac{m_{z}}{m_{x}}\right)^{2}\right]\
,$ $\displaystyle|\langle
K^{*}|\bar{s}_{L}b_{R}|B\rangle|^{2}=\frac{1}{4}\frac{m_{B}^{4}\lambda_{B,K^{*}a}}{(m_{b}+m_{s})^{2}}[{\cal
A}_{K^{*}}(m_{a}^{2})]^{2}\ ,$ (36) $\displaystyle{\cal
A}_{K^{*}}(q^{2})=\frac{1.36}{1-q^{2}/27.9\,{\rm
GeV}^{2}}-\frac{0.99}{1-q^{2}/36.8\,{\rm GeV}^{2}}\ .$

The pion-phobic QCD axion variant, which by an ad-hoc model construction only
couples to the first-generation quarks at the tree level, when considered
together with the gluon loop generated coupling to the top quark, can mediate
the $B$-decay with a contribution much larger than that of the SM. For
numerical illustration, taking $m_{a}=10$ MeV, we find

$\frac{{\cal B}^{\rm BSM}(B^{0}\to K^{*0}a(\to e^{+}e^{-}))}{{\cal B}^{\rm
SM}(B^{0}\to K^{*0}e^{+}e^{-})}\simeq 10^{3}\ln{\left(\frac{\Lambda_{\rm
UV}}{m_{t}}\right)}^{2}\ ,$ (37)

where we have used the QCD relation in Eq. (14) together with $v_{a}=N_{\rm
DW}f_{a}$ for this model, which fixes $v_{a}\simeq 1.7$ GeV for $m_{a}=10$
MeV. This BSM contribution to the rare meson decay increases with increasing
$m_{a}$, as the corresponding $v_{a}$ decreases. Hence, in light of the LHCb
analysis Aaij  _et al._ (2013), the pion-phobic QCD axion variant is ruled out
for $m_{a}>30$ MeV mass range.

The current most precise measurement of the branching ratio for $B\to
K^{(*)}e^{+}e^{-}$ comes from Belle Wei  _et al._ (2009) and BaBar Aubert  _et
al._ (2009) collaborations. However, they used a veto of $q^{2}\lesssim 0.1$
GeV${}^{2}/c^{4}$ to reduce the background. Therefore, we cannot use their
result to constrain the axion mass below $\sim$ 300 MeV. We note that below 30
MeV dilepton invariant mass, multiple scatterings obscure the angular analysis
of the outgoing electrons in $B^{0}\to K^{*0}e^{+}e^{-}$, however, if one is
only interested in searching for this mode without the angular correlation,
the background contamination only comes from the photon converting to
$e^{+}e^{-}$ in the mode $B^{0}\to K^{*0}\gamma$. Hence, there is a strong
motivation for the LHCb and Belle II experiments to make an updated analysis
and push the search for this mode towards the lower end of the dilepton
invariant mass from 30 MeV down to a few MeV that will determine the fate of
the model conclusively.

##  VI Conclusions

We have considered new theoretical and phenomenological constraints on the
${\cal O}$(10) MeV pion-phobic QCD axion variant that has been argued to be
still viable Alves and Weiner (2018). Due to the very low PQ breaking scale,
the phase transition takes place after inflation, and the string-domain wall
network is formed. We have investigated the interplay of the lifetime of the
network and the quality issue of the PQ symmetry due to a necessary bias term.
It was found that an appropriate choice of a higher-dimensional operator by
some model-building can solve the domain wall problem with the high-quality of
the PQ symmetry. From a UV completion perspective, as the PQ charges of the up
and down quarks and electron have to be different in the model, one has to
introduce multiple new scalar fields. This brings the model a plethora of new
CP-violating terms that manifest themselves as EDMs. We have identified the
dominant electron EDM contribution and analyzed the resulting constraints.
Finally, we have pointed out that the gluon loop-induced top-quark coupling of
the axion contributes to rare $B$-decay modes, and in particular to $B\to
K^{(*)}a(\to e^{+}e^{-})$ for the current axion model. Experimental data for
this mode already exclude the model for the axion mass larger than 30 MeV.
There is a strong motivation for pushing the experimental analysis to a lower
dilepton mass window by the LHCb and Belle collaborations, which will dictate
the fate of the model.

## Acknowledgments

We thank Chia-Wei Liu and Yoshihiro Shigekami for useful discussions.

## References

  * Ellis and Gaillard (1979) J. R. Ellis and M. K. Gaillard, Nucl. Phys. B 150, 141 (1979).
  * Workman  _et al._ (2022) R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Alexandrou  _et al._ (2020) C. Alexandrou, J. Finkenrath, L. Funcke, K. Jansen, B. Kostrzewa, F. Pittler, and C. Urbach, Phys. Rev. Lett. 125, 232001 (2020), arXiv:2002.07802 [hep-lat] .
  * Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Kim (1979) J. E. Kim, Phys. Rev. Lett. 43, 103 (1979).
  * Dine  _et al._ (1981) M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B 104, 199 (1981).
  * Alves and Weiner (2018) D. S. M. Alves and N. Weiner, JHEP 07, 092 (2018), arXiv:1710.03764 [hep-ph] .
  * Krasznahorkay  _et al._ (2016) A. J. Krasznahorkay _et al._ , Phys. Rev. Lett. 116, 042501 (2016), arXiv:1504.01527 [nucl-ex] .
  * Liu  _et al._ (2021) J. Liu, N. McGinnis, C. E. M. Wagner, and X.-P. Wang, JHEP 05, 138 (2021), arXiv:2102.10118 [hep-ph] .
  * Krauss and Nash (1988) L. M. Krauss and D. J. Nash, Phys. Lett. B 202, 560 (1988).
  * Alves (2021) D. S. M. Alves, Phys. Rev. D 103, 055018 (2021), arXiv:2009.05578 [hep-ph] .
  * Aaij  _et al._ (2013) R. Aaij _et al._ (LHCb), JHEP 05, 159 (2013), arXiv:1304.3035 [hep-ex] .
  * Donnelly  _et al._ (1978) T. W. Donnelly, S. J. Freedman, R. S. Lytel, R. D. Peccei, and M. Schwartz, Phys. Rev. D 18, 1607 (1978).
  * Wilczek (1977) F. Wilczek, Phys. Rev. Lett. 39, 1304 (1977).
  * Hall and Wise (1981) L. J. Hall and M. B. Wise, Nucl. Phys. B 187, 397 (1981).
  * Albrecht  _et al._ (1986) H. Albrecht _et al._ (ARGUS), Phys. Lett. B 179, 403 (1986).
  * Hsueh and Palestini (1992) S. Y. Hsueh and S. Palestini, Phys. Rev. D 45, R2181 (1992).
  * Armstrong  _et al._ (1993) T. A. Armstrong _et al._ (E760), Phys. Rev. D 47, 772 (1993).
  * Bai  _et al._ (1995) J. Z. Bai _et al._ (BES), Phys. Lett. B 355, 374 (1995), [Erratum: Phys.Lett.B 363, 267 (1995)].
  * Adams  _et al._ (2006) G. S. Adams _et al._ (CLEO), Phys. Rev. D 73, 051103 (2006), arXiv:hep-ex/0512046 .
  * Bergsma  _et al._ (1985) F. Bergsma _et al._ (CHARM), Phys. Lett. B 157, 458 (1985).
  * Bjorken  _et al._ (1988) J. D. Bjorken, S. Ecklund, W. R. Nelson, A. Abashian, C. Church, B. Lu, L. W. Mo, T. A. Nunamaker, and P. Rassmann, Phys. Rev. D 38, 3375 (1988).
  * Davier and Nguyen Ngoc (1989) M. Davier and H. Nguyen Ngoc, Phys. Lett. B 229, 150 (1989).
  * Blumlein  _et al._ (1991) J. Blumlein _et al._ , Z. Phys. C 51, 341 (1991).
  * Eichler  _et al._ (1986) R. Eichler _et al._ (SINDRUM), Phys. Lett. B 175, 101 (1986).
  * Parker  _et al._ (2018) R. H. Parker, C. Yu, W. Zhong, B. Estey, and H. Müller, Science 360, 191 (2018), arXiv:1812.04130 [physics.atom-ph] .
  * Morel  _et al._ (2020) L. Morel, Z. Yao, P. Cladé, and S. Guellati-Khélifa, Nature 588, 61 (2020).
  * Fan  _et al._ (2023) X. Fan, T. G. Myers, B. A. D. Sukra, and G. Gabrielse, Phys. Rev. Lett. 130, 071801 (2023), arXiv:2209.13084 [physics.atom-ph] .
  * Anastasi  _et al._ (2015) A. Anastasi _et al._ , Phys. Lett. B 750, 633 (2015), arXiv:1509.00740 [hep-ex] .
  * Banerjee  _et al._ (2020) D. Banerjee _et al._ (NA64), Phys. Rev. D 101, 071101 (2020), arXiv:1912.11389 [hep-ex] .
  * Depero  _et al._ (2020) E. Depero _et al._ (NA64), Eur. Phys. J. C 80, 1159 (2020), arXiv:2009.02756 [hep-ex] .
  * Lees  _et al._ (2014) J. P. Lees _et al._ (BaBar), Phys. Rev. Lett. 113, 201801 (2014), arXiv:1406.2980 [hep-ex] .
  * Krasznahorkay  _et al._ (2019) A. J. Krasznahorkay _et al._ , (2019), arXiv:1910.10459 [nucl-ex] .
  * Dine and Seiberg (1986) M. Dine and N. Seiberg, Nucl. Phys. B 273, 109 (1986).
  * Barr and Seckel (1992) S. M. Barr and D. Seckel, Phys. Rev. D 46, 539 (1992).
  * Kamionkowski and March-Russell (1992a) M. Kamionkowski and J. March-Russell, Phys. Lett. B 282, 137 (1992a), arXiv:hep-th/9202003 .
  * Kamionkowski and March-Russell (1992b) M. Kamionkowski and J. March-Russell, Phys. Rev. Lett. 69, 1485 (1992b), arXiv:hep-th/9201063 .
  * Holman  _et al._ (1992) R. Holman, S. D. H. Hsu, T. W. Kephart, E. W. Kolb, R. Watkins, and L. M. Widrow, Phys. Lett. B 282, 132 (1992), arXiv:hep-ph/9203206 .
  * Kallosh  _et al._ (1995) R. Kallosh, A. D. Linde, D. A. Linde, and L. Susskind, Phys. Rev. D 52, 912 (1995), arXiv:hep-th/9502069 .
  * Carpenter  _et al._ (2009a) L. M. Carpenter, M. Dine, and G. Festuccia, Phys. Rev. D 80, 125017 (2009a), arXiv:0906.1273 [hep-th] .
  * Carpenter  _et al._ (2009b) L. M. Carpenter, M. Dine, G. Festuccia, and L. Ubaldi, Phys. Rev. D 80, 125023 (2009b), arXiv:0906.5015 [hep-th] .
  * Baker  _et al._ (2006) C. A. Baker _et al._ , Phys. Rev. Lett. 97, 131801 (2006), arXiv:hep-ex/0602020 .
  * Pendlebury  _et al._ (2015) J. M. Pendlebury _et al._ , Phys. Rev. D 92, 092003 (2015), arXiv:1509.04411 [hep-ex] .
  * Riordan  _et al._ (1987) E. M. Riordan _et al._ , Phys. Rev. Lett. 59, 755 (1987).
  * Kibble (1976) T. W. B. Kibble, J. Phys. A 9, 1387 (1976).
  * Kibble  _et al._ (1982) T. W. B. Kibble, G. Lazarides, and Q. Shafi, Phys. Rev. D 26, 435 (1982).
  * Beyer and Sarkar (2023) K. A. Beyer and S. Sarkar, SciPost Phys. 15, 003 (2023), arXiv:2211.14635 [hep-ph] .
  * Dine (2023) M. Dine, (2023), arXiv:2307.04710 [hep-ph] .
  * Vilenkin (1981) A. Vilenkin, Phys. Rev. D 23, 852 (1981).
  * Gelmini  _et al._ (1989) G. B. Gelmini, M. Gleiser, and E. W. Kolb, Phys. Rev. D 39, 1558 (1989).
  * Larsson  _et al._ (1997) S. E. Larsson, S. Sarkar, and P. L. White, Phys. Rev. D 55, 5129 (1997), arXiv:hep-ph/9608319 .
  * Kawasaki  _et al._ (2015) M. Kawasaki, K. Saikawa, and T. Sekiguchi, Phys. Rev. D 91, 065014 (2015), arXiv:1412.0789 [hep-ph] .
  * Huang and Sikivie (1985) M. C. Huang and P. Sikivie, Phys. Rev. D 32, 1560 (1985).
  * Hiramatsu  _et al._ (2013) T. Hiramatsu, M. Kawasaki, K. Saikawa, and T. Sekiguchi, JCAP 01, 001 (2013), arXiv:1207.3166 [hep-ph] .
  * Roussy  _et al._ (2023) T. S. Roussy, L. Caldwell, T. Wright, W. B. Cairncross, Y. Shagam, K. B. Ng, N. Schlossberger, S. Y. Park, A. Wang, J. Ye, and E. A. Cornell, Science 381, 46–50 (2023).
  * Barr and Zee (1990) S. M. Barr and A. Zee, Phys. Rev. Lett. 65, 21 (1990), [Erratum: Phys.Rev.Lett. 65, 2920 (1990)].
  * Nakai and Reece (2017) Y. Nakai and M. Reece, JHEP 08, 031 (2017), arXiv:1612.08090 [hep-ph] .
  * Cesarotti  _et al._ (2019) C. Cesarotti, Q. Lu, Y. Nakai, A. Parikh, and M. Reece, JHEP 05, 059 (2019), arXiv:1810.07736 [hep-ph] .
  * Graner  _et al._ (2016) B. Graner, Y. Chen, E. Lindahl, and B. Heckel, Physical Review Letters 116 (2016), 10.1103/physrevlett.116.161601.
  * Batell  _et al._ (2011) B. Batell, M. Pospelov, and A. Ritz, Phys. Rev. D 83, 054005 (2011), arXiv:0911.4938 [hep-ph] .
  * Jäger and Martin Camalich (2013) S. Jäger and J. Martin Camalich, JHEP 05, 043 (2013), arXiv:1212.2263 [hep-ph] .
  * Wei  _et al._ (2009) J. T. Wei _et al._ (Belle), Phys. Rev. Lett. 103, 171801 (2009), arXiv:0904.0770 [hep-ex] .
  * Aubert  _et al._ (2009) B. Aubert _et al._ (BaBar), Phys. Rev. Lett. 102, 091803 (2009), arXiv:0807.4119 [hep-ex] .