aainstitutetext:  Department of Physics, Nagoya University, Furo-cho Chikusa-
ku, Nagoya 464-8602 Japan bbinstitutetext:  Kobayashi-Maskawa Institute for
the Origin of Particles and the Universe, Nagoya University, Furo-cho Chikusa-
ku, Nagoya 464-8602 Japan ccinstitutetext:  Kavli IPMU (WPI), UTIAS, The
University of Tokyo, Kashiwa 277-8584, Japan ddinstitutetext:  Institute for
Advanced Research, Nagoya University, Furo-cho Chikusa-ku, Nagoya 464-8601,
Japaneeinstitutetext:  KEK Theory Center, IPNS, KEK, Tsukuba 305–0801,
Japanffinstitutetext:  CAS Key Laboratory of Theoretical Physics, Institute of
Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China

# Novel loop-diagrammatic approach to QCD $\theta$ parameter and application
to the left-right model

Junji Hisano b,d,e,f Teppei Kitahara a Naohiro Osamura a and Atsuyuki Yamada
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

When the QCD axion is absent in full theory, the strong $CP$ problem has to be
explained by an additional mechanism, e.g., the left-right symmetry. Even
though tree-level QCD $\bar{\theta}$ parameter is restricted by the mechanism,
radiative corrections to $\bar{\theta}$ are mostly generated, which leads to a
dangerous neutron electric dipole moment (EDM). The ordinary method for
calculating the radiative $\bar{\theta}$ utilizes an equation
$\bar{\theta}=-\text{arg}\,\text{det}\,m_{q}^{\rm loop}$ based on the chiral
rotations of complex quark masses. In this paper, we point out that when full
theory includes extra heavy quarks, the ordinary method is unsettled for the
extra quark contributions and does not contain its full radiative corrections.
We formulate a novel method to calculate the radiative corrections to
$\bar{\theta}$ through a direct loop-diagrammatic approach, which should be
more robust than the ordinary one. As an application, we investigate the
radiative $\bar{\theta}$ in the minimal left-right symmetric model. We first
confirm a seminal result that two-loop level radiative $\bar{\theta}$
completely vanishes (corresponding to one-loop corrections to the quark mass
matrices). Furthermore, we estimate the size of a non-vanishing radiative
$\bar{\theta}$ at three-loop level. It is found that the resultant induced
neutron EDM is comparable to the current experimental bound, and the expected
size is restricted by the perturbative unitarity bound in the minimal left-
right model.

###### Keywords:

$CP$ violation, Electric Dipole Moments, Left-Right Models

††preprint: IPMU23-0002, KEK-TH-2494

## 1 Introduction

The QCD $\theta$ term is $P$\- and $T$-odd, and then $CP$-odd under the $CPT$
invariance. Because it is identical to the total derivative, it never locally
affects physics at the classical level (as long as the momentum conservation
holds), while its effect occurs only via nonperturbative processes [1, 2, 3,
4, 5]. It is known that this interaction induces the neutron electric dipole
moment (EDM) [6, 7]. Measurement of the neutron EDM by the nEDM collaboration
has set the severe upper bound: $|d_{n}|_{\rm exp}<1.8\times
10^{-26}\leavevmode\nobreak\ e\,\text{cm}$ (90% CL) [8]. Using the latest
lattice result $d_{n}=-0.00148(34)\,\bar{\theta}\,e\,\text{fm}$ [9]#1#1#1The
first nonzero calculation by using the lattice QCD simulation was achieved in
Ref. [10], then the first statistically significant result was obtained in
Ref. [9]. and assuming that $\bar{\theta}$ is the only source of $CP$
violation, one obtains the upper bound on the angle,

$\displaystyle|\bar{\theta}|\lesssim 1.2\times
10^{-10}\quad(90\%\leavevmode\nobreak\ \text{CL})\,,$ (1.1)

where $\bar{\theta}$ is the physical $CP$-violating angle in the QCD
Lagrangian which will be defined explicitly in the next section.

Although the experimental bound requires that $\bar{\theta}$ must be around
zero, such a $CP$-violating phase is not restricted in the Standard Model
(SM). In fact, the $CP$-violating phase in the Cabibbo-Kobayashi-Maskawa (CKM)
matrix [11, 12] is $\mathcal{O}(1)$; $\delta_{\rm CKM}\simeq 66^{\circ}=1.2$
rad [13] (in the standard parameterization [14, 15]). If there is no trick in
the full theory, $\bar{\theta}\ll 1$, or equivalently
$\bar{\theta}\ll\delta_{\rm CKM}=\mathcal{O}(1)$, requires a fine-tuning at
$\mathcal{O}(10^{-10})$ level. This is known as the strong $CP$ problem.

The massless up quark could be a solution to the strong $CP$ problem if the
observed hadron masses are explained by the nonperturbative effect [16, 17,
18, 19, 20]. However, some lattice studies ruled out this solution [21]. Thus,
the strong $CP$ problem would suggest that the SM has to be extended to
suppress $\bar{\theta}$ without the fine-tuning. Axion is the simplest
solution to the strong $CP$ problem [22, 23, 24], though it suffers from
another fine-tuning in the quantum gravity sector (axion quality problem) [25,
26, 27].

Alternatively, one may resolve the strong $CP$ problem by extended parity
symmetry [28, 29, 30, 31].#2#2#2Other possibilities are spontaneous $CP$
violation referred to as the Nelson-Barr mechanism [32, 33, 34]. In such
scenarios, the extended parity involves the left-right (LR) gauge symmetry.
The parity symmetry forbids the bare $\bar{\theta}$ parameter, while
$\mathcal{O}(1)$ of $\delta_{\rm CKM}$ is allowed. It is known that even
though the bare $\bar{\theta}$ parameter is strictly forbidden by the parity
symmetry, radiative correction to $\bar{\theta}$ is regenerated since the
parity symmetry must be softly broken in nature. Eventually, one has to
consider the experimental bound on the model from the neutron EDM measurements
in Eq. (1.1) through the radiatively regenerated $\bar{\theta}$.

The ordinary method for calculating the radiatively generated $\bar{\theta}$
parameter, adopted in many papers, utilizes

$\displaystyle\bar{\theta}=-\text{arg}\,\text{det}\,m_{u}^{\rm
loop}-\text{arg}\,\text{det}\,m_{d}^{\rm loop}\,.$ (1.2)

Here, $m_{u,d}^{\rm loop}$ are the up- and down-type quark mass matrices
including the radiative corrections. This relation is based on the chiral
rotations for the complex quark masses and an anomalous divergence of the
axial-vector current, known as the Adler-Bell-Jackiw anomaly [35, 36]. Or, it
is also derived using the path-integral formalism referred to as the Fujikawa
method [37].

The ordinary method is simple though it may not be accurate. For instance,
within the SM, even if one sets the bare $\bar{\theta}$ parameter to be zero,
it is radiatively produced via the $CP$-violating phase in the CKM matrix. It
is shown with the above method in Ref. [38] that the contribution via the
radiative corrections to quark masses is of
$\mathcal{O}(G_{F}^{2}\alpha_{s}^{3})$ ($G_{F}$ is the Fermi coupling constant
and $\alpha_{s}$ is the QCD coupling constant). On the other hand, the direct
loop calculation of the correction to $\bar{\theta}$ shows that it is derived
at four-loop order ($\mathcal{O}(G_{F}^{2}\alpha_{s})$) [39].#3#3#3A hadronic
long-distance evaluation has confirmed that the radiative $\bar{\theta}$
parameter is produced at order of $\mathcal{O}(G_{F}^{2}\alpha_{s})$ [40]. It
is consistent with the fact that the EDMs (and also the chromo-EDMs) of quarks
are induced at three-loop order ($\mathcal{O}(G_{F}^{2}\alpha_{s})$) [41]. In
fact, the ordinary method corresponds to the diagrams where the external
gluons are attached in the same fermion line in loop diagrams contributing to
$\bar{\theta}$. It implies that the leading $\bar{\theta}$ in the SM [39]
comes from diagrams with external gluons attached to different fermion lines.

In this paper, we formulate a novel approach to evaluate the radiative
corrections to $\bar{\theta}$ through a direct loop-diagrammatic calculation,
which should be more robust than the ordinary one. In Ref. [39], the external
gluon field is introduced to calculate the correction to $\bar{\theta}$ from
the CKM matrix, while details of the technique are not written. We introduce
the Fock-Swinger gauge method to directly calculate the radiative corrections
to $\bar{\theta}$ under the gluon field-strength background.

As an application, we investigate the radiative $\bar{\theta}$ in the minimal
LR symmetric model [30, 31], in which the bare and one-loop level
$\bar{\theta}$ parameters are strictly forbidden by the LR symmetry. Although
the extra heavy quarks whose Yukawa interactions violate $CP$ symmetry are
introduced, the $CP$-violating Yukawa interactions do not contribute to the
$\bar{\theta}$ parameters at one-loop level. Furthermore, it is known that
two-loop level $\bar{\theta}$ parameter also vanishes, which corresponds to
one-loop corrections to the quark masses in the ordinary method [30, 42, 43].
However, the ordinary method is unsettled for the extra quark contributions to
$\bar{\theta}$ and indeed does not contain its full radiative corrections,
like the SM calculations [39]. We first confirm this seminal result by using
the proposed method. While new type diagrams contribute to $\bar{\theta}$ at
two-loop level, the sum of the diagrams still gives no contribution to
$\bar{\theta}$. Next, we estimate the size of a non-vanishing radiative
$\bar{\theta}$ at three-loop level. It will be found that the resultant
induced neutron EDM is comparable to the current experimental bound. We will
also investigate a relation between the radiative three-loop level
$\bar{\theta}$ and the perturbative unitarity bounds of the Yukawa couplings
in the minimal LR symmetric model.

This paper is organized as follows. In Sec. 2, we discuss methods of direct
calculation of the loop diagrams contributing to $\bar{\theta}$. We show that
the operator Schwinger method and the Fock-Schwinger gauge method are
applicable, though the latter method has a merit to extend the calculation to
the higher-loop diagrams. In Sec. 3, the minimal LR symmetric model is briefly
summarized. We also derive the parameterization by the physical parameters
based on the seesaw mechanism in the LR symmetric model. We confirm that the
two-loop level radiative $\bar{\theta}$ vanishes by using the proposed method
in Sec. 4. In Sec. 5, we investigate the numerical size of the non-vanishing
radiative $\bar{\theta}$ and compare both experimental (from neutron EDM) and
theoretical bounds (from the perturbative unitarity bound). Section 6 is
devoted to conclusions and discussion. Details of the loop calculations are
given in the Appendix.

## 2 Loop-diagrammatic evaluation of QCD $\theta$ parameter

In the QCD Lagrangian, imaginary parts of the quark masses and the QCD
$\theta$ term are $P$-odd and $T$-odd interactions that are not restricted
from the $SU(3)_{C}$ gauge symmetry,

$\displaystyle\mathcal{L}_{{\ooalign{\hfil/\hfil\crcr$P$}},\,{\ooalign{\hfil/\hfil\crcr$T$}}}=-\sum_{q=\text{all}}\text{Im}(m_{q})\bar{q}i\gamma_{5}q+\theta_{G}\frac{\alpha_{s}}{8\pi}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}\,,$
(2.5)

where $m_{q}$ stands for the complex quark masses with
$m_{q}\equiv|m_{q}|\exp(i\theta_{q})$, $G_{\mu\nu}^{a}$ is the gluon field-
strength tensor,
$\tilde{G}^{a\mu\nu}\equiv\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}G^{a}_{\rho\sigma}$
with $\epsilon^{0123}=+1$, and $\alpha_{s}=g_{s}^{2}/(4\pi)$ is the
$SU(3)_{C}$ coupling constant. It is well-known that the axial rotation of
quarks

$\displaystyle q\to
q^{\prime}=\exp{\left(-\frac{i}{2}\theta_{q}\gamma_{5}\right)}\,q\,,\qquad\bar{q}\to\bar{q}^{\prime}=\bar{q}\exp{\left(-\frac{i}{2}\theta_{q}\gamma_{5}\right)}\,,$
(2.6)

turns off the imaginary part of the quark masses and generates an additional
QCD $\theta$ term,

$\displaystyle\mathcal{L}_{{\ooalign{\hfil/\hfil\crcr$P$}},\,{\ooalign{\hfil/\hfil\crcr$T$}}}={\bar{\theta}}\,\frac{\alpha_{s}}{8\pi}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}\,,$
(2.11)

where

$\displaystyle\bar{\theta}\equiv\theta_{G}-\sum_{q}\theta_{q}\,,$ (2.12)

is a physical $\bar{\theta}$ parameter, if all quarks are massive. This is
derived using the path-integral formalism referred to as the Fujikawa method
[37] or with the Adler-Bell-Jackiw anomaly [35, 36, 44, 45]

$\displaystyle\partial_{\mu}(\bar{q}\gamma^{\mu}\gamma_{5}q)=2i\mathop{\mathrm{Re}}(m_{q})\bar{q}\gamma_{5}q-2\mathop{\mathrm{Im}}(m_{q})\bar{q}q-\frac{\alpha_{s}}{4\pi}G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}\,.$
(2.13)

The contribution of the quark mass phase to the QCD $\theta$ term should be
able to be directly evaluated by loop-diagrammatically integrating out quarks,
not via the Adler-Bell-Jackiw anomaly or transformation of measure in the path
integral (the Fujikawa method). However, it does not generate the QCD $\theta$
term if the momenta in the diagrams are conserved. It is because the $\theta$
term is equivalent to total-derivative and the total momentum has to be zero.
Thus, we have to abandon the momentum conservation or equivalently the
translation invariance in order to evaluate the QCD $\theta$ term with the
loop-diagrammatic calculation.

It can be realized by introducing the gluon field strength background. In this
section, we evaluate the QCD $\theta$ term with two different methods, 1) the
operator Schwinger method and 2) the Fock-Schwinger gauge method. We show that
they produce consistent results with Eq. (2.12).#4#4#4An alternative way to
evaluate the QCD $\theta$ term is the $CP$-odd spurion trick [46]. Equation
(2.12) can be derived by introducing a spurion, whose vacuum expectation value
produces the $CP$-violating phase of the quark mass, with an external momentum
injection via the spurion.

### 2.1 Operator Schwinger method

First, we consider the operator Schwinger method.#5#5#5See Ref. [47] for the
review about the operator Schwinger method. The effective action $\Delta S$
induced by the integration of quarks at one-loop level is given by the log-
determinant (or trace-log) of the Dirac operator. Now we introduce the complex
mass parameters for quarks. In this case, the effective action is given as

$\displaystyle\Delta S$ $\displaystyle=$ $\displaystyle-i{\rm
Tr}\,\log\left|\left|{\ooalign{\hfil/\hfil\crcr$P$}}-(m_{q}^{\ast}P_{L}+m_{q}P_{R})\right|\right|\,,$
(2.16)

where $P_{\mu}=iD_{\mu}$, $D_{\mu}\equiv\partial_{\mu}+ig_{s}T^{a}G_{\mu}^{a}$
is the QCD covariant derivative, and $P_{L/R}=(1\mp\gamma_{5})/2$. In this
method, the following basic commutation relation is used,

$\displaystyle[{P}_{\mu},{P}_{\nu}]$ $\displaystyle=$ $\displaystyle-
ig_{s}T^{a}G^{a}_{\mu\nu}(\equiv-ig_{s}G_{\mu\nu})\,,$ (2.17)

since the gluon field-strength tensor appears from it. Then, the derivative of
$\Delta S$ over the fermion mass $m_{q}$ for $P_{R}$ is

$\displaystyle\frac{d}{dm_{q}}\Delta S$ $\displaystyle=$ $\displaystyle i{\rm
Tr}\left[\frac{1}{{\ooalign{\hfil/\hfil\crcr${P}$}}-(m_{q}^{\ast}P_{L}+m_{q}P_{R})}P_{R}\right]$
(2.20) $\displaystyle=$ $\displaystyle i{\rm
Tr}\left[\frac{1}{{P}^{2}-|m_{q}|^{2}-\frac{1}{2}g_{s}\sigma_{\mu\nu}G^{\mu\nu}}m_{q}^{\ast}P_{R}\right]\,,$
(2.21)

where $\sigma_{\mu\nu}=\frac{i}{2}[\gamma_{\mu},\gamma_{\nu}]$. Since the
Levi-Civita tensor appears from the trace of a product of four $\gamma$’s and
$\gamma_{5}$, the second order of $\frac{1}{2}g_{s}\sigma_{\mu\nu}G^{\mu\nu}$
leads to the $G\tilde{G}$ term as

$\displaystyle\frac{d}{dm_{q}}\Delta S\supset$ $\displaystyle i{\rm
Tr}\left[\frac{1}{{P}^{2}-|m_{q}|^{2}}\left(\frac{1}{2}g_{s}\sigma_{\mu\nu}G^{\mu\nu}\right)\frac{1}{{P}^{2}-|m_{q}|^{2}}\left(\frac{1}{2}g_{s}\sigma_{\rho\sigma}G^{\rho\sigma}\right)\frac{1}{{P}^{2}-|m_{q}|^{2}}m_{q}^{\ast}P_{R}\right]$
$\displaystyle=$ $\displaystyle\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\frac{g_{s}^{2}}{2}\frac{-m_{q}^{\ast}}{(p^{2}-|m_{q}|^{2})^{3}}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}+\cdots$
$\displaystyle\supset$ $\displaystyle\int
d^{4}x\frac{i}{32\pi^{2}}\frac{g_{s}^{2}}{2}\frac{1}{m_{q}}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}\,,$
(2.22)

where $\text{Tr}(T^{a}T^{b})=(1/2)\delta^{ab}$ is used. Here, $P^{2}$ is
replaced by $-\partial^{2}$, and it is integrated in momentum space.
Similarly, $d\Delta S/dm_{q}^{\ast}$ leads to the ${G}\tilde{G}$ term. By
Integrating $d\Delta S/dm_{q}$ with $m_{q}$ and $d\Delta S/dm_{q}^{\ast}$ with
$m_{q}^{\ast}$, we get

$\displaystyle\Delta\mathcal{L}=$
$\displaystyle\frac{-i}{32\pi^{2}}\frac{g_{s}^{2}}{2}\log\frac{m_{q}^{\ast}}{m_{q}}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}$
$\displaystyle=$
$\displaystyle-{\theta_{q}}\frac{\alpha_{s}}{8\pi}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}\,.$
(2.23)

Now we obtain the contribution from a quark with complex mass to the QCD
$\theta$ term by integrating out the quark, which is consistent with the axial
rotation (2.11).

### 2.2 Fock-Schwinger gauge method

In the previous section, we showed that the operator Schwinger method enables
us to derive the physical QCD $\theta$ term by the diagrammatic evaluation.
However, the operator Schwinger method is not suitable to calculate effective
operators induced at higher loops because it is the method to obtain an
effective action by integrating fermions out with the log-determinant of the
Dirac operator. Here alternatively, we introduce the Fock-Schwinger gauge
method, which is more applicable to diagrammatic calculation.#6#6#6See Ref.
[47] for the review about the Fock-Schwinger gauge method.

The Fock-Schwinger gauge is to take such a gauge

$\displaystyle(x^{\mu}-x^{\mu}_{0})G_{\mu}^{a}(x)=0\,,$ (2.24)

which violates the translation symmetry. Because of breaking the translation
symmetry, we can derive perturbatively the QCD $\theta$ term as will be shown
below. While the Fock-Schwinger gauge violates the translation symmetry, the
physical observables do not depend on it. In the below argument, we take gauge
dependence parameter $x_{0}$ as $x_{0}=0$ for simplicity.

The gluon field $G_{\mu}^{a}$ can be expanded under this gauge around $x=0$
and it is given with the gluon field-strength tensor at $x=0$,
$G^{a}_{\mu\nu}(0)$, as [47]

$\displaystyle G_{\mu}^{a}(x)$
$\displaystyle=\frac{1}{2}x^{\nu}G_{\nu\mu}^{a}(0)+\cdots$ $\displaystyle=\int
d^{4}ke^{-ik\cdot
x}\left(-\frac{i}{2}G^{a}_{\nu\mu}(0)\frac{\partial}{\partial
k_{\nu}}\delta^{(4)}(k)\right)+\cdots\,.$ (2.25)

Here, the discarded terms are covariant derivatives of the background gluon
field-strength tensor, which are irrelevant to the calculation of the QCD
$\theta$ term. We can systematically evaluate the interaction of the
propagating quarks with the background gluon field-strength tensor in this
gauge fixing. However, we found that the effective gluon operators such as the
QCD $\theta$ term cannot be evaluated from the simple quark bubble diagrams.
The background gluon fields bring momenta, $k$ in Eq. (2.25), which are taken
to be zero in the last step of the calculation due to $\delta^{(4)}(k)$. Thus,
the quark momentum is not constant due to interaction with the background
field and the quark line cannot be closed without violating momentum
conservation.

Figure 1: Feynman amplitude $i\Pi_{X}^{q}$ for the loop-diagrammatic
evaluation of the $\bar{\theta}$ parameter.

In order to fix this problem, we introduce an auxiliary (dimensionless)
background field $X$, and it is coupled to the $CP$-odd quark mass terms as

$\displaystyle\mathcal{L}_{{\ooalign{\hfil/\hfil\crcr$P$}},\,{\ooalign{\hfil/\hfil\crcr$T$}}}$
$\displaystyle=-\sum_{q=\text{all}}\text{Im}(m_{q})\bar{q}i\gamma_{5}q$ (2.30)
$\displaystyle\Rightarrow-\sum_{q=\text{all}}\text{Im}(m_{q})\left(\bar{q}i\gamma_{5}q\right)X\,.$
(2.31)

We evaluate the leading contribution of $\text{Im}(m_{q})$ to the QCD $\theta$
term in perturbative way, assuming $\text{Im}(m_{q})\ll\text{Re}(m_{q})$. The
field $X$ is taken to be $1$ in the last step of the calculation.#7#7#7This
technique has been applied for evaluation of the Weinberg operator in the QCD
[48].

The radiative QCD $\theta$ term comes from a bubble diagram. The Feynman
diagram in Fig. 1 shows the leading contribution, which is realized by
integrating the delta function in Eq. (2.25) as

$\begin{split}&i\Pi_{X}^{q}=-\int
d^{4}k_{1}d^{4}k_{2}\frac{d^{4}p}{(2\pi)^{4}}\tilde{X}(-k_{1}-k_{2})\\\
&\quad\times{\rm
Tr}\left[(-ig_{s}\gamma^{\mu}T^{a})\left(-\frac{i}{2}G^{a}_{\rho\mu}(0)\frac{\partial}{\partial
k_{1;\rho}}\delta^{(4)}(k_{1})\right)\frac{i}{\not{p}+\not{k}_{1}-{\rm
Re}(m_{q})}{\rm Im}(m_{q})\gamma_{5}\right.\\\
&\quad\times\left.\frac{i}{\not{p}-\not{k}_{2}-{\rm
Re}(m_{q})}(-ig_{s}\gamma^{\nu}T^{b})\left(-\frac{i}{2}G^{b}_{\sigma\nu}(0)\frac{\partial}{\partial
k_{2;\sigma}}\delta^{(4)}(k_{2})\right)\frac{i}{\not{p}-{\rm
Re}(m_{q})}\right]\,,\end{split}$ (2.32)

where $\tilde{X}(k)=\int d^{4}xe^{ik\cdot x}X(x)$ and $p$ is the loop
momentum. We followed the Feynman rules under the Fock-Schwinger gauge, which
includes the gluon field $G^{a}_{\mu}$ expressed as Eq. (2.25) and the
modified $CP$-odd quark mass term in Eq. (2.31). Until integration of the
delta functions, two independent momenta $k_{1},k_{2}$ flow into the vertices
with the background field-strength tensors $G_{\rho\mu}^{a}$ and
$G_{\sigma\nu}^{b}$, respectively, and the artificial background field $X$
brings a momentum $-k_{1}-k_{2}$ (see Fig. 1). After some calculation, we get

$\begin{split}i\Pi_{X}^{q}&=i\,{\rm
Im}(m_{q})\frac{g_{s}^{2}}{8}G^{a}_{\rho\mu}(0)G^{a}_{\sigma\nu}(0)\int\frac{d^{4}p}{(2\pi)^{4}}\frac{4i{\rm
Re}(m_{q})\epsilon^{\mu\nu\rho\sigma}}{[p^{2}-({\rm
Re}(m_{q}))^{2}]^{3}}\tilde{X}(0)\\\
&=i\frac{\alpha_{s}}{8\pi}\left(-\frac{{\rm Im}(m_{q})}{{\rm
Re}(m_{q})}\right)G^{a}_{\mu\nu}(0)\tilde{G}^{a\mu\nu}(0)\\\
&\simeq-i\theta_{q}\frac{\alpha_{s}}{8\pi}G^{a}_{\mu\nu}(0)\tilde{G}^{a\mu\nu}(0)\,.\end{split}$
(2.33)

Here, we take $\tilde{X}(0)=1$.

After integrating the quark $q$ out in the full theory,
$\Delta\mathcal{L}=\Pi_{X}^{q}$ is obtained in the effective action of the
gluon. Eventually, one can derive the QCD $\theta$ term in the Fock-Schwinger
gauge method. This result is consistent with that of the chiral rotation,
$q\to q^{\prime}=\exp(-\frac{i}{2}\theta_{q}\gamma_{5})\,q$, in Eq. (2.11) and
also the operator Schwinger method in Eq. (2.23). Hence, we reached a
clarification of the equivalence among the Fock-Schwinger gauge method, the
operator Schwinger method, and the ordinary chiral rotation, and we noticed
that the Fock-Schwinger gauge method is more intuitive than the operator
Schwinger method for higher-loop order calculations.

It might be concerned that the diagrammatic evaluations of the light quark
contribution to the QCD $\theta$ term is not justified from the viewpoint of
perturbation, since the loop momentum around the quark mass dominates the
integrals in Eqs. (2.22) and (2.33). It might be healthy to evaluate the light
quark mass phases above the $\Lambda_{\rm QCD}$ scale and derive the QCD
$\theta$ parameter by the chiral rotation. However, since the diagrammatic
evaluations are consistent with those of the chiral rotation, we may forget
the problem in practical cases.

In this paper, to evaluate the QCD $\theta$ term diagrammatically, we will use
the Fock-Schwinger gauge method. Note that the auxiliary background field $X$
should be attached to any perturbative interactions, but we suppress them in
the following calculations for the sake of clarity.

## 3 The minimal left-right symmetric model

### 3.1 Model

From this section, we introduce the minimal LR symmetric model that can solve
the strong $CP$ problem. The LR symmetry, which is formed by introducing a new
$SU(2)_{R}$ gauge symmetry, with spatial parity symmetry is motivated to
forbid the QCD $\bar{\theta}$ term at tree level. In particular, we focus on
the minimal LR symmetric model, which embeds the $SU(2)_{L}$ singlet right-
handed quarks, $u_{R}$ and $d_{R}$, to the $SU(2)_{R}$ doublets,
$Q_{R}\equiv(u_{R},d_{R})^{T}$. Furthermore, a $SU(2)_{R}$ doublet Higgs,
$H^{\prime}$, and three flavors of the up-type and down-type vector-like
quarks, $U_{L},U_{R},D_{L}$ and $D_{R}$, have to be introduced. The matter
contents are listed in Table 1.

---
| $SU(3)_{C}$ | $SU(2)_{L}$ | $SU(2)_{R}$ | $U(1)_{B-L}$ | $U(1)_{Y}$
$Q_{L}^{i}\equiv(u_{L}^{i},d_{L}^{i})^{T}$ | $\bf 3$ | $\bf 2$ | $\bf 1$ | $1/6$ | $(1/6,1/6)$
$Q_{R}^{i}\equiv(u_{R}^{i},d_{R}^{i})^{T}$ | $\bf 3$ | $\bf 1$ | $\bf 2$ | $1/6$ | $(2/3,-1/3)$
$H$ | $\bf 1$ | $\bf 2$ | $\bf 1$ | $1/2$ | $(1/2,1/2)$
$H^{\prime}$ | $\bf 1$ | $\bf 1$ | $\bf 2$ | $1/2$ | $(1,0)$
$U_{L}^{a}$ | $\bf 3$ | $\bf 1$ | $\bf 1$ | $2/3$ | $2/3$
$U_{R}^{a}$ | $\bf 3$ | $\bf 1$ | $\bf 1$ | $2/3$ | $2/3$
$D_{L}^{a}$ | $\bf 3$ | $\bf 1$ | $\bf 1$ | $-1/3$ | $-1/3$
$D_{R}^{a}$ | $\bf 3$ | $\bf 1$ | $\bf 1$ | $-1/3$ | $-1/3$
Table 1: The matter contents and their gauge charges in the minimal LR
symmetric model, where $U(1)_{Y}=T^{R}_{3}\,+\,U(1)_{B-L}$. The indices $i$
and $a$ represent the flavors for the doublet and singlet quarks,
respectively.

To solve the strong $CP$ problem, the spatial parity symmetry has to be
extended to symmetrize the left-handed and right-handed sectors as well as the
$SU(2)_{L}$ and $SU(2)_{R}$ gauge bosons,

$\displaystyle\vec{x}$ $\displaystyle\longleftrightarrow$
$\displaystyle-\vec{x}\,,$ $\displaystyle W_{\mu}$
$\displaystyle\longleftrightarrow$ $\displaystyle W^{\prime\mu}\,,$ (3.1)
$\displaystyle Q_{L},\leavevmode\nobreak\ U_{L},\leavevmode\nobreak\
D_{L},\leavevmode\nobreak\ H$ $\displaystyle\longleftrightarrow$
$\displaystyle Q_{R},\leavevmode\nobreak\ U_{R},\leavevmode\nobreak\
D_{R},\leavevmode\nobreak\ H^{\prime}\,,$

while the other gauge bosons are invariant. The spontaneous violation of the
extended parity symmetry, $SU(2)_{R}\times U(1)_{B-L}\to U(1)_{Y}$, is caused
by the vacuum expectation value (VEV) of $H^{\prime\,0}$, $\langle
H^{\prime}\rangle=(0,v^{\prime})$. After the symmetry breaking, the $U(1)_{Y}$
gauge symmetry is generated with the gauge charge of
$U(1)_{Y}=T^{R}_{3}\,+\,U(1)_{B-L}$. The $SU(2)_{R}$ gauge bosons absorb the
Nambu-Goldstone (NG) bosons in the doublet $H^{\prime}$ ($\varphi^{\prime\,+}$
and $\varphi^{\prime\,0}$) to become massive states ($W^{\prime\,+}$ and
$Z^{\prime}$). The physical neutral Higgs boson associated with this symmetry
breaking is denoted as $h^{\prime}$. Then, the $SU(2)_{L}\times U(1)_{Y}$
gauge symmetry is broken to $U(1)_{\rm EM}$ by the VEV of $H^{0}$, $\langle
H\rangle=(0,v)$. The $W^{+}$ and $Z$ bosons absorb the NG bosons in the
doublet $H$ ($\varphi^{+}$ and $\varphi^{0}$), and the physical (SM) neutral
Higgs boson with this symmetry breaking is denoted as $h$.#8#8#8 The
$h$–$h^{\prime}$ and $Z$–$Z^{\prime}$ mixings are induced in the model at tree
level, though the mixings are suppressed by $v/v^{\prime}$ and
$(v/v^{\prime})^{2}$, respectively [30]. Since we take
$v/v^{\prime}\rightarrow 0$ in the calculation of the $\bar{\theta}$
parameter, they are ignored.

The resultant parity violation in nature comes from $v^{\prime}\neq v$.
Namely, we assume that soft parity breaking terms are contained in the $H$ and
$H^{\prime}$ Higgs potentials which lead to $v^{\prime}\gg v\neq 0$. These two
VEVs can be chosen as real and positive without loss of generality. Since the
extended parity is a discrete symmetry, its spontaneous breaking leads to the
formation of the domain walls, which dominate the energy density in the
Universe. This domain wall problem can be naturally solved by the Planck
suppressed higher-dimensional operators, which explicitly violate the parity
symmetry [43].

The Yukawa interactions and Dirac mass terms for the vector-like quarks are
represented as

$\begin{split}-{\cal
L}_{Y}=&\overline{Q}_{L}^{i}x^{ia}_{u}U^{a}_{R}\tilde{H}+\overline{Q}_{R}^{i}x^{ia}_{u}U^{a}_{L}\tilde{H}^{\prime}+M^{a}_{u}\overline{U}_{L}^{a}U_{R}^{a}\\\
&+\overline{Q}_{L}^{i}x^{ia}_{d}D^{a}_{R}{H}+\overline{Q}_{R}^{i}x^{ia}_{d}D^{a}_{L}{H}^{\prime}+M^{a}_{d}\overline{D}_{L}^{a}D_{R}^{a}+{\rm
h.c.}\,,\end{split}$ (3.2)

where $i=1$–$3$ is a flavor index for $SU(2)_{L/R}$ doublets, $a=1$–$3$ is
that for the singlets (vector-like quarks), and $\tilde{H}^{(\prime)}=\epsilon
H^{(\prime)\ast}$ ($\epsilon_{12}=1$). The LR symmetry requires that the
Yukawa $x^{ia}_{u/d}$ in the first two terms (in both lines) must be the same
complex matrices, and the Dirac mass terms $M_{u}$ and $M_{d}$ must be
Hermitian. The Dirac mass terms $M_{u}^{a}$ and $M_{d}^{a}$ in Eq. (3.2) are
diagonalized to real and positive eigenvalues by the field redefinitions of
the vector-like quarks.#9#9#9One can also consider non-Hermitian vector-like
quark mass matrices which correspond to soft parity breaking terms [30].
However, such contributions produce large quark EDM and radiative
$\bar{\theta}$, and thus they are severely constrained from the EDM bounds
[43, 49].  The SM quark masses are realized by the seesaw mechanism such as
higher dimensional operators induced by integrating out the vector-like
quarks.#10#10#10One can also extend the lepton sector that is insensitive to
the QCD $\bar{\theta}$ term. If one considers $SU(5)_{L}\times SU(5)_{R}$
grand unification [50, 51, 52], vector-like neutral leptons are absent and the
neutrinos keep massless at the tree level. Interestingly, suitable Dirac
neutrino masses are generated from the two-loop radiative corrections [53]
with predicting a nonzero $\Delta N_{\rm eff}$ [54]. Furthermore, an
$\mathcal{O}(10)$ keV sterile neutrino dark matter with the leptogenesis
mechanism can be incorporated [55, 56].

Before discussing the mass matrices in detail in the next section, let us
count on the number of physical $CP$ phases in this model. The Yukawa
couplings $x_{u/d}^{ia}$ are $3\times 3$ complex matrices. Nine real
parameters are removed from the Yukawa matrices by field redefinition of
$Q_{L/R}^{i}$ as $Q_{L/R}^{i}\rightarrow U^{ij}Q_{L/R}^{j}$ with a unitary
matrix $U$. Furthermore, phase redefinition of $U_{L/R}^{a}$ and $D_{L/R}^{a}$
removes five phases in the total. A remaining phase rotation corresponds to
the baryon number conservation, and it does not change $x_{u}$ and $x_{d}$.
Thus, $x_{u}$ and $x_{d}$ have a total of 22 physical real parameters. We
parametrize these 22 parameters as

$\displaystyle\begin{aligned}
x_{u}&=\overline{\Phi}^{\dagger}(\theta_{d3},\theta_{d8})\,V_{Q}\,\overline{\Phi}(\theta_{u3},\theta_{u8})\,\bar{x}_{u}V_{U}\,,\\\
x_{d}&=\bar{x}_{d}V_{D}\,,\end{aligned}$ (3.3)

with

$\displaystyle\overline{\Phi}(\theta_{3},\theta_{8})$ $\displaystyle\equiv$
$\displaystyle\exp(i\tau_{3}\theta_{3})\exp(i\tau_{8}\theta_{8})\,,$ (3.4)

where $\tau_{3}$ and $\tau_{8}$ are the third and eighth Gell-Mann matrices.
Here, $\bar{x}_{u/d}$ are real diagonal matrices and $V_{Q}$, $V_{U}$, and
$V_{D}$ are CKM-like unitary matrices which have three rotation angles and one
$CP$-violating phase. It is found that there are seven $CP$-violating phases
in this model ($\theta_{u3/u8}$, $\theta_{d3/d8}$, and three phases in
$V_{Q/U/D}$). When the Dirac masses are assumed to be universal such as
$M_{u}^{a}=M_{u}$ and $M_{d}^{a}=M_{d}$ for $a=1$–$3$, the parameters
$\theta_{u3/u8}$, $\theta_{d3/d8}$, $V_{U/D}$ become unphysical and only
$V_{Q}$ remains physical.

In this paper, we assume that
$v^{\prime}\mathop{}_{\textstyle\sim}^{\textstyle<}M_{q}^{a}$ and will utilize
expansions by $v^{\prime}/M_{q}^{a}$. This inequality is motivated because the
seesaw mechanism may explain the SM fermion mass hierarchy naturally. On the
other hand, if $v^{\prime}\gg M_{q}^{a}$, a copy of the SM fermions has a mass
spectrum similar to the SM fermions, which spread over five orders of
magnitude. A new naturalness problem might appear in such a model, but the QCD
$\bar{\theta}$ term is suppressed by $M_{q}^{a}/v^{\prime}$ since only the CKM
phase survives in a limit of $M_{q}^{a}\rightarrow 0$. In the following, we
will consider the case of
$v^{\prime}\mathop{}_{\textstyle\sim}^{\textstyle<}M_{q}^{a}$.

### 3.2 Parametrization of Yukawa coupling constants

In this section, we show the quark mass matrices and define the mass
eigenstates. In the mass matrices the Yukawa coupling constants, $x_{u}$ and
$x_{d}$, appear, and we have to determine them from the observed quark masses
and CKM matrix in order to evaluate the radiative $\bar{\theta}$ parameter. We
give the parameterization of the Yukawa coupling constants assuming the SM
quark masses are given by the seesaw mechanism with $v^{\prime}\lesssim
M_{q}^{a}$.

From Eq. (3.2), the quark mass matrices in the flavor eigenstates are given as

$\displaystyle-{\cal L}_{M}$
$\displaystyle=\left(\overline{u}_{L}^{i},\overline{U}_{L}^{a}\right)\begin{pmatrix}0&x^{ib}_{u}v\\\
x^{\dagger
aj}_{u}v^{\prime}&M^{a}_{u}\delta^{ab}\end{pmatrix}\begin{pmatrix}u_{R}^{j}\\\
U_{R}^{b}\end{pmatrix}+\left(\overline{d}_{L}^{i},\overline{D}_{L}^{a}\right)\begin{pmatrix}0&x^{ib}_{d}v\\\
x^{\dagger
aj}_{d}v^{\prime}&M^{a}_{d}\delta^{ab}\end{pmatrix}\begin{pmatrix}d_{R}^{j}\\\
D_{R}^{b}\end{pmatrix}+{\rm h.c.}$
$\displaystyle\equiv\overline{\mathcal{U}}_{L}^{p}\mathcal{M}_{u}^{(0)pq}\mathcal{U}_{R}^{q}+\overline{\mathcal{D}}_{L}^{p}\mathcal{M}_{d}^{(0)pq}\mathcal{D}_{R}^{q}+{\rm
h.c.}\,,$ (3.5)

where ${\mathcal{U}}_{L/R}^{p}$ and ${\mathcal{D}}_{L/R}^{p}$
($p,q=1,\cdots,6$) are the up- and down-type flavor eigenstates, and $v\simeq
174.1\,$GeV. Here, $M^{a}_{u/d}$ are real diagonal, while $x^{ia}_{u/d}$ are
complex matrices. It is obvious that $\text{arg}\,\text{det}\,{\cal
M}^{(0)}_{u/d}=0$. Then, the $6\times 6$ fermion mass matrices are
diagonalized by bi-unitary matrices, $V_{qL}$ and $V_{qR}$, as

$\displaystyle{\cal M}_{q}^{(0)pq}$ $\displaystyle=$ $\displaystyle
V_{qL}^{\dagger
pP}\bar{M}_{q}^{P}V_{qR}^{Pq}\quad\text{for}\leavevmode\nobreak\
q=u\text{\leavevmode\nobreak\ and\leavevmode\nobreak\ }d\,,$ (3.6)

with diagonal mass matrices $\bar{M}_{q}$. The mass eigenstates,
${\mathcal{U}}_{ML/R}^{P}$ and ${\mathcal{D}}_{ML/R}^{P}$ for $P=1$–$6$, are
given as

It is difficult to reconstruct the model parameters from the experimental data
in general. Here we assume that
$v^{\prime}\mathop{}_{\textstyle\sim}^{\textstyle<}M_{q}^{a}$ and we take
leading terms in the expansion of $v^{\prime}/M^{a}_{q}$ for the quark mass
eigenvalues. In this expansion, the SM quark masses are given by the following
seesaw relation

$V_{q}^{\dagger
iI}\frac{m_{q}^{I}}{v}V_{q}^{Ij}=x_{q}^{ia}\frac{v^{\prime}}{M_{q}^{a}}x_{q}^{\dagger
aj}\,,$ (3.7)

with a $3\times 3$ unitary matrix $V_{q}$ and $I=1$–$3$, while the heavy quark
masses are given by $M_{q}^{a}$, for $q=u$ and $d$.

Now let us rewrite Eq. (3.7) as

$\displaystyle{\bm{I}}_{3}$ $\displaystyle=$
$\displaystyle\left(\frac{\sqrt{v}}{\sqrt{m_{q}}}V_{q}x_{q}\frac{\sqrt{v^{\prime}}}{\sqrt{M_{q}}}\right)\left(\frac{\sqrt{v^{\prime}}}{\sqrt{M_{q}}}x^{\dagger}_{q}V^{\dagger}_{q}\frac{\sqrt{v}}{\sqrt{m_{q}}}\right)\equiv
U_{q}U^{\dagger}_{q}\,,$ (3.8)

where ${\bm{I}}_{3}$ is a unit matrix in the three-dimensional space and we
ignore the indices of matrices.#11#11#11 A similar parameterization technique,
referred to as the Casas-Ibarra parameterization, is applied in the minimal
seesaw model [57, 58].  Thus, the Yukawa matrices $x_{q}$ are given with a
unitary matrix $U_{q}$ by

$\displaystyle
x_{q}=V^{\dagger}_{q}\frac{\sqrt{m_{q}}}{\sqrt{v}}U_{q}\frac{\sqrt{M_{q}}}{\sqrt{v^{\prime}}}\,.$
(3.9)

According to the previous section, one can remove some unphysical parameters
in $U_{q}$ and $V_{q}$ by the field redefinitions. Then, we get

$\displaystyle\begin{aligned} x_{u}&=V_{\rm
CKM}^{\dagger}\frac{\sqrt{m_{u}}}{\sqrt{v}}\overline{\Phi}(\theta_{u3},\theta_{u8})V_{U}\frac{\sqrt{M_{u}}}{\sqrt{v^{\prime}}}\,,\\\
x_{d}&=\frac{\sqrt{m_{d}}}{\sqrt{v}}\overline{\Phi}(\theta_{d3},\theta_{d8})V_{D}\frac{\sqrt{M_{d}}}{\sqrt{v^{\prime}}}\,,\end{aligned}$
(3.10)

where $V_{\rm CKM}(\equiv V_{u}V_{d}^{\dagger})$ corresponds to the CKM matrix
in the SM. Here, $V_{U/D}$ are CKM-like unitary matrices with three mixing
angles and one $CP$-violating phase, though they are different from those in
Eq. (3.3). Now we have seven physical $CP$-violating phases ($\theta_{u3/u8}$,
$\theta_{d3/d8}$, and three phases in $V_{U/D}$ and $V_{\rm CKM}$), which is
consistent with our previous counting, and all phases can be $\mathcal{O}(1)$
under the extended parity symmetry.

Since we assume that
$v^{\prime}\mathop{}_{\textstyle\sim}^{\textstyle<}M_{q}^{a}$, the $6\times 6$
diagonalization matrices in Eq. (LABEL:eq:diagonalization) are given of
leading terms in the expansion of $v^{\prime}/M^{a}_{q}$ as

$\displaystyle\begin{aligned} V_{uL}&\simeq\begin{pmatrix}-V_{\rm CKM}&V_{\rm
CKM}x_{u}\frac{v}{M_{u}}\\\
\frac{v}{M_{u}}x_{u}^{\dagger}&{\bm{I}}_{3}\end{pmatrix}\,,\\\
V_{uR}&\simeq\begin{pmatrix}V_{\rm CKM}&-V_{\rm
CKM}x_{u}\frac{v^{\prime}}{M_{u}}\\\
\frac{v^{\prime}}{M_{u}}x_{u}^{\dagger}&{\bm{I}}_{3}\end{pmatrix}\,,\\\
V_{dL}&\simeq\begin{pmatrix}-{\bm{I}}_{3}&x_{d}\frac{v}{M_{d}}\\\
\frac{v}{M_{d}}x_{d}^{\dagger}&{\bm{I}}_{3}\end{pmatrix}\,,\\\
V_{dR}&\simeq\begin{pmatrix}{\bm{I}}_{3}&-x_{d}\frac{v^{\prime}}{M_{d}}\\\
\frac{v^{\prime}}{M_{d}}x_{d}^{\dagger}&{\bm{I}}_{3}\end{pmatrix}\,,\end{aligned}$
(3.11)

where $x_{q}$ are given by Eq. (3.10). Here, the diagonal eigenvalue matrices
$\bar{M}_{q}$ are given as,

$\displaystyle\bar{M}_{q}^{P}=\text{diag}(m_{q}^{I},M_{q}^{a})\quad\text{for}\leavevmode\nobreak\
q=u\text{\leavevmode\nobreak\ and\leavevmode\nobreak\ }d\,.$ (3.12)

### 3.3 Quark EDMs

Before discussing the radiative $\bar{\theta}$ parameter in the minimal LR
symmetric model, let us comment on contributions to quark (chromo) EDMs. As
long as the vector-like mass matrices are Hermitian, the quark EDMs vanish
completely at one-loop level. The $W^{(\prime)\pm}$ contribution at one-loop
level vanishes trivially since the chirality is conserved in the diagrams. On
the other hand, the one-loop quark EDM contributions from neutral Higgses and
$Z^{(\prime)}$ may have a chirality flip in the diagrams. Nevertheless, they
also vanish because the extended parity symmetry restricts the product of two
vertices to be strictly real, as shown in Ref. [43].

## 4 Confirmation of vanishing QCD $\theta$ parameter in two-loop order

The parity symmetry is spontaneously broken by the $\langle
H^{\prime}\rangle\leavevmode\nobreak\ (\gg\langle H\rangle)$ in the LR
symmetric model. Since $\text{arg}\,\text{det}\,{\cal M}^{(0)}_{u/d}=0$ holds,
fermion one-loop contributions to the QCD $\bar{\theta}$ term remain zero.
However, it is expected that fermion-loop diagrams at a higher than one-loop
level would generate it. In this section, we show fermion two-loop
contributions to the QCD $\bar{\theta}$ term still vanish.

Integrating out quarks, the following higher-dimensional operators are
expected to be generated,

$\displaystyle{\cal L}_{\rm eff}$ $\displaystyle=$
$\displaystyle\sum_{n=1}C_{n}\frac{|H^{\prime}|^{2n}-|H|^{2n}}{M_{q}^{2n}}\frac{\alpha_{s}}{8\pi}G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}\,,$
(4.1)

with $M_{q}$ as the scale of vector-like quark masses, and the QCD
$\bar{\theta}$ term are induced by the spontaneous parity symmetry breaking
$\langle H^{\prime}\rangle\leavevmode\nobreak\ (\gg\langle H\rangle)$,

$\displaystyle\bar{\theta}$ $\displaystyle=$
$\displaystyle\sum_{n=1}C_{n}\left(\frac{\langle
H^{\prime}\rangle}{M_{q}}\right)^{2n}.$ (4.2)

We evaluate the Wilson coefficients of the operators $C_{n}$ in the following.

First, we consider the contributions to the QCD $\bar{\theta}$ term at two-
loop level coming from an exchange of the $W^{\prime\,\pm}$ boson. The two-
loop fermion bubble diagrams mediated by the $W^{\prime\,\pm}$ boson under the
gluon field-strength background conserve chirality in the fermion line, and
then it is proportional to

$\displaystyle V_{dR}^{Pi}V_{uR}^{\dagger iQ}V_{uR}^{Qj}V_{dR}^{\dagger
jP}f\left[(\bar{M}_{d}^{P})^{2},(\bar{M}_{u}^{Q})^{2},m_{W^{\prime}}^{2}\right]\,,$
(4.3)

where $i$ and $j$ run $1$–$3$ as the flavor index for the $SU(2)_{R}$ doublet,
while $P$ and $Q$ run $1$–$6$ for the quark mass eigenstates. Here, a two-loop
function $f[(\bar{M}_{d}^{P})^{2},(\bar{M}_{u}^{Q})^{2},m_{W^{\prime}}^{2}]$
is a real function. Since $\left(V_{dR}^{Pi}V_{uR}^{\dagger
iQ}V_{uR}^{Qj}V_{dR}^{\dagger jP}\right)^{\ast}=V_{dR}^{Pj}V_{uR}^{\dagger
jQ}V_{uR}^{Qi}V_{dR}^{\dagger iP}$ and it corresponds to Eq. (4.3) by an
exchange of $i\leftrightarrow j$, the contribution is real so that it does not
generate the QCD $\bar{\theta}$ term.

The reason why the exchange of the $W^{\prime\,\pm}$ boson does not contribute
to the QCD $\bar{\theta}$ term at two-loop level is clear. However, the above
discussion is based on the structure of the mixing matrices in the
contribution, not on the structure of Lagrangian parameters, such as $x_{u/d}$
and $M_{u/d}$, and then, it is unclear what is required to generate the QCD
$\bar{\theta}$ term in higher-order diagrams. We make it clear by explicit
calculation of the loop diagrams in the following.

In the unitary gauge, the lowest dimension operator ($n=1$) in Eq. (4.1) might
come from diagrams which include the longitudinal mode of the
$W^{\prime\,\pm}$ boson. It is because the propagator is proportional to
$k^{\mu}k^{\nu}/m_{W^{\prime}}^{2}$ ($k^{\nu}$ the momentum of the
$W^{\prime\,\pm}$ boson), and it could give the lowest order contribution with
regard to $v^{\prime\,2}$. The Yukawa coupling constants $x_{u}$ and $x_{d}$
are multiplied with the Higgs VEVs in the mixing matrices as in Eq. (3.11).

In our calculation, we adopt the $R_{\xi}$ gauge with the Feynman-’t Hooft
gauge $\xi=1$ (for $SU(2)_{L}\times SU(2)_{R}\times U(1)_{B-L}$ gauge), in
order to avoid the messy calculation in the unitary gauge. The lowest
dimension operator ($n=1$) in Eq. (4.1) could arise the charged NG boson
exchange $\varphi^{\prime\,\pm}$ in this gauge. The charged NG boson is
absorbed by $W^{\prime\,\pm}$ boson in the Higgs mechanism, and its mass,
$m_{\varphi^{\prime}}$, is equal to the $W^{\prime\,\pm}$ mass. The charged NG
boson interactions are given as

$\displaystyle-{\cal L}_{\varphi^{\prime\pm}}$ $\displaystyle=$
$\displaystyle\overline{u}_{R}^{i}x_{d}^{ia}D^{a}_{L}\varphi^{\prime+}-\overline{d}_{R}^{i}x_{u}^{ia}U^{a}_{L}\varphi^{\prime\,-}+{\rm
h.c.}$ (4.4) $\displaystyle=$ $\displaystyle(x_{d}^{ia}V_{uR}^{Pi}V_{dL}^{\ast
Qa})\,\overline{\cal U}_{MR}^{P}{\cal
D}^{Q}_{ML}\varphi^{\prime+}-(x_{u}^{\ast ia}V_{uL}^{Pa}V_{dR}^{\ast
Qi})\,\overline{\cal U}_{ML}^{P}{\cal D}^{Q}_{MR}\varphi^{\prime+}+{\rm
h.c.}\,.$

Both the left- and right-handed quarks are coupled with the charged NG boson.
We will show that the charged NG boson diagrams at two-loop level do not
contribute to the QCD $\bar{\theta}$ term.

(a) diagram $A$ (b) diagram $B$ (c) diagram $C$ Figure 2: The charged NG boson
contributions to the QCD $\bar{\theta}$ term at two-loop level.

Three diagrams (dubbed as diagrams $A$, $B$ and $C$) in Fig. 2 could give
contributions to the $\bar{\theta}$ parameter. By using the Fock-Schwinger
gauge method (for $SU(3)_{C}$ gauge) in Sec. 2.2, we obtain

$\displaystyle\delta\theta|_{A}$ $\displaystyle=\frac{1}{8\pi^{2}}{\rm
Im}\left(x_{u}^{ia}V_{uL}^{\ast
Pa}V_{dR}^{Qi}x_{d}^{jb}V_{uR}^{Pj}V_{dL}^{\ast
Qb}\right)\bar{M}_{u}^{P}\bar{M}_{d}^{Q}\,\bar{I}_{(1;3)}\left[(\bar{M}_{u}^{P})^{2};(\bar{M}_{d}^{Q})^{2};m_{\varphi^{\prime}}^{2}\right]\,,$
(4.5) $\displaystyle\delta\theta|_{B}$ $\displaystyle=\frac{1}{8\pi^{2}}{\rm
Im}\left(x_{u}^{ia}V_{uL}^{\ast
Pa}V_{dR}^{Qi}x_{d}^{jb}V_{uR}^{Pj}V_{dL}^{\ast
Qb}\right)\bar{M}_{u}^{P}\bar{M}_{d}^{Q}\,\bar{I}_{(3;1)}\left[(\bar{M}_{u}^{P})^{2};(\bar{M}_{d}^{Q})^{2};m_{\varphi^{\prime}}^{2}\right]\,,$
(4.6) $\displaystyle\delta\theta|_{C}$ $\displaystyle=\frac{1}{8\pi^{2}}{\rm
Im}\left(x_{u}^{ia}V_{uL}^{\ast
Pa}V_{dR}^{Qi}x_{d}^{jb}V_{uR}^{Pj}V_{dL}^{\ast
Qb}\right)\bar{M}_{u}^{P}\bar{M}_{d}^{Q}\,I_{(2;2)}\left[(\bar{M}_{u}^{P})^{2};(\bar{M}_{d}^{Q})^{2};m_{\varphi^{\prime}}^{2}\right]\,,$
(4.7)

where the two-loop functions $\bar{I}_{(1;3)}$, $\bar{I}_{(3;1)}$, and
$I_{(2;2)}$ are defined in Appendix A. In the above evaluation, we pick up
contributions proportional to both $x_{u}$ and $x_{d}$, which are also
proportional to the quark masses in the mass eigenstate propagators. The terms
proportional to $x_{u}$ and $x_{u}^{\ast}$ (or $x_{d}$ and $x_{d}^{\ast}$) is
real.

Here, we use the dimensional regularization $(d=4-2\epsilon)$ for loop
momentum integrals and the partial diagrams produce UV divergence. However,
the contributions from diagrams $A$ and $B$, proportional to
$\bar{I}_{\epsilon{(1;3)}}$ and $\bar{I}_{\epsilon{(3;1)}}$ in Eq. (A.13),
respectively, vanish, so that the correction to the $\bar{\theta}$ parameter
is finite and scale-independent. The UV divergent parts ($1/\epsilon$ terms)
in $\bar{I}_{\epsilon{(1;3)}}$ and $\bar{I}_{\epsilon{(3;1)}}$ cancel out
since

$\displaystyle{\rm Im}\left(x_{u}^{ia}V_{uL}^{\ast
Pa}V_{dR}^{Qi}x_{d}^{jb}V_{uR}^{Pj}V_{dL}^{\ast
Qb}\right)\frac{\bar{M}_{u}^{P}}{\bar{M}_{d}^{Q}}$ $\displaystyle=$
$\displaystyle{\rm Im}\left(x_{u}^{\ast ia}x_{u}^{ia}\right)=0\,,$ (4.8)
$\displaystyle{\rm Im}\left(x_{u}^{ia}V_{uL}^{\ast
Pa}V_{dR}^{Qi}x_{d}^{jb}V_{uR}^{Pj}V_{dL}^{\ast
Qb}\right)\frac{\bar{M}_{d}^{Q}}{\bar{M}_{u}^{P}}$ $\displaystyle=$
$\displaystyle{\rm Im}\left(x_{d}^{\ast ia}x_{d}^{ia}\right)=0\,.$ (4.9)

To derive the above equations we use the following equations,

$\displaystyle[({\cal M}_{q}^{(0)})^{-1}]^{pq}$ $\displaystyle=$
$\displaystyle V_{qR}^{\dagger pP}(\bar{M}_{q}^{-1})^{P}V_{qL}^{Pq}$ (4.10)
$\displaystyle=$
$\displaystyle\begin{pmatrix}-\frac{1}{vv^{\prime}}(x_{q}^{\dagger})^{-1}_{ic}M^{c}_{q}(x_{q})^{-1}_{cj}&\frac{1}{v^{\prime}}(x_{q}^{\dagger})^{-1}_{ib}\\\
\frac{1}{v}(x_{q})^{-1}_{aj}&0\end{pmatrix}\quad\text{for}\leavevmode\nobreak\
q=u\text{\leavevmode\nobreak\ and\leavevmode\nobreak\ }d\,,$

in addition to Eq. (3.6) with $(\bar{M}_{q}^{P})^{\ast}=\bar{M}_{q}^{P}$.
Furthermore, we observed that the following combinations also
vanish#12#12#12The factors $1/M\log M$ ($M$: quark mass) in Eqs. (4.11) and
(4.12) correspond to the $O(\epsilon)$ term in Eq. (2.33) when changing
$d^{4}p$ to $d^{d}p$ $(d=4-2\epsilon)$. They become $O(\epsilon^{0})$ since
$1/\epsilon$ comes from the quark self-energy subdiagrams in diagrams $A$ and
$B$. Equations (4.11) and (4.12) can be perturbatively proved by assuming the
off-diagonal terms in the quark mass matrices are small. The similar trick is
also used around Eq. (4.26). We also checked Eqs. (4.11) and (4.12)
numerically.

$\displaystyle{\rm Im}\left(x_{u}^{ia}V_{uL}^{\ast
Pa}V_{dR}^{Qi}x_{d}^{jb}V_{uR}^{Pj}V_{dL}^{\ast
Qb}\right)\frac{\bar{M}_{u}^{P}}{\bar{M}_{d}^{Q}}\log\bar{M}_{d}^{Q}$
$\displaystyle=$ $\displaystyle 0\,,$ (4.11) $\displaystyle{\rm
Im}\left(x_{u}^{ia}V_{uL}^{\ast
Pa}V_{dR}^{Qi}x_{d}^{jb}V_{uR}^{Pj}V_{dL}^{\ast
Qb}\right)\frac{\bar{M}_{d}^{Q}}{\bar{M}_{u}^{P}}\log\bar{M}_{u}^{P}$
$\displaystyle=$ $\displaystyle 0\,.$ (4.12)

Therefore, the second terms of $\bar{I}_{\epsilon{(1;3)}}$ and
$\bar{I}_{\epsilon{(3;1)}}$ in Eq. (A.13) also do not affect the
$\bar{\theta}$ parameter. On the other hand, the loop function $I_{(2;2)}$ in
Eq. (4.7) is UV finite.

Similar to the $\varphi^{\prime\,\pm}$ contribution, the contribution from the
SM charged NG boson $\varphi^{\pm}$, absorbed into $W^{\pm}$, is derived by
replacing $L(R)$ with $R(L)$, and $m_{\varphi^{\prime}}^{2}$ with
$m_{\varphi}^{2}$ in the above formulae. Furthermore, $(-1)$ is multiplied
since the chiralities of circulating fermions are opposite to diagrams of Fig.
2. It means that, if one sets $v=v^{\prime}$ corresponding to the LR symmetric
limit, those two contributions of $\varphi^{\prime\,\pm}$ and $\varphi^{\pm}$
cancel each other.

The diagrams $A$ and $B$ correspond to the one-loop correction to the fermion
mass terms. The two-loop function $\bar{I}_{(3;1)}(x_{1};x_{2};x_{3})$ is
expressed as

$\displaystyle\bar{I}_{(3;1)}(x_{1};x_{2};x_{3})$ $\displaystyle=$
$\displaystyle(16\pi^{2}\mu^{2\epsilon})\int\frac{d^{d}p}{i(2\pi)^{d}}\frac{1}{(p^{2}-x_{1})^{3}}F_{0}(p^{2},x_{2},x_{3})\,,$
(4.13)

where $F_{0}(p^{2},x_{2},x_{3})$ is a loop function of one-loop diagrams for
the fermion mass correction,

$\displaystyle F_{0}(p^{2},x_{2},x_{3})$ $\displaystyle=$
$\displaystyle\int^{1}_{0}dz\log\frac{-z(1-z)p^{2}+zx_{2}+(1-z)x_{3}}{Q^{2}}\,,$
(4.14)

with $Q^{2}\equiv 4\pi\mu^{2}{\rm e}^{-\gamma_{E}}$ and $\mu$ is the
renormalization scale. ($\bar{I}_{(1;3)}(x_{1};x_{2};x_{3})$ also has a
similar expression, see Appendix A.) $\bar{I}_{(3;1)}(x_{1};x_{2};x_{3})$ has
an IR-singular behavior when $x_{1}\ll x_{3}$ as

$\displaystyle\bar{I}_{(3;1)}(x_{1};x_{2};x_{3})$ $\displaystyle\simeq$
$\displaystyle-\frac{1}{2x_{1}}F_{0}(0,x_{2},x_{3})\,,$ (4.15)

while small $x_{2}$ and $x_{3}$ do not lead to IR singularities. This behavior
is expected. It is because if a fermion with real mass $m_{f}$ gets a constant
radiative correction to the fermion mass $m_{f}+\delta m_{f}$, the correction
to the $\bar{\theta}$ parameter is given by $\delta\theta\simeq-{\rm
Im}(\delta m_{f})/m_{f}$, see Eq. (2.33).#13#13#13 In Eq. (2.33), the loop
function and the chirality flip lead to $\sim m_{f}/m_{f}^{2}=1/m_{f}$, and
then $\delta\theta\simeq-{\rm Im}(\delta m_{f})/m_{f}$.  However, this
evaluation of $\delta\theta$ is justified only when the correction to the
fermion mass term is independent of fermion momentum.

The IR-singular behaviors of the SM fermion masses in $\bar{I}_{(3;1)}$ and
$\bar{I}_{(1;3)}$ are not physical in $\delta\theta|_{A}$ and
$\delta\theta|_{B}$ in Eqs. (4.5) and (4.6), and they can be removed indeed
using Eq. (4.10) as

$\displaystyle\begin{aligned} \delta\theta|_{A}\simeq&\,\frac{1}{8\pi^{2}}{\rm
Im}\left[x_{u}^{ia}x_{d}^{jb}\left(V_{uL}^{\ast
Aa}\bar{M}_{u}^{A}V_{uR}^{Aj}\right)\left(V_{dR}^{Bi}\frac{1}{\bar{M}_{d}^{B}}V^{\ast
Bb}_{dL}\right)\right]\\\
&\times\left[(\bar{M}_{d}^{B})^{2}\bar{I}_{(1;3)}\left((\bar{M}_{u}^{A})^{2};(\bar{M}_{d}^{B})^{2};m^{2}_{\varphi^{\prime}}\right)-(\bar{M}_{d}^{B})^{2}\bar{I}_{(1;3)}\left(0;(\bar{M}_{d}^{B})^{2};m^{2}_{\varphi^{\prime}}\right)\right.\\\
&\qquad\left.+\frac{1}{2}F_{0}\left(0,(\bar{M}_{u}^{A})^{2},m^{2}_{\varphi^{\prime}}\right)-\frac{1}{2}F_{0}\left(0,0,m^{2}_{\varphi^{\prime}}\right)\right]\\\
&+\frac{1}{8\pi^{2}}{\rm Im}\left[x_{u}^{ia}x_{u}^{\ast
ja}x_{d}^{jb}v^{\prime}\left(V_{dR}^{Ai}\frac{1}{\bar{M}_{d}^{A}}V^{\ast
Ab}_{dL}\right)\right]\\\
&\times\left[(\bar{M}_{d}^{A})^{2}\bar{I}_{(1;3)}\left(0;(\bar{M}_{d}^{A})^{2};m^{2}_{\varphi^{\prime}}\right)+\frac{1}{2}F_{0}\left(0,0,m^{2}_{\varphi^{\prime}}\right)\right]\,,\end{aligned}$
(4.16) $\displaystyle\begin{aligned}
\delta\theta|_{B}\simeq&\,\frac{1}{8\pi^{2}}{\rm
Im}\left[x_{u}^{ia}x_{d}^{jb}\left(V_{uR}^{Aj}\frac{1}{\bar{M}_{u}^{A}}V^{\ast
Aa}_{uL}\right)\left(V^{\ast
Bb}_{dL}\bar{M}_{d}^{B}V_{dR}^{Bi}\right)\right]\\\
&\times\left[(\bar{M}_{u}^{A})^{2}\bar{I}_{(3;1)}\left((\bar{M}_{u}^{A})^{2};(\bar{M}_{d}^{B})^{2};m^{2}_{\varphi^{\prime}}\right)-(\bar{M}_{u}^{A})^{2}\bar{I}_{(3;1)}\left((\bar{M}_{u}^{A})^{2};0;m^{2}_{\varphi^{\prime}}\right)\right.\\\
&\qquad\left.+\frac{1}{2}F_{0}\left(0,(\bar{M}_{d}^{B})^{2},m^{2}_{\varphi^{\prime}}\right)-\frac{1}{2}F_{0}\left(0,0,m^{2}_{\varphi^{\prime}}\right)\right]\\\
&+\frac{1}{8\pi^{2}}{\rm Im}\left[x_{u}^{ia}x_{d}^{jb}x_{d}^{\ast
ib}v^{\prime}\left(V_{uR}^{Aj}\frac{1}{\bar{M}_{u}^{A}}V_{uL}^{{}^{\ast}Aa}\right)\right]\\\
&\times\left[(\bar{M}_{u}^{A})^{2}\bar{I}_{(3;1)}\left((\bar{M}_{u}^{A})^{2};0;m^{2}_{\varphi^{\prime}}\right)+\frac{1}{2}F_{0}\left(0,0,m^{2}_{\varphi^{\prime}}\right)\right]\,,\end{aligned}$
(4.17)

where $A$ and $B$ run $4$–$6$ as the heavy quark mass eigenstates, see Eq.
(3.12). Here, the SM quark masses in the loop function are taken to be zero.

On the other hand, the contribution of diagram $C$ is not associated with the
correction to the quark masses, and then it is a new type contribution to the
$\bar{\theta}$ parameter. It is suppressed by the heavier fermion or
$\varphi^{\prime\,\pm}$ masses. By taking the SM quark masses to be zero in
the loop function (see Appendix A), it is given as

$\displaystyle\delta\theta|_{C}\simeq\frac{1}{8\pi^{2}}{\rm
Im}\left[x_{u}^{ia}x_{d}^{jb}\left(V_{uL}^{\ast
Aa}\bar{M}_{u}^{A}V_{uR}^{Aj}\right)\left(V_{dL}^{\ast
Bb}\bar{M}_{d}^{B}V_{dR}^{Bi}\right)\right]{I}_{(2;2)}\left((\bar{M}_{u}^{A})^{2};(\bar{M}_{d}^{B})^{2};m^{2}_{\varphi^{\prime}}\right)\,,$
(4.18)

where $A$ and $B$ run $4$–$6$ as the heavy quark mass eigenstates. It is found
that the diagram $C$ does not contain any IR-singular behavior, unlike the
diagrams $A$ and $B$.

When $v^{\prime}\mathop{}_{\textstyle\sim}^{\textstyle<}M_{q}^{a}$, the
leading contributions of $\mathcal{O}(v^{\prime 2}/(M_{q}^{a})^{2})$ are given
as

$\displaystyle\delta\theta|_{A}$ $\displaystyle\simeq\frac{1}{8\pi^{2}}{\rm
Im}\left[(A_{u}^{a})^{ij}(A_{d}^{b})^{ji}\right]\left\\{v^{\prime
2}\,\bar{I}_{(1;3)}\left[(M_{u}^{a})^{2};(M_{d}^{b})^{2};m^{2}_{\varphi^{\prime}}\right]+\frac{v^{\prime
2}}{2(M_{d}^{b})^{2}}F_{0}\left[0,(M_{u}^{a})^{2},m^{2}_{\varphi^{\prime}}\right]\right\\}\,,$
(4.19) $\displaystyle\delta\theta|_{B}$
$\displaystyle\simeq\frac{1}{8\pi^{2}}{\rm
Im}\left[(A_{u}^{a})^{ij}(A_{d}^{b})^{ji}\right]\left\\{v^{\prime
2}\,\bar{I}_{(3;1)}\left[(M_{u}^{a})^{2};(M_{d}^{b})^{2};m^{2}_{\varphi^{\prime}}\right]+\frac{v^{\prime
2}}{2(M_{u}^{a})^{2}}F_{0}\left[0,(M_{d}^{b})^{2},m^{2}_{\varphi^{\prime}}\right]\right\\}\,,$
(4.20) $\displaystyle\delta\theta|_{C}$
$\displaystyle\simeq\frac{1}{8\pi^{2}}{\rm
Im}\left[(A_{u}^{a})^{ij}(A_{d}^{b})^{ji}\right]v^{\prime
2}\,{I}_{(2;2)}\left[(M_{u}^{a})^{2};(M_{d}^{b})^{2};m^{2}_{\varphi^{\prime}}\right]\,,$
(4.21)

with a Hermitian matrix $A_{q}^{a},$

$\displaystyle(A_{q}^{a})^{ij}\equiv x_{q}^{ia}x_{q}^{\ast
ja}\quad\text{for}\leavevmode\nobreak\ q=u,d\text{ and\leavevmode\nobreak\
not\leavevmode\nobreak\ sum\leavevmode\nobreak\ }a\text{\leavevmode\nobreak\
index}\,.$ (4.22)

This can be derived from the above formulae by

$\displaystyle(V^{\dagger}_{qL})^{aA}\frac{\bar{M}_{q}^{A}}{p^{2}-(\bar{M}_{q}^{A})^{2}}(V_{qR})^{Ai}$
$\displaystyle\rightarrow$ $\displaystyle\frac{x_{q}^{\dagger
ai}v^{\prime}}{p^{2}-(M_{q}^{a})^{2}}\,,$ (4.23)
$\displaystyle(V^{\dagger}_{qL})^{aA}\frac{\bar{M}_{q}^{A}}{[p^{2}-(\bar{M}_{q}^{A})^{2}]^{2}}(V_{qR})^{Ai}$
$\displaystyle\rightarrow$ $\displaystyle\frac{x_{q}^{\dagger
ai}v^{\prime}}{[p^{2}-(M_{q}^{a})^{2}]^{2}}\,,$ (4.24)
$\displaystyle(V^{\dagger}_{qL})^{aA}\frac{\bar{M}_{q}^{A}}{[p^{2}-(\bar{M}_{q}^{A})^{2}]^{3}}(V_{qR})^{Ai}$
$\displaystyle\rightarrow$ $\displaystyle\frac{x_{q}^{\dagger
ai}v^{\prime}}{[p^{2}-(M_{q}^{a})^{2}]^{3}}\,.$ (4.25)

It is found that these radiative corrections to the $\bar{\theta}$ parameter
vanish. For example, $\delta\theta|_{A}$ is given as

$\displaystyle\delta\theta|_{A}$ $\displaystyle=$ $\displaystyle{\rm
Im\,Tr}\left(A^{a}_{u}A^{b}_{d}\right)f\left[(M^{a}_{u})^{2},(M^{b}_{d})^{2},m_{\varphi^{\prime}}^{2}\right]$
(4.26) $\displaystyle=$ $\displaystyle\frac{1}{2}{\rm
Im\,Tr}\left([A^{a}_{u},A^{b}_{d}]\right)f\left[(M^{a}_{u})^{2},(M^{b}_{d})^{2},m_{\varphi^{\prime}}^{2}\right]=0\,,$

where $f[(M_{u}^{a})^{2},(M_{d}^{b})^{2},m_{\varphi^{\prime}}^{2}]$ is the
real function, and the Hermitian property of the matrix $A_{q}^{a}$ is used.
The same conclusions are applicable to $\delta\theta|_{B}$ and
$\delta\theta|_{C}$ at this order.

Now we showed that the charged NG boson contribution to the $\bar{\theta}$
parameter at two-loop level vanishes at the leading order of $v^{\prime}$
($n=1$ in Eq. (4.1)). It comes from the fact that the contributions are
proportional to the fourth power of $x_{u/d}$. We have also checked that the
contributions of the sixth power of $x_{u/d}$, corresponding to
$\mathcal{O}(v^{\prime 4}/(M_{q}^{a})^{4})$ contributions, also vanish. The
contributions are derived from the above formulae with the mass-insertion
approximation,

$\displaystyle(V^{\dagger}_{qL})^{aA}\frac{i\bar{M}_{q}^{A}}{p^{2}-(\bar{M}_{q}^{A})^{2}}(V_{qR})^{Ai}$
$\displaystyle\rightarrow$ $\displaystyle i\frac{x_{q}^{\dagger
ai}v^{\prime}}{p^{2}-(M_{q}^{a})^{2}}$ (4.27)
$\displaystyle+i\frac{1}{p^{2}-(M_{q}^{a})^{2}}(x_{q}^{\dagger}x_{q})^{ab}v^{\prime
2}\frac{x_{q}^{\dagger bi}v^{\prime}}{p^{2}-(M_{q}^{b})^{2}}\,,$
$\displaystyle(V^{\dagger}_{qL})^{aA}\frac{i\bar{M}_{q}^{A}}{[p^{2}-(\bar{M}_{q}^{A})^{2}]^{2}}(V_{qR})^{Ai}$
$\displaystyle\rightarrow$ $\displaystyle i\frac{x_{q}^{\dagger
ai}v^{\prime}}{[p^{2}-(M_{q}^{a})^{2}]^{2}}$ (4.28)
$\displaystyle+i\frac{1}{p^{2}-(M_{q}^{a})^{2}}(x_{q}^{\dagger}x_{q})^{ab}v^{\prime
2}\frac{x_{q}^{\dagger bi}v^{\prime}}{[p^{2}-(M_{q}^{b})^{2}]^{2}}$
$\displaystyle+i\frac{1}{[p^{2}-(M_{q}^{a})^{2}]^{2}}(x_{q}^{\dagger}x_{q})^{ab}v^{\prime
2}\frac{x_{q}^{\dagger bi}v^{\prime}}{p^{2}-(M_{q}^{b})^{2}}\,,$
$\displaystyle(V^{\dagger}_{qL})^{aA}\frac{i\bar{M}_{q}^{A}}{[p^{2}-(\bar{M}_{q}^{A})^{2}]^{3}}(V_{qR})^{Ai}$
$\displaystyle\rightarrow$ $\displaystyle i\frac{x_{q}^{\dagger
ai}v^{\prime}}{[p^{2}-(M_{q}^{a})^{2}]^{3}}$ (4.29)
$\displaystyle+i\frac{1}{p^{2}-(M_{q}^{a})^{2}}(x_{q}^{\dagger}x_{q})^{ab}v^{\prime
2}\frac{x_{q}^{\dagger bi}v^{\prime}}{[p^{2}-(M_{q}^{b})^{2}]^{3}}$
$\displaystyle+i\frac{1}{[p^{2}-(M_{q}^{a})^{2}]^{2}}(x_{q}^{\dagger}x_{q})^{ab}v^{\prime
2}\frac{x_{q}^{\dagger bi}v^{\prime}}{[p^{2}-(M_{q}^{b})^{2}]^{2}}$
$\displaystyle+i\frac{1}{[p^{2}-(M_{q}^{a})^{2}]^{3}}(x_{q}^{\dagger}x_{q})^{ab}v^{\prime
2}\frac{x_{q}^{\dagger bi}v^{\prime}}{p^{2}-(M_{q}^{b})^{2}}\,.$

Each first term is the aforementioned leading contribution, which vanishes
(see Eq. (4.26)). The next-to-leading contributions are

$\displaystyle\delta\theta|_{A}\simeq$ $\displaystyle\frac{1}{8\pi^{2}}{\rm
Im}\left[(A_{u}^{a})^{ij}(A_{u}^{b})^{jk}(A_{d}^{c})^{ki}\right]$
$\displaystyle\times\left\\{v^{\prime
4}I_{(1,1;3)}\left[(M^{a}_{u})^{2},(M^{b}_{u})^{2};(M^{c}_{d})^{2};m^{2}_{\varphi^{\prime}}\right]-\frac{v^{\prime
4}}{(M^{c}_{d})^{2}}B_{(1,1)}\left[0,(M^{a}_{u})^{2},(M^{b}_{u})^{2};m^{2}_{\varphi^{\prime}}\right]\right\\}$
$\displaystyle+\frac{1}{8\pi^{2}}{\rm
Im}\left[(A_{u}^{a})^{ij}(A_{d}^{b})^{jk}(A_{d}^{c})^{ki}\right]v^{\prime
4}\left(I_{(1;1,3)}+I_{(1;2,2)}+I_{(1;3,1)}\right)\left[(M^{a}_{u})^{2};(M^{b}_{d})^{2},(M^{c}_{d})^{2};m^{2}_{\varphi^{\prime}}\right]\,,$
(4.30) $\displaystyle\delta\theta|_{B}\simeq$
$\displaystyle\frac{1}{8\pi^{2}}{\rm
Im}\left[(A_{u}^{a})^{ij}(A_{d}^{b})^{jk}(A_{d}^{c})^{ki}\right]$
$\displaystyle\times\left\\{v^{\prime
4}I_{(3;1,1)}\left[(M^{a}_{u})^{2};(M^{b}_{d})^{2},(M^{c}_{d})^{2};m^{2}_{\varphi^{\prime}}\right]-\frac{v^{\prime
4}}{(M^{a}_{u})^{2}}B_{(1,1)}\left[0,(M^{b}_{d})^{2},(M^{c}_{d})^{2};m^{2}_{\varphi^{\prime}}\right]\right\\}$
$\displaystyle+\frac{1}{8\pi^{2}}{\rm
Im}\left[(A_{u}^{a})^{ij}(A_{u}^{b})^{jk}(A_{d}^{c})^{ki}\right]v^{\prime
4}\left(I_{(1,3;1)}+I_{(2,2;1)}+I_{(3,1;1)}\right)\left[(M^{a}_{u})^{2},(M^{b}_{u})^{2};(M^{c}_{d})^{2};m^{2}_{\varphi^{\prime}}\right]\,,$
(4.31) $\displaystyle\delta\theta|_{C}\simeq$
$\displaystyle\frac{1}{8\pi^{2}}{\rm
Im}\left[(A_{u}^{a})^{ij}(A_{u}^{b})^{jk}(A_{d}^{c})^{ki}\right]v^{\prime
4}\left({I}_{(1,2;2)}+{I}_{(2,1;2)}\right)\left[(M_{u}^{a})^{2},(M_{u}^{b})^{2};(M_{d}^{c})^{2};m^{2}_{\varphi^{\prime}}\right]$
$\displaystyle+\frac{1}{8\pi^{2}}{\rm
Im}\left[(A_{u}^{a})^{ij}(A_{d}^{b})^{jk}(A_{d}^{c})^{ki}\right]v^{\prime
4}\left({I}_{(2:1,2)}+{I}_{(2:2,1)}\right)\left[(M_{u}^{a})^{2};(M_{d}^{b})^{2},(M_{d}^{c})^{2};m^{2}_{\varphi^{\prime}}\right]\,.$
(4.32)

These loop functions, which come from the mass-insertion approximation, are
given in Appendix A. The sequences of masses connected by commas in the loop
functions are introduced by the mass insertion.

Again, it is found that these contributions are zero. For example, the first
term in $\delta\theta|_{A}$ is given as

$\displaystyle\delta\theta|_{A}$ $\displaystyle=$ $\displaystyle{\rm
Im\,Tr}\left(A^{a}_{u}A^{b}_{u}A^{c}_{d}\right)f\left[(M^{a}_{u})^{2},(M^{b}_{u})^{2};(M^{c}_{d})^{2};m_{\varphi^{\prime}}^{2}\right]$
(4.33) $\displaystyle+{\rm
Im\,Tr}\left(A^{a}_{u}A^{b}_{d}A^{c}_{d}\right)g\left[(M^{a}_{u})^{2};(M^{b}_{d})^{2},(M^{c}_{d})^{2};m_{\varphi^{\prime}}^{2}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{2}{\rm
Im\,Tr}\left(A^{c}_{d}\,[A^{a}_{u},A^{b}_{u}]\right)f\left[(M^{a}_{u})^{2},(M^{b}_{u})^{2};(M^{c}_{d})^{2};m_{\varphi^{\prime}}^{2}\right]$
$\displaystyle+\frac{1}{2}{\rm
Im\,Tr}\left(A^{a}_{u}\,[A^{b}_{d},A^{c}_{d}]\right)g\left[(M^{a}_{u})^{2};(M^{b}_{d})^{2},(M^{c}_{d})^{2};m_{\varphi^{\prime}}^{2}\right]$
$\displaystyle=$ $\displaystyle 0\,.$

Here, above two real loop functions $f$ and $g$ are symmetric under exchanges
of $(M^{a}_{u})^{2}\leftrightarrow(M^{b}_{u})^{2}$ and
$(M^{b}_{d})^{2}\leftrightarrow(M^{c}_{d})^{2}$, respectively. Then, the above
equation vanishes, see Appendix A. The symmetry comes from the mass-insertion
approximation. Even if we include the higher-order contributions of $x_{q}$ in
the mass-insertion approximation, they still vanish since the loop function is
real and symmetric for the exchange of the heavy fermion masses.

Now we found that the charged NG boson does not give a contribution to the
$\bar{\theta}$ parameter at two-loop level. We also numerically checked this
fact by using Eqs. (4.16)–(4.18).

The $W^{\prime\,\pm}$ contributions to the $\bar{\theta}$ parameter in the
Feynman-’t Hooft gauge at two-loop level vanish. The Yukawa coupling
dependence comes from only the mixing matrices, and then the leading
contributions, which are proportional to the fourth power of $x_{u/d}$ at
most, vanish. The higher-order contributions, coming from mass-insertion
approximation, also vanish due to the symmetry of heavy fermion masses in loop
functions.

A similar discussion is applicable for the other contribution, such as
$Z^{\prime}$, $h^{\prime}$, and $\varphi^{\prime\,0}$ at two-loop level. Then,
we confirmed the two-loop contribution to the $\bar{\theta}$ parameter
vanishes as far as $M_{q}\,\mathop{}_{\textstyle\sim}^{\textstyle>}\langle
H^{\prime}\rangle\,\gg\langle H\rangle$.

## 5 Non-vanishing contribution to QCD $\theta$ parameter in three-loop order

In the previous section, we confirmed that the QCD $\bar{\theta}$ term is not
generated in the two-loop level contribution, i.e., up to the fourth order of
the Yukawa interaction $x_{q}$. We also found that it is valid even if one
considers the higher-order contributions of $x_{q}$ by using the mass-
insertion approximation. In order to give non-vanishing contributions to the
$\bar{\theta}$ parameter, the commutation relation $[A^{a}_{q},A^{b}_{q}]$
must be nonzero, see Eq. (4.33). It implies that non-vanishing contribution
should be proportional to ${\rm
Im\,Tr}(A^{a}_{q^{\prime}}\,[A^{b}_{q},A^{c}_{q}])$ for $q,q^{\prime}=u$
and/or $d$ rather than ${\rm Im\,Tr}([A^{b}_{q},A^{c}_{q}])$, and the loop
function has to be asymmetric under exchange between $(M_{q}^{b})^{2}$ and
$(M_{q}^{c})^{2}$. Thus, the contributions of the following form might be
leading if they are non-vanishing in the three-loop order,

$\displaystyle\delta\theta_{uuu}$
$\displaystyle\approx\frac{1}{(16\pi^{2})^{2}}\frac{v^{\prime
2}}{\widetilde{M}^{2}}\,{\rm
Im\,Tr}\left(A^{a}_{u}\,[A^{b}_{u},A^{c}_{u}]\right)f_{uuu}^{abc}\,,$ (5.1)
$\displaystyle\delta\theta_{duu}$
$\displaystyle\approx\frac{1}{(16\pi^{2})^{2}}\frac{v^{\prime
2}}{\widetilde{M}^{2}}{\rm
Im\,Tr}\left(A^{a}_{d}\,[A^{b}_{u},A^{c}_{u}]\right)f_{duu}^{abc}\,,$ (5.2)

where $\widetilde{M}$ is the heaviest quark mass in the loop diagrams. Here,
$f_{uuu}^{abc}$ is a dimensionless three-loop function which is totally
antisymmetric under permutation of $(M_{u}^{a})^{2},(M_{u}^{b})^{2}$ and
$(M_{u}^{c})^{2}$, while $f_{duu}^{abc}$ is antisymmetric under permutation of
$(M_{u}^{b})^{2}$ and $(M_{u}^{c})^{2}$. Although the other types such as
${\rm Im\,Tr}(A_{u}^{a}\,[A_{d}^{b},A_{d}^{c}])$ and ${\rm
Im\,Tr}(A_{d}^{a}\,[A_{d}^{b},A_{d}^{c}])$ can also contribute to the
$\bar{\theta}$ parameter, they are suppressed by the SM down-type quark
masses, so we do not take them into account in this paper.

We found that $\delta\theta_{uuu}$ in Eq. (5.1) is not generated from diagrams
in the three-loop order in the minimal LR model. The corresponding three-loop
diagrams do not have an asymmetric structure for three Dirac fermion masses
when the internal scalar lines respect the $U(1)_{B-L}$. If the neutral scalar
lines break the $U(1)_{B-L}$ (or the scalar lines pick $v^{\prime}$ using the
four-point Higgs interaction), the diagrams may have the asymmetric structure.
However, in the case, the contribution the $\bar{\theta}$ parameter is
proportional to $v^{\prime 4}/{\widetilde{M}}^{4}$ ($n=2$ in Eq. (4.1)), not
$v^{\prime 2}/{\widetilde{M}}^{2}$. This situation is not changed even in the
four-loop order. Thus, we conclude that $\delta\theta_{uuu}$ is not the
leading contribution and the largest non-vanishing contribution to the
$\bar{\theta}$ parameter comes from $\delta\theta_{duu}$ in Eq. (5.2) in the
minimal LR symmetric model.

### 5.1 Leading contribution: probability density function of
$\delta\theta_{duu}$

Figure 3: Diagrams that contribute to Eq. (5.2) and have asymmetric structure
in the loop function. The neutral scalar propagator corresponds to the neutral
NG boson $\varphi^{\prime\,0}$ and the neutral Higgs boson $h^{\prime}$.

In this section, we estimate the size of $\delta\theta_{duu}$ in Eq. (5.2). We
find that diagrams in Fig. 3 would provide the antisymmetric three-loop
function and produce the non-vanishing $\delta\theta_{duu}$.

When one considers a universal down-type vector-like quark mass
$M_{d}^{1}=M_{d}^{2}=M_{d}^{3}\equiv\widetilde{M}_{d}$ for simplicity,
$\overline{\Phi}(\theta_{d3},\theta_{d8})V_{D}$ in the down-type Yukawa
$x_{d}$ in Eq. (3.10) become unphysical parameters, because these are removed
by changing the basis of $D_{L/R}$. Then, the contribution from
$\delta\theta_{duu}$ is simplified as

$\displaystyle\delta\theta_{duu}\approx$
$\displaystyle\frac{1}{(16\pi^{2})^{2}}\frac{v^{\prime
2}}{\widetilde{M}^{2}}\frac{\widetilde{M}_{d}M_{u}^{b}M_{u}^{c}}{v^{\prime
3}}\frac{m_{d}^{i}m_{u}^{k}\sqrt{m_{u}^{j}m_{u}^{l}}}{v^{3}}$
$\displaystyle\times{\rm Im\,Tr}\left[V_{\rm CKM}^{\dagger ij}V_{U}^{\prime
jb}V_{U}^{\prime\dagger bk}V_{U}^{\prime kc}V_{U}^{\prime\dagger cl}V_{\rm
CKM}^{li}-{(b\leftrightarrow c)}\right]\tilde{f}_{duu}^{bc}$
$\displaystyle\approx$ $\displaystyle\frac{4}{(16\pi^{2})^{2}}\frac{v^{\prime
2}}{\widetilde{M}^{2}}\frac{\widetilde{M}_{d}M_{u}^{b}M_{u}^{c}}{v^{\prime
3}}\frac{m_{b}m_{t}^{\frac{3}{2}}\sqrt{m_{c}}}{v^{3}}\operatorname{Im}\left(V_{\rm
CKM}^{\dagger 33}V_{U}^{\prime 3b}V_{U}^{\prime\dagger b3}V_{U}^{\prime
3c}V_{U}^{\prime\dagger c2}V_{\rm CKM}^{23}\right)\tilde{f}_{duu}^{bc}\,,$
(5.3)

where a three-loop function $\tilde{f}^{bc}_{duu}$ is antisymmetric under the
permutation of $b$ and $c$, and
$\tilde{f}^{bc}_{duu}=f^{1bc}_{duu}=f^{2bc}_{duu}=f^{3bc}_{duu}$ for the
universal down-type Dirac quark mass. Here,
$V_{U}^{\prime}\equiv\overline{\Phi}(\theta_{u3},\theta_{u8})V_{U}$. A term
proportional to $m_{b}m_{t}^{2}$ (corresponding to $i=3$ and $j=k=l=3$)
vanishes by $\text{Im}(V_{\rm CKM}^{\dagger 33}V_{U}^{\prime
3b}V_{U}^{\prime\dagger b3}V_{U}^{\prime 3c}V_{U}^{\prime\dagger c3}V_{\rm
CKM}^{33})=0$. Although the above contribution ($i=3$, $j=k=3$ and $l=2$) is
suppressed by $V_{\rm CKM}^{23}\simeq 0.04$, we found that it can provide a
larger contribution than a term proportional to $m_{b}m_{t}m_{c}$ ($i=3$,
$j=l=3$ and $k=2$).

We considered the benchmark points where
$\widetilde{M}=M_{d}^{1}=M_{d}^{2}=M_{d}^{3}=M_{u}^{1}=M_{u}^{2}=10^{3}M_{u}^{3}$
and
$\widetilde{M}=M_{d}^{1}=M_{d}^{2}=M_{d}^{3}=M_{u}^{1}=M_{u}^{2}=10^{2}M_{u}^{3}$.
Both benchmark points have the degenerate down-type vector-like quark masses
and the partially degenerate up-type masses. This is because we would like to
focus on the case that the hierarchy in the SM down-type quark masses is
explained by not one in the down-type vector-like quark masses $M_{d}^{a}$ but
one in the components of the Yukawa coupling $x_{d}^{ia}$ in the seesaw
mechanism Eq. (3.7). This is motivated by the fact that the SM down-type quark
masses have a moderate hierarchy compared to the up-type ones.

We estimate the size of the three-loop function as $v^{\prime
2}/\widetilde{M}^{2}\tilde{f}_{duu}^{bc}$, where $\widetilde{M}$ is the
heaviest quark mass. It is naively expected that when a loop function is made
up of $\mathcal{O}(1)$ mass-ratio parameter, its size is maximized. In this
case a mass-ratio $m_{\varphi^{\prime}}^{2}/(M_{u}^{3})^{2}$ can become
$\mathcal{O}(1)$, so that the three-loop function would be maximized when $b$
or $c=3$. Therefore, the leading contributions in Eq. (5.3) would be dominated
by $(b,c)=(1,3)$ and $(2,3)$.

(a) $M_{u}^{1}=10^{3}M_{u}^{3}$

(b) $M_{u}^{1}=10^{2}M_{u}^{3}$

Figure 4:  The probability density functions (PDFs) of $|\delta\theta_{duu}|$
in Eq. (5.2) with only $b=1$ and $c=3$ normalized by a three-loop function
$|\bar{f}_{duu}^{13}|$ (solid line). We take
$\widetilde{M}=\widetilde{M}_{d}=M_{u}^{1}=10^{3}M_{u}^{3}$ in (a) and
$\widetilde{M}=\widetilde{M}_{d}=M_{u}^{1}=10^{2}M_{u}^{3}$ in (b). After
imposing the perturbative unitarity bound of
$\textrm{max}(x_{u})<4\pi\,(\sqrt{4\pi})$, the PDFs are expected by the dashed
(dotted) lines. The areas of dashed and dotted lines are not normalized by one
(see text). Note that the solid and dashed lines overlap in (b). If
$\bar{f}_{duu}^{13}=1$, the right area of the vertical line is excluded at 90%
CL by the neutron EDM measurement (1.1).

(a) $M_{u}^{1}=M_{u}^{2}=10^{3}M_{u}^{3}$

(b) $M_{u}^{1}=M_{u}^{2}=10^{2}M_{u}^{3}$

Figure 5:  The PDFs of $|\delta\theta_{duu}|$ in Eq. (5.2) with $b=1,\,2$ and
$c=3$ normalized by $|\bar{f}_{duu}^{13}|=|\bar{f}_{duu}^{23}|$ (solid line).
We take $\widetilde{M}=\widetilde{M}_{d}=M_{u}^{1}=M_{u}^{2}=10^{3}M_{u}^{3}$
in (a) and
$\widetilde{M}=\widetilde{M}_{d}=M_{u}^{1}=M_{u}^{2}=10^{2}M_{u}^{3}$ in (b).
After imposing the perturbative unitarity bound of
$\textrm{max}(x_{u})<4\pi\,(\sqrt{4\pi})$, the PDFs are expected by the dashed
(dotted) lines. Note that the solid and dashed lines overlap in (b). If
$\bar{f}_{duu}^{13}=\bar{f}_{duu}^{23}=1$, the right area of the vertical line
is excluded at 90% CL by the neutron EDM measurement (1.1).

In Figs. 4(a) and 4(b), we show the probability density functions (PDFs) of
the absolute value of $\delta\theta_{duu}$ in Eq. (5.2) with $b=1$ and $c=3$
(also included $b=3$ and $c=1$), by the solid lines. The total integrals of
the solid line PDFs are $1$ in all plots. When $M_{u}^{1}\gg M_{u}^{2}$
(reflecting the fact that $m_{u}\ll m_{c})$, this contribution would be
dominant in the non-vanishing $\delta\theta_{duu}$, see Eq. (5.3). In the
PDFs, the $\delta\theta_{duu}$ is normalized by $|\bar{f}_{duu}^{13}|$ where
we define
$\bar{f}_{duu}^{13}=\tilde{f}^{13}_{duu}-\tilde{f}^{31}_{duu}=2\tilde{f}^{13}_{duu}$.
In Fig. 4(a) and 4(b), we take
$\widetilde{M}=\widetilde{M}_{d}=M_{u}^{1}=10^{3}M_{u}^{3}$ and
$\widetilde{M}=\widetilde{M}_{d}=M_{u}^{1}=10^{2}M_{u}^{3}$, respectively,
with $v^{\prime}=M_{u}^{3}$, by solid lines, whose hierarchy is concerned with
the top quark masses in the SM and the seesaw mechanism in Eq. (3.7). After
fixing the ratios of the mass parameters, and assuming that all light quark
masses and the CKM components reproduce the SM one, the remaining free
parameters are the three rotation angles and one $CP$-violating phase in the
$V_{U}$ matrix and two additional $CP$-violating phases in
$\overline{\Phi}(\theta_{u3},\theta_{u8})$ in Eq. (3.10). In order to obtain
the PDFs, we varied these six angles from $0$ to $2\pi$ with equal
probability. This provides the solid line PDFs in Fig. 4(a) and 4(b).

Although the radiative $\bar{\theta}$ parameter should be renormalization
scale invariant,#14#14#14Recently, there is a discussion of a perturbative
running of a renormalized $\theta$ parameter within renormalizable theories
[59]. In our setup, the bare $\theta$ parameter (hence the renormalized one)
is forbidden. in order to investigate the perturbative unitarity bounds for
the Yukawa interactions, we used the SM (running) quark masses at
$\mu=1\,$TeV.#15#15#15We take $m_{t}=143\,{\rm GeV},m_{b}=2.41\,{\rm
GeV},m_{c}=528\,{\rm MeV}$, and $m_{s}=45.4\,{\rm MeV}$ at $\mu=1$ TeV [60].
Furthermore, by assuming $\bar{f}_{duu}^{13}=1$, the vertical line in the
figures stands for the experimental upper bound from the neutron EDM
measurement in Eq. (1.1) and the right area is excluded at 90% CL.

Given the fixed values of $\overline{\Phi}(\theta_{u3},\theta_{u8})V_{U}$, one
can obtain the eigenvalues of the up-type Yukawa matrix $x_{u}$. In the PDF
analysis, we find that the eigenvalues of the Yukawa matrix $x_{u}$ can be
larger than $O(1)$ easily depending on
$\overline{\Phi}(\theta_{u3},\theta_{u8})V_{U}$. Therefore, we impose the
maximal eigenvalue to be smaller than $4\pi$ or $\sqrt{4\pi}$ for the dashed
or dotted lines, respectively, as the perturbative unitarity bound. When the
maximal eigenvalue exceeds the unitary bound, we discarded these points in the
PDFs. In order to display the reduction in statistics as a result of setting
the perturbative bound, we do not normalize the dashed and dotted line PDFs by
$1$, and hence their total integration is less than $1$.

---
| $M_{u}^{3}/M_{u}^{1}=10^{-3}$ | $M_{u}^{3}/M_{u}^{1}=10^{-2}$
$\textrm{max}(x_{u})<4\pi$ | $42.1\%\,(58.9\%)$ | $32.2\%\,(46.2\%)$
$\textrm{max}(x_{u})<\sqrt{4\pi}$ | $3.65\%\,(10.6\%)$ | $38.8\%\,(55.9\%)$
Table 2:  Ratios of excluded parameter regions for the parameter sets in Fig.
4 from the current neutron EDM measurement, under an assumption of
$\bar{f}_{duu}^{13}=1$. The numbers in parentheses are those in Fig. 5 under
an assumption of $\bar{f}_{duu}^{13}=\bar{f}_{duu}^{23}=1$. The parameter sets
are restricted by the perturbative bound of the Yukawa coupling $x_{u}$ as
$4\pi$ or $\sqrt{4\pi}$. Here,
$\widetilde{M}=\widetilde{M}_{d}=M_{u}^{1}=M_{u}^{2}$.

The ratios of the excluded parameter regions by the neutron EDM measurement in
each of the PDFs are shown in Table 2 under an assumption of
$\bar{f}_{duu}^{13}=1$. Here, the ratios of the excluded parameter regions are
obtained by comparing the areas of the dashed or dotted line PDFs; right area
of the vertical line over the total area where the unitarity bound is imposed.
We found that some fractions of parameter regions have already been excluded
even if the perturbative bound is imposed. In the case of the perturbative
bound of $4\pi$ to the eigenvalues of $x_{u}$ matrix, about 30–40% of the
whole parameter region is already excluded. On the other hand, in the case of
$\sqrt{4\pi}$, the dependence of $M_{u}^{3}$ appears obviously. In particular,
only 3.65% is excluded in the case of $M_{u}^{1}=10^{3}M_{u}^{3}$. These
results show that this model is sensitive to the current bound from the
neutron EDM experiments and has a possibility to be explored by future
improvement of the experiments.

Moreover, in Figs. 5(a) and 5(b), we show the PDFs of the absolute value of
$\delta\theta_{duu}$ with $(b,c)=(1,3)$ plus $(2,3)$ (also included $(3,1)$
and $(3,2)$), with assuming $M_{u}^{1}=M_{u}^{2}$ which should be a somewhat
aggressive parameter choice in light of $m_{u}\ll m_{c}$. We observed that the
PDFs in Fig. 5 are slightly larger than the PDFs in Fig. 4. The ratios of
excluded parameter regions by the neutron EDM measurement are shown in Table 2
at the numbers in parentheses.

Above estimation of the leading contribution to the radiative $\bar{\theta}$
parameter is numerically consistent with the latest analysis Eq. (28) of Ref.
[49], where the ordinary calculation method is used. In this paper, we showed
for the first time that an effect from the perturbative unitarity bound is
important and it reduces the radiative corrections to $\bar{\theta}$.

### 5.2 Case for the GIM by universal vector-like mass

Next, we investigate a case that all Dirac quark masses are degenerate as
$\widetilde{M}=M_{u}^{a}=M_{d}^{a}$. In this case, the GIM-like mechanism
occurs and then the CKM matrix becomes a unique source of the $CP$-violating
phase. As a consequence, the radiative $\bar{\theta}$ parameter would be
significantly suppressed, and it is expected as the minimum value of the
$\bar{\theta}$ parameter in the minimal LR symmetric model. The induced
$\bar{\theta}$ parameter should be proportional to the Jarlskog invariant of
the CKM matrix $J_{\rm CKM}$, which is given by $\text{Im}(V_{\rm
CKM}^{ij}V_{\rm CKM}^{{\dagger}jk}V_{\rm CKM}^{kl}V_{\rm
CKM}^{{\dagger}li})=J_{\rm CKM}\sum_{m,n}\epsilon_{ikm}\epsilon_{jln}$ [61]
and $J_{\rm CKM}=(3.08{}^{+0.15}_{-0.13})\times 10^{-5}$ [15].

If the $\bar{\theta}$ parameter is induced at three-loop level, the
contribution would be given as

$\displaystyle\delta\theta_{\rm minimum}$
$\displaystyle\approx\frac{1}{\left(16\pi^{2}\right)^{2}}\frac{v^{\prime
8}}{\widetilde{M}^{8}}\operatorname{Im}\operatorname{Tr}\left(A_{u}^{2}A_{d}^{2}A_{u}A_{d}\right){f_{\rm
CKM}}$ $\displaystyle=\frac{1}{\left(16\pi^{2}\right)^{2}}\frac{v^{\prime
2}}{\widetilde{M}^{2}}J_{\rm CKM}{f_{\rm CKM}}$
$\displaystyle\quad\times\frac{1}{v^{6}}\left(m_{t}-m_{c}\right)\left(m_{t}-m_{u}\right)\left(m_{c}-m_{u}\right)\left(m_{b}-m_{s}\right)\left(m_{b}-m_{d}\right)\left(m_{s}-m_{d}\right),$
$\displaystyle\simeq 1\times 10^{-{19}}\frac{v^{\prime
2}}{\widetilde{M}^{2}}{f_{\rm CKM}}\,,$ (5.7)

where $f_{\rm CKM}$ is a dimensionless loop function of $O(1)$. It is assumed
that in the three-loop diagrams, two scalar bosons are exchanged inside the
fermion loop, and the other scalar lines are replaced by $v^{\prime}$. The
above contribution is proportional to $v^{\prime 8}$. The perturbativity of
the top Yukawa requires $\widetilde{M}\simeq{v^{\prime}}$. Then,
$\delta\theta_{\rm minimum}$ is expected as at most $10^{-{19}}$ or smaller
when all vector-like quark masses are degenerate. This size is comparable to
or smaller than the CKM phase contribution to the $\bar{\theta}$ parameter in
the SM at four-loop level, evaluated in Ref. [39] (though the top quark had
been assumed to be lighter than the $W$ boson). This is because the quark mass
suppression of the contribution in the SM is milder than Eq. (5.7). In
addition, the neutron EDM induced by the CKM phase via long-distance hadronic
contributions in the SM [62, 63] is much larger than by the contribution to
the $\bar{\theta}$ parameter in the above benchmark point. Then, a more
precise evaluation of the $\bar{\theta}$ parameter in the benchmark point is
too academic and beyond our scope.

## 6 Conclusions and discussion

When the QCD axion is absent, one has to solve the strong $CP$ problem by an
additional discrete symmetry. The extended parity with the LR gauge symmetry
can solve it with generating the SM as the low-energy theory. However, it is
known that the radiative $\bar{\theta}$ parameter is induced from the soft
symmetry breaking of the parity. In this paper, we first formulated a novel
method of direct loop-diagrammatic calculation of the radiative $\bar{\theta}$
parameter by using Fock-Schwinger gauge. This approach should be more robust
than the the ordinary calculation method based on the chiral rotations. By
using the Fock-Schwinger gauge method, we confirmed a seminal result that two-
loop level $\bar{\theta}$ vanishes completely in the minimal LR symmetric
model.

Furthermore, we estimated the size of the leading contributions to the non-
vanishing radiative $\bar{\theta}$ parameter at three-loop level. We derive
the parameterization by the physical parameters based on the seesaw mechanism
in the LR symmetric model, and we obtained the probability density functions
of the radiative $\bar{\theta}$ parameter by varying all free but physical
parameters. Here, we also investigated the impact of the perturbative
unitarity conditions for the LR symmetric Yukawa matrices. It is found that
the resultant $\bar{\theta}$ parameters are partially excluded by the current
neutron EDM bound. It implies that this model has a possibility to be explored
by future improvement of the experiments. One should note that the minimal LR
symmetric model predicts that all hadronic EDMs are dominated by the radiative
$\bar{\theta}$ parameter. Therefore, this model can predict a distinctive and
non-vanishing correlation between the neutron, proton, nuclei (2H, 3He),
diamagnetic atoms (Hg, Ra), paramagnetic atoms and molecules (YbF, HfF, ThO)
EDMs [49, 64, 65, 66].

The large mass-scale difference between $v$ and $v^{\prime}$ can be explained
by a Higgs parity mechanism that predicts $\lambda_{\rm
SM}(\mu=v^{\prime})\simeq 0$ and $v^{\prime}=\mathcal{O}(10^{10})$ GeV [42,
67, 68, 69]. Since the above estimation of the radiative $\bar{\theta}$
parameter is insensitive to the mass scale itself of the right-handed sector,
one could predict the neutron EDM for the case of
$v^{\prime}=\mathcal{O}(10^{10})$ GeV.

We also comment on the radiative $\bar{\theta}$ parameter in the spontaneous
$CP$ violation (Nelson-Barr) model [32, 33, 34], which is another model that
can explain the strong $CP$ problem by the discrete symmetry. In Ref. [70],
the radiative $\bar{\theta}$ parameter has been investigated in detail. It is
found that although the reducible $\bar{\theta}$, which comes from quartic
couplings of scalar fields responsible for the spontaneous $CP$ violation with
the SM Higgs one, is induced at two-loop level, it can be numerically
neglected if the couplings are small enough. On the other hand, the
irreducible $\bar{\theta}$, which is related to the CKM phase, is induced at
three-loop level and is safely below the current experimental bounds. The
loop-diagrammatic approach we proposed would provide a more robust estimation
if possible.

It would be an interesting direction to investigate correlations between the
radiative $\bar{\theta}$ parameter and other (flavor) observables in the
minimal LR symmetric model: The lepton flavor universality violation in $B\to
D^{(\ast)}l\nu$ ($R(D^{(\ast)})$ anomaly) (the recent review [71, 72]), the
Cabibbo angle anomaly (the recent review [73]) and the $W$-boson mass anomaly
[74] could be explained [75, 76, 77]. Furthermore, the investigation of a
correlation with electroweak-like baryogenesis in the right-handed sector
should be an attractive prospect [78, 79].

###### Acknowledgements.

We would like to acknowledge Keisuke Harigaya for a collaboration at the early
stage. This work is supported by the JSPS Grant-in-Aid for Scientific Research
Grant No. 20H01895 (J.H.), No. 21K03572 (J.H.) and for Early-Career Scientists
Grant No. 19K14706 (T.K.). The work of J.H. is also supported by World Premier
International Research Center Initiative (WPI Initiative), MEXT, Japan. This
work is also supported by JSPS Core-to-Core Program Grant No. JPJSCCA20200002.
This work was financially supported by JST SPRING, Grant Number JPMJSP2125.
The author (A.Y.) would like to take this opportunity to thank the
“Interdisciplinary Frontier Next-Generation Researcher Program of the Tokai
Higher Education and Research System.”

##  Appendix A Loop functions

We define the loop functions in this section. The two-loop functions used in
this paper are given as

$\displaystyle
I_{(n_{1},\cdots;m_{1},\cdots)}(x_{1},\cdots;x_{2},\cdots;x_{3})$
$\displaystyle\quad=(16\pi^{2}\mu^{2\epsilon})^{2}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{d^{d}q}{(2\pi)^{d}}\frac{1}{[p^{2}-x_{1}]^{n_{1}}\cdots[q^{2}-x_{2}]^{m_{1}}\cdots[(p+q)^{2}-x_{3}]}\,,$
(A.1)

where the dimensional regularization is used on $d=4-2\epsilon$ dimension and
$\mu$ is the renormalization scale. The functions can also be derived by
derivative or finite difference of $I(x_{1};x_{2};x_{3})$ ($\equiv
I_{(1;1)}(x_{1};x_{2};x_{3})$) as

$\displaystyle I_{(n;m)}(x_{1};x_{2};x_{3})$ $\displaystyle=$
$\displaystyle\frac{1}{(n-1)!(m-1)!}\frac{d^{n-1}}{dx_{1}^{n-1}}\frac{d^{m-1}}{dx_{2}^{m-1}}I(x_{1};x_{2};x_{3})\,,$
(A.2) $\displaystyle I_{(1,1;1)}(x_{1},x_{1}^{\prime};x_{2};x_{3})$
$\displaystyle=$
$\displaystyle\frac{1}{x_{1}-x_{1}^{\prime}}\left[I(x_{1};x_{2};x_{3})-I(x_{1}^{\prime};x_{2};x_{3})\right]\,,$
(A.3)

where

$\displaystyle I(x_{1};x_{2};x_{3})$ $\displaystyle\equiv$ $\displaystyle
I_{(1;1)}(x_{1};x_{2};x_{3})$ (A.4) $\displaystyle=$
$\displaystyle(16\pi^{2}\mu^{2\epsilon})^{2}\int\frac{d^{d}pd^{d}q}{(2\pi)^{2d}}\frac{1}{[p^{2}-x_{1}][q^{2}-x_{2}][(p+q)^{2}-x_{3}]}\,.$

The explicit form of $I(x_{1};x_{2};x_{3})$ is given as

$\displaystyle I(x_{1};x_{2};x_{3})$ $\displaystyle=$
$\displaystyle\bar{I}_{\epsilon}(x_{1};x_{2};x_{3})+\bar{I}(x_{1};x_{2};x_{3})\,,$
(A.5)

with the UV divergent part

$\displaystyle\bar{I}_{\epsilon}(x_{1};x_{2};x_{3})=-\sum_{i=1,2,3}x_{i}\left[\frac{1}{2\epsilon^{2}}-\frac{1}{\epsilon}\left(\log\frac{x_{i}}{Q^{2}}-\frac{3}{2}\right)+\left(1+\frac{\pi^{2}}{12}-\log\frac{x_{i}}{Q^{2}}+\frac{1}{2}\log^{2}\frac{x_{i}}{Q^{2}}\right)\right]\,,$
(A.6)

while the finite part

$\displaystyle\bar{I}(x_{1};x_{2};x_{3})=-\frac{1}{2}\biggl{[}(-x_{1}+x_{2}+x_{3})\log\frac{x_{2}}{Q^{2}}\log\frac{x_{3}}{Q^{2}}+(x_{1}-x_{2}+x_{3})\log\frac{x_{1}}{Q^{2}}\log\frac{x_{3}}{Q^{2}}$
$\displaystyle\quad+(x_{1}+x_{2}-x_{3})\log\frac{x_{1}}{Q^{2}}\log\frac{x_{2}}{Q^{2}}-4\left(x_{1}\log\frac{x_{1}}{Q^{2}}+x_{2}\log\frac{x_{2}}{Q^{2}}+x_{3}\log\frac{x_{3}}{Q^{2}}\right)$
$\displaystyle\quad+5(x_{1}+x_{2}+x_{3})+\xi(x_{1},x_{2},x_{3})\biggr{]}\,,$
(A.7)
$\displaystyle\xi(x_{1},x_{2},x_{3})=R\left[2\log\left(\frac{x_{3}+x_{1}-x_{2}-R}{2x_{3}}\right)\log\left(\frac{x_{3}-x_{1}+x_{2}-R}{2x_{3}}\right)-\log\frac{x_{1}}{x_{3}}\log\frac{x_{2}}{x_{3}}\right.$
$\displaystyle\quad\left.-2{\rm
Li}_{2}\left(\frac{x_{3}+x_{1}-x_{2}-R}{2x_{3}}\right)-2{\rm
Li}_{2}\left(\frac{x_{3}-x_{1}+x_{2}-R}{2x_{3}}\right)+\frac{\pi^{2}}{3}\right]\,,$
(A.8)

where
$R=\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{3}-2x_{1}x_{2}-2x_{2}x_{3}-2x_{3}x_{1}}$
and $Q^{2}\equiv 4\pi\mu^{2}{\rm e}^{-\gamma_{E}}$ [80, 81, 82]. Note that
although the last terms of $\bar{I}_{\epsilon}$ are UV finite, they do not
affect any physical quantity [83]. These terms are suppressed by
$\mathcal{O}(\epsilon)$ at one-loop level, while they are uplifted to
$\mathcal{O}(\epsilon^{0})$ at two-loop level.

### A.1 $I_{(1;3)}$ and $I_{(3;1)}$

The two-loop function $I_{(3;1)}$ is also given by the two-point one-loop
function $B_{0}$ as

$\displaystyle I_{(3;1)}(x_{1};x_{2};x_{3})$ $\displaystyle=$
$\displaystyle-(16\pi^{2}\mu^{2\epsilon})\int\frac{d^{d}p}{i(2\pi)^{d}}\frac{1}{(p^{2}-x_{1})^{3}}B_{0}(p^{2},x_{2},x_{3})\,,$
(A.9)

where

$\displaystyle B_{0}(p^{2},x_{2},x_{3})$
$\displaystyle\equiv(16\pi^{2}\mu^{2\epsilon})\int\frac{d^{d}q}{i(2\pi)^{d}}\frac{1}{[q^{2}-x_{2}][(p+q)^{2}-x_{3}]}$
$\displaystyle=\frac{1}{\epsilon}-F_{0}(p^{2},x_{2},x_{3})\,,$ (A.10)
$\displaystyle F_{0}(p^{2},x_{2},x_{3})$
$\displaystyle=\int^{1}_{0}dz\log\frac{-z(1-z)p^{2}+zx_{2}+(1-z)x_{3}}{Q^{2}}\,.$
(A.11)

Similarly, $\bar{I}_{(1;3)}$ is obtained by replacements of
$x_{1}\leftrightarrow x_{2}$ in the above equations of $\bar{I}_{(3;1)}$.

The finite part of ${I}_{(3;1)}$ ($=\bar{I}_{(3;1)}$) is given as

$\displaystyle\bar{I}_{(3;1)}(x_{1};x_{2};x_{3})$ $\displaystyle=$
$\displaystyle(16\pi^{2}\mu^{2\epsilon})\int\frac{d^{d}p}{i(2\pi)^{d}}\frac{1}{(p^{2}-x_{1})^{3}}F_{0}(p^{2},x_{2},x_{3})\,,$
(A.12)

while the divergent part is

$\displaystyle\bar{I}_{\epsilon{(3;1)}}(x_{1};x_{2};x_{3})=\frac{1}{2x_{1}}\left(\frac{1}{\epsilon}-\log\frac{x_{1}}{Q^{2}}\right)\,.$
(A.13)

In the case of $x_{1}\ll x_{3}$, $\bar{I}_{(3;1)}$ is given as

$\displaystyle\bar{I}_{(3;1)}(x_{1};x_{2};x_{3})$ $\displaystyle\simeq$
$\displaystyle-\frac{1}{2x_{1}}F_{0}(0,x_{2},x_{3})\,,$ (A.14)

where

$\displaystyle F_{0}(0,x_{2},x_{3})$ $\displaystyle=$
$\displaystyle\int^{1}_{0}dz\log\frac{zx_{2}+(1-z)x_{3}}{Q^{2}}$ (A.15)
$\displaystyle=$
$\displaystyle-\frac{x_{2}-x_{3}-x_{2}\log(x_{2}/Q^{2})+x_{3}\log(x_{3}/Q^{2})}{x_{2}-x_{3}}\,.$

On the other hand, in the case of $x_{2}\ll x_{3}$, $\bar{I}_{(3;1)}$ is given
as

$\displaystyle\bar{I}_{(3;1)}(x_{1};x_{2};x_{3})$ $\displaystyle\simeq$
$\displaystyle\frac{1}{2x_{1}}\frac{x_{1}-x_{3}-x_{1}\log(x_{1}/Q^{2})+x_{3}\log(x_{3}/Q^{2})}{x_{1}-x_{3}}\,.$
(A.16)

Furthermore, we define the loop function at one-loop level as

$\displaystyle B_{(n_{1},\cdots)}(p^{2},x_{1},\cdots;x_{3})$ $\displaystyle=$
$\displaystyle(16\pi^{2}\mu^{2\epsilon})\int\frac{d^{d}q}{i(2\pi)^{d}}\frac{1}{[q^{2}-x_{1}]^{n_{1}}\cdots[(p+q)^{2}-x_{3}]}\,.$
(A.17)

### A.2 $I_{(2;2)}$

The two-loop function $I_{(2;2)}$ does not contain the UV divergence.
$I_{(2;2)}(x_{1};x_{2};x_{3})$ is a symmetric function in variables $x_{1}$
and $x_{2}$. In the case of $x_{1},\,x_{2}\ll x_{3}$, $I_{(2;2)}$ is given as

$\displaystyle I_{(2;2)}(x_{1};x_{2};x_{3})$ $\displaystyle\simeq$
$\displaystyle\frac{1}{x_{3}}\left(\frac{\pi^{2}}{3}+\log\frac{x_{1}x_{2}}{x_{3}^{2}}+\log\frac{x_{1}}{x_{3}}\log\frac{x_{2}}{x_{3}}\right)\,.$
(A.18)

In the case of $x_{1}\ll x_{2}<x_{3}$, we find

$\displaystyle I_{(2;2)}(x_{1};x_{2};x_{3})$ $\displaystyle\simeq$
$\displaystyle\frac{x_{3}}{(x_{3}-x_{2})^{2}}\left\\{\frac{\pi^{2}}{3}-\frac{x_{2}}{x_{3}}\log\frac{x_{1}}{x_{2}}+\log\frac{x_{1}x_{2}}{x_{3}^{2}}+\log\frac{x_{1}}{x_{3}}\log\frac{x_{2}}{x_{3}}\right.$
(A.19)
$\displaystyle\qquad\left.-2\left[\log\frac{x_{2}}{x_{3}}\log\left(1-\frac{x_{2}}{x_{3}}\right)+{\rm
Li}_{2}\left(\frac{x_{2}}{x_{3}}\right)\right]\right\\}\,,$

while for $x_{1}\ll x_{3}<x_{2}$

$\displaystyle
I_{(2;2)}(x_{1};x_{2};x_{3})\simeq\frac{x_{3}}{(x_{3}-x_{2})^{2}}\left[-\frac{x_{2}}{x_{3}}\log\frac{x_{1}}{x_{2}}+\log\frac{x_{1}x_{2}}{x_{3}^{2}}+\log\frac{x_{1}}{x_{3}}\log\frac{x_{2}}{x_{3}}+2{\rm
Li}_{2}\left(1-\frac{x_{2}}{x_{3}}\right)\right]\,.$ (A.20)

## References

  * [1] A. A. Belavin, A. M. Polyakov, A. S. Schwartz, and Y. S. Tyupkin, “Pseudoparticle Solutions of the Yang-Mills Equations,” Phys. Lett. B 59 (1975) 85–87.
  * [2] G. ’t Hooft, “Symmetry Breaking Through Bell-Jackiw Anomalies,” Phys. Rev. Lett. 37 (1976) 8–11.
  * [3] C. G. Callan, Jr., R. F. Dashen, and D. J. Gross, “The Structure of the Gauge Theory Vacuum,” Phys. Lett. B 63 (1976) 334–340.
  * [4] R. Jackiw and C. Rebbi, “Vacuum Periodicity in a Yang-Mills Quantum Theory,” Phys. Rev. Lett. 37 (1976) 172–175.
  * [5] E. Witten, “Dyons of Charge e theta/2 pi,” Phys. Lett. B 86 (1979) 283–287.
  * [6] V. Baluni, “CP Violating Effects in QCD,” Phys. Rev. D 19 (1979) 2227–2230.
  * [7] R. J. Crewther, P. Di Vecchia, G. Veneziano, and E. Witten, “Chiral Estimate of the Electric Dipole Moment of the Neutron in Quantum Chromodynamics,” Phys. Lett. B 88 (1979) 123. [Erratum: Phys.Lett.B 91, 487 (1980)].
  * [8] C. Abel et al., “Measurement of the Permanent Electric Dipole Moment of the Neutron,” Phys. Rev. Lett. 124 (2020) 081803 [arXiv:2001.11966].
  * [9] J. Liang, et al., “Nucleon Electric Dipole Moment from the $\theta$ Term with Lattice Chiral Fermions.” arXiv:2301.04331.
  * [10] J. Dragos, T. Luu, A. Shindler, J. de Vries, and A. Yousif, “Confirming the Existence of the strong CP Problem in Lattice QCD with the Gradient Flow,” Phys. Rev. C 103 (2021) 015202 [arXiv:1902.03254].
  * [11] N. Cabibbo, “Unitary Symmetry and Leptonic Decays,” Phys. Rev. Lett. 10 (1963) 531–533.
  * [12] M. Kobayashi and T. Maskawa, “CP Violation in the Renormalizable Theory of Weak Interaction,” Prog. Theor. Phys. 49 (1973) 652–657.
  * [13] CKMfitter Group Collaboration, “CP violation and the CKM matrix: Assessing the impact of the asymmetric $B$ factories,” Eur. Phys. J. C 41 (2005) 1–131 [hep-ph/0406184]. updated results and plots available at: http://ckmfitter.in2p3.fr.
  * [14] L.-L. Chau and W.-Y. Keung, “Comments on the Parametrization of the Kobayashi-Maskawa Matrix,” Phys. Rev. Lett. 53 (1984) 1802.
  * [15] Particle Data Group Collaboration, “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * [16] H. Georgi and I. N. McArthur, “INSTANTONS AND THE mu QUARK MASS.”.
  * [17] D. B. Kaplan and A. V. Manohar, “Current Mass Ratios of the Light Quarks,” Phys. Rev. Lett. 56 (1986) 2004.
  * [18] K. Choi, C. W. Kim, and W. K. Sze, “Mass Renormalization by Instantons and the Strong CP Problem,” Phys. Rev. Lett. 61 (1988) 794.
  * [19] T. Banks, Y. Nir, and N. Seiberg, “Missing (up) mass, accidental anomalous symmetries, and the strong CP problem,” in 2nd IFT Workshop on Yukawa Couplings and the Origins of Mass, pp. 26–41.  1994\.  hep-ph/9403203.
  * [20] M. Srednicki, “Comment on ”Ambiguities in the up-quark mass”,” Phys. Rev. Lett. 95 (2005) 059101 [hep-ph/0503051].
  * [21] C. Alexandrou, et al., “Ruling Out the Massless Up-Quark Solution to the Strong $\boldsymbol{CP}$ Problem by Computing the Topological Mass Contribution with Lattice QCD,” Phys. Rev. Lett. 125 (2020) 232001 [arXiv:2002.07802].
  * [22] R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38 (1977) 1440–1443.
  * [23] S. Weinberg, “A New Light Boson?” Phys. Rev. Lett. 40 (1978) 223–226.
  * [24] F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978) 279–282.
  * [25] M. Kamionkowski and J. March-Russell, “Planck scale physics and the Peccei-Quinn mechanism,” Phys. Lett. B 282 (1992) 137–141 [hep-th/9202003].
  * [26] R. Holman, et al., “Solutions to the strong CP problem in a world with gravity,” Phys. Lett. B 282 (1992) 132–136 [hep-ph/9203206].
  * [27] S. M. Barr and D. Seckel, “Planck scale corrections to axion models,” Phys. Rev. D 46 (1992) 539–549.
  * [28] M. A. B. Beg and H. S. Tsao, “Strong P, T Noninvariances in a Superweak Theory,” Phys. Rev. Lett. 41 (1978) 278.
  * [29] R. N. Mohapatra and G. Senjanovic, “Natural Suppression of Strong p and t Noninvariance,” Phys. Lett. B 79 (1978) 283–286.
  * [30] K. S. Babu and R. N. Mohapatra, “A Solution to the Strong CP Problem Without an Axion,” Phys. Rev. D 41 (1990) 1286.
  * [31] S. M. Barr, D. Chang, and G. Senjanovic, “Strong CP problem and parity,” Phys. Rev. Lett. 67 (1991) 2765–2768.
  * [32] A. E. Nelson, “Naturally Weak CP Violation,” Phys. Lett. B 136 (1984) 387–391.
  * [33] S. M. Barr, “Solving the Strong CP Problem Without the Peccei-Quinn Symmetry,” Phys. Rev. Lett. 53 (1984) 329.
  * [34] S. M. Barr, “A Natural Class of Nonpeccei-quinn Models,” Phys. Rev. D 30 (1984) 1805.
  * [35] S. L. Adler, “Axial vector vertex in spinor electrodynamics,” Phys. Rev. 177 (1969) 2426–2438.
  * [36] J. S. Bell and R. Jackiw, “A PCAC puzzle: $\pi^{0}\to\gamma\gamma$ in the $\sigma$ model,” Nuovo Cim. A 60 (1969) 47–61.
  * [37] K. Fujikawa, “Path Integral Measure for Gauge Invariant Fermion Theories,” Phys. Rev. Lett. 42 (1979) 1195–1198.
  * [38] J. R. Ellis and M. K. Gaillard, “Strong and Weak CP Violation,” Nucl. Phys. B 150 (1979) 141–162.
  * [39] I. B. Khriplovich, “Quark Electric Dipole Moment and Induced $\theta$ Term in the Kobayashi-Maskawa Model,” Phys. Lett. B 173 (1986) 193–196.
  * [40] J.-M. Gérard and P. Mertens, “Weakly-induced strong CP-violation,” Phys. Lett. B 716 (2012) 316–321 [arXiv:1206.0914].
  * [41] A. Czarnecki and B. Krause, “Neutron Electric Dipole Moment in the Standard Model: Complete Three-Loop Calculation of the Valence Quark Contributions,” Phys. Rev. Lett. 78 (1997) 4339–4342 [hep-ph/9704355].
  * [42] L. J. Hall and K. Harigaya, “Implications of Higgs Discovery for the Strong CP Problem and Unification,” JHEP 10 (2018) 130 [arXiv:1803.08119].
  * [43] N. Craig, I. Garcia Garcia, G. Koszegi, and A. McCune, “P not PQ,” JHEP 09 (2021) 130 [arXiv:2012.13416].
  * [44] S. L. Adler and W. A. Bardeen, “Absence of higher order corrections in the anomalous axial vector divergence equation,” Phys. Rev. 182 (1969) 1517–1536.
  * [45] G. ’t Hooft and M. J. G. Veltman, “Regularization and Renormalization of Gauge Fields,” Nucl. Phys. B 44 (1972) 189–213.
  * [46] H. Georgi, T. Tomaras, and A. Pais, “$\theta_{\rm QCD}$ without instantons,” Phys. Rev. D 23 (1981) 469–472.
  * [47] V. A. Novikov, M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, “Calculations in External Fields in Quantum Chromodynamics. Technical Review,” Fortsch. Phys. 32 (1984) 585.
  * [48] T. Abe, J. Hisano, and R. Nagai, “Model independent evaluation of the Wilson coefficient of the Weinberg operator in QCD,” JHEP 03 (2018) 175 [arXiv:1712.09503]. [Erratum: JHEP 09, 020 (2018)].
  * [49] J. de Vries, P. Draper, and H. H. Patel, “Do Minimal Parity Solutions to the Strong $CP$ Problem Work?” arXiv:2109.01630.
  * [50] A. Davidson and K. C. Wali, “$SU(5)_{L}\otimes SU(5)_{R}$ Hybrid Unification,” Phys. Rev. Lett. 58 (1987) 2623.
  * [51] A. Davidson and K. C. Wali, “Universal Seesaw Mechanism?” Phys. Rev. Lett. 59 (1987) 393.
  * [52] C.-H. Lee and R. N. Mohapatra, “Vector-Like Quarks and Leptons, SU(5) $\otimes$ SU(5) Grand Unification, and Proton Decay,” JHEP 02 (2017) 080 [arXiv:1611.05478].
  * [53] K. S. Babu and X. G. He, “Dirac Neutrino Masses as Two-Loop Radiative Corrections,” Mod. Phys. Lett. A 4 (1989) 61.
  * [54] K. S. Babu, X.-G. He, M. Su, and A. Thapa, “Naturally light Dirac and pseudo-Dirac neutrinos from left-right symmetry,” JHEP 08 (2022) 140 [arXiv:2205.09127].
  * [55] J. A. Dror, D. Dunsky, L. J. Hall, and K. Harigaya, “Sterile Neutrino Dark Matter in Left-Right Theories,” JHEP 07 (2020) 168 [arXiv:2004.09511].
  * [56] D. Dunsky, L. J. Hall, and K. Harigaya, “Sterile Neutrino Dark Matter and Leptogenesis in Left-Right Higgs Parity,” JHEP 01 (2021) 125 [arXiv:2007.12711].
  * [57] J. A. Casas and A. Ibarra, “Oscillating neutrinos and $\mu\to e,\gamma$,” Nucl. Phys. B 618 (2001) 171–204 [hep-ph/0103065].
  * [58] J. R. Ellis, J. Hisano, M. Raidal, and Y. Shimizu, “A New parametrization of the seesaw mechanism and applications in supersymmetric models,” Phys. Rev. D 66 (2002) 115013 [hep-ph/0206110].
  * [59] A. Valenti and L. Vecchi, “Perturbative running of the topological angles,” JHEP 01 (2023) 131 [arXiv:2210.09328].
  * [60] K. G. Chetyrkin, J. H. Kuhn, and M. Steinhauser, “RunDec: A Mathematica package for running and decoupling of the strong coupling and quark masses,” Comput. Phys. Commun. 133 (2000) 43–65 [hep-ph/0004189].
  * [61] C. Jarlskog, “Commutator of the Quark Mass Matrices in the Standard Electroweak Model and a Measure of Maximal $CP$ Nonconservation,” Phys. Rev. Lett. 55 (1985) 1039.
  * [62] S. Dar, “The Neutron EDM in the SM: A Review.” hep-ph/0008248.
  * [63] T. Mannel and N. Uraltsev, “Loop-Less Electric Dipole Moment of the Nucleon in the Standard Model,” Phys. Rev. D 85 (2012) 096002 [arXiv:1202.6270].
  * [64] W. Dekens, et al., “Unraveling models of CP violation through electric dipole moments of light nuclei,” JHEP 07 (2014) 069 [arXiv:1404.6082].
  * [65] J. de Vries, P. Draper, K. Fuyuto, J. Kozaczuk, and D. Sutherland, “Indirect Signs of the Peccei-Quinn Mechanism,” Phys. Rev. D 99 (2019) 015042 [arXiv:1809.10143].
  * [66] J. de Vries, P. Draper, K. Fuyuto, J. Kozaczuk, and B. Lillard, “Uncovering an axion mechanism with the EDM portfolio,” Phys. Rev. D 104 (2021) 055039 [arXiv:2107.04046].
  * [67] D. Dunsky, L. J. Hall, and K. Harigaya, “Higgs Parity, Strong CP, and Dark Matter,” JHEP 07 (2019) 016 [arXiv:1902.07726].
  * [68] L. J. Hall and K. Harigaya, “Higgs Parity Grand Unification,” JHEP 11 (2019) 033 [arXiv:1905.12722].
  * [69] D. Dunsky, L. J. Hall, and K. Harigaya, “Dark Matter, Dark Radiation and Gravitational Waves from Mirror Higgs Parity,” JHEP 02 (2020) 078 [arXiv:1908.02756].
  * [70] A. Valenti and L. Vecchi, “The CKM phase and $\overline{\theta}$ in Nelson-Barr models,” JHEP 07 (2021) 203 [arXiv:2105.09122].
  * [71] HFLAV Collaboration, “Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of 2021.” arXiv:2206.07501.
  * [72] S. Iguro, T. Kitahara, and R. Watanabe, “Global fit to $b\to c\tau\nu$ anomalies 2022 mid-autumn.” arXiv:2210.10751.
  * [73] A. Crivellin, M. Kirk, T. Kitahara, and F. Mescia, “Global Fit of Modified Quark Couplings to EW Gauge Bosons and Vector-Like Quarks in Light of the Cabibbo Angle Anomaly.” arXiv:2212.06862.
  * [74] CDF Collaboration, “High-precision measurement of the $W$ boson mass with the CDF II detector,” Science 376 (2022) 170–176.
  * [75] K. S. Babu, B. Dutta, and R. N. Mohapatra, “A theory of R(D∗, D) anomaly with right-handed currents,” JHEP 01 (2019) 168 [arXiv:1811.04496].
  * [76] K. S. Babu and R. Dcruz, “Resolving $W$ Boson Mass Shift and CKM Unitarity Violation in Left-Right Symmetric Models with Universal Seesaw.” arXiv:2212.09697.
  * [77] R. Dcruz, “Flavor Physics Constraints on Left-Right Symmetric Models with Universal Seesaw.” arXiv:2301.10786.
  * [78] K. Fujikura, K. Harigaya, Y. Nakai, and I. R. Wang, “Electroweak-like baryogenesis with new chiral matter,” JHEP 07 (2021) 224 [arXiv:2103.05005]. [Erratum: JHEP 12, 192 (2021), Erratum: JHEP 1, 156 (2022), Erratum: JHEP 01, 156 (2022)].
  * [79] K. Harigaya and I. R. Wang, “Baryogenesis in a Parity Solution to the Strong CP Problem.” arXiv:2210.16207.
  * [80] C. Ford, I. Jack, and D. R. T. Jones, “The Standard model effective potential at two loops,” Nucl. Phys. B 387 (1992) 373–390 [hep-ph/0111190]. [Erratum: Nucl.Phys.B 504, 551–552 (1997)].
  * [81] J. R. Espinosa and R.-J. Zhang, “Complete two loop dominant corrections to the mass of the lightest CP even Higgs boson in the minimal supersymmetric standard model,” Nucl. Phys. B 586 (2000) 3–38 [hep-ph/0003246].
  * [82] S. P. Martin, “Two Loop Effective Potential for a General Renormalizable Theory and Softly Broken Supersymmetry,” Phys. Rev. D 65 (2002) 116003 [hep-ph/0111209].
  * [83] S. P. Martin and H. H. Patel, “Two-loop effective potential for generalized gauge fixing,” Phys. Rev. D 98 (2018) 076008 [arXiv:1808.07615].