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:  CAS Key
Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese
Academy of Sciences, Beijing 100190, China

# Closer look at the matching condition for radiative QCD $\theta$ parameter

Tatsuya Banno a,b,c Junji Hisano b,d Teppei Kitahara a and Naohiro Osamura
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we scrutinize a radiatively generated QCD $\theta$ parameter at
the two-loop level based on both full analytical loop functions with the Fock-
Schwinger gauge method and the effective field theory approach, using
simplified models. We observe that the radiatively generated $\theta$
parameters at the low energy scale precisely match between them. It provides
validity to perturbative loop calculations of the QCD $\theta$ parameter with
the Fock-Schwinger gauge method. Furthermore, it is also shown that the
ordinary Fujikawa method for the radiative $\theta$ parameter by using
$\bar{\theta}=-\text{arg}\,\text{det}\,\mathcal{M}_{q}^{\rm loop}$ does not
cover all contributions in the simplified models. But, we also find that when
there is a scale hierarchy in $CP$-violating sector, evaluation of the
Fujikawa method is numerically sufficient. As an application, we calculate the
radiative $\theta$ parameter at the two-loop level in a slightly extended
Nelson-Barr model, where the spontaneous $CP$ violation occurs to solve the
strong $CP$ problem. It is found a part of the radiative $\theta$ parameters
cannot be described by the Fujikawa method.

###### Keywords:

$CP$ violation, Electric Dipole Moments

††preprint: IPMU23-0042

## 1 Introduction

In the renormalizable QCD Lagrangian, the following terms are $P$-odd and
$T$-odd, and then $CP$-odd terms under the $CPT$ invariance,

$\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}\,.$
(1.5)

Here, $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 known that one can turn off the first
terms with the axial rotation of quarks in order to get the physical QCD
$\theta$ parameter, $\bar{\theta}\,(\equiv\theta_{G}-\sum_{q}\theta_{q})$,
defined as

$\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}\,.$
(1.10)

This QCD $\theta$ parameter would be expected to be $O(1)$ naturally while it
is constrained by the neutron electric dipole moments (EDM) measurements as
$|\bar{\theta}|\lesssim 10^{-10}$ Abel:2020pzs ; Liang:2023jfj
.#1#1#1Recently, EDMs of paramagnetic molecules such as ThOACME:2018yjb  or
HfF+Roussy:2022cmp  are significantly improved, and the constraint on the QCD
$\theta$ parameter from $CP$-odd semileptonic interactions by these
experiments reaches $|\bar{\theta}|<10^{-8}$ Flambaum:2019ejc  if the source
of $CP$-violation is dominated by the hadronic sector. This serious
naturalness problem is called the strong $CP$ problem.

Some models have been proposed to solve the strong $CP$ problem. One of them
is the Peccei-Quinn (PQ) mechanism Peccei:1977hh ; Weinberg:1977ma ;
Wilczek:1977pj . The mechanism works well to solve the problem since the
strong dynamics tunes the QCD $\theta$ parameter to be zero with the axion
field, while an introduction of the PQ global symmetry suffers from the
quality problem of the global symmetry in quantum gravity Kamionkowski:1992mf
; Holman:1992us ; Barr:1992qq . The second option is to forbid the QCD
$\theta$ term with the $CP$ Nelson:1983zb ; Barr:1984qx ; Barr:1984fh  or $P$
symmetry Beg:1978mt ; Mohapatra:1978fy ; Babu:1989rb ; Barr:1991qx . Those
discrete symmetries are broken in the Standard Model (SM), so that the
extensions of the SM with those symmetries are mandatory. In the left-right
symmetric extension of the SM, the generalized $P$ symmetry forbids the QCD
$\theta$ term (and also the imaginary parts of quark masses) at the tree
level. After the spontaneous breaking of the $P$ symmetry to match the SM, the
QCD $\theta$ parameter is radiatively generated. Also, it is known that there
is no quality problem in the mechanism Craig:2020bnv . On the other hand, in
the Nelson-Barr models, the $CP$ symmetry is imposed so that the QCD $\theta$
term is forbidden at the tree level. Then, the $CP$ symmetry is spontaneously
broken to generate the $CP$ phase in the CKM matrix, and the model is tuned to
keep the QCD $\theta$ parameter zero at the tree level, but it is generated
radiatively.

In those models of the second option to solve the strong $CP$ problem, the QCD
$\theta$ parameter is radiatively generated and a finite neutron EDM is
predicted even if the energy scale of the symmetry breaking is much higher
than the weak scale Hall:2018let . Thus, the precise evaluation of the
radiative $\theta$ parameter is important to judge whether the models are
viable under the current constraint from the neutron EDM measurements and also
whether they are testable in future EDM experiments for the neutron and other
developed systems. In past literature, e.g., Refs. Babu:1989rb ; Hall:2018let
; Craig:2020bnv ; deVries:2021pzl , the radiative corrections to the QCD
$\theta$ parameter have been evaluated by using the radiative corrections to
the imaginary parts of the colored particle masses Ellis:1978hq ,

$\displaystyle\bar{\theta}=-\textrm{arg}\,\textrm{det}\,\mathcal{M}_{q}^{\rm
loop}\,,$ (1.11)

where $\mathcal{M}_{q}^{\rm loop}$ is the renormalized (up- and down-type)
colored fermion mass matrix. This relation is based on the Fujikawa method
Fujikawa:1979ay  or the Adler-Bell-Jackiw anomaly with Adler-Bardeen theorem
Adler:1969gk ; Bell:1969ts ; Adler:1969er ; Fujikawa:1979ay ; tHooft:1972tcz .

However, it is shown in Ref. Hisano:2023izx  that the radiative correction to
the QCD $\theta$ parameter can be evaluated directly and perturbatively, not
via the imaginary parts of quark masses, with the Feynman diagrammatic
calculations using the Fock-Schwinger gauge method Fock:1937dy ;
Schwinger:1951nm ; key86364857 ; Cronstrom:1980hj ; Shifman:1980ui ;
Dubovikov:1981bf ; Novikov:1983gd . The QCD $\theta$ term is equivalent to a
total derivative term, and hence the radiative correction to the QCD $\theta$
parameter is not perturbatively generated if the translational invariance is
maintained. The operator-Schwinger method Novikov:1983gd  using the log-
determinant of the Dirac operators works for the evaluation of the radiative
corrections to the action with keeping the translational symmetry. However,
the operator-Schwinger method makes it difficult to go beyond the one-loop
corrections to the QCD $\theta$ parameter. On the other hand, the Fock-
Schwinger gauge method is applicable to the higher-loop corrections
Abe:2017sam ; Hisano:2023izx . In fact, the Fock-Schwinger gauge violates the
translational invariance, while it is shown that the invariance is recovered
in the physical overvalues Nikolaev:1981ff ; Nikolaev:1982rq . Some of the
authors applied the Fock-Schwinger gauge method for evaluation of the QCD
$\theta$ parameter in the minimal left-right model Hisano:2023izx . It is
found that the two-loop contribution to the QCD $\theta$ parameter vanishes in
the model. The order of magnitude of the three-loop contribution from the
diagrams is evaluated and it is found that it can be as large as the current
experimental bound of the neutron EDM measurement.

In this paper, we derive the structure of the radiative correction to the QCD
$\theta$ parameter at two-loop level. We compare the full two-loop corrections
with those evaluated by the Fujikawa method taking into account only the
radiative correction to the colored fermion masses at one-loop level in Eq.
(1.11). In Sec. 2, we review several calculations related to the QCD $\theta$
parameter at one-loop level by using $CP$-violating low-energy effective
theory. In Sec. 3, we study minimal models at two-loop level, in which colored
fermion(s) and a real scalar have $CP$-violating Yukawa coupling(s) inducing
the radiative QCD $\theta$ parameter. We find that if the colored fermions and
the real scalar have masses comparable to each others, the Fujikawa method in
Eq. (1.11) does not work and we need to calculate the two-loop contribution to
the QCD $\theta$ parameter directly.

In Sec. 4, we discuss the two-loop corrections to the QCD $\theta$ parameter
in an extension of the Nelson-Barr model, as a concrete example of UV model.
In the model, we introduce the four-point couplings of the SM Higgs ($H$) and
the spontaneous $CP$-violating scalar ($\Sigma_{a}$, $a=1,\cdots,N_{\Sigma}$),
and a real scalar ($S$) fields, responsible to the real vector-like fermion
mass. Both interactions of $|H|^{2}\Sigma_{a}\Sigma_{b}^{\ast}$ and
$S^{2}\Sigma_{a}\Sigma_{b}^{\ast}$, induce the QCD $\theta$ parameter at two-
loop level. The former is determined by the radiative correction to the SM
quark masses. On the other hand, the latter depends on the heavy particle mass
spectrum, and it may not be determined by the radiative correction to the
fermion mass. Section 5 is devoted to conclusions. Details of the two-loop
calculations are given in the Appendix.

## 2 One-loop calculations for matching condition of QCD $\theta$ parameter

In this section, we evaluate the QCD $\theta$ parameter by integrating out a
colored fermion $q$, whose mass and interactions are $CP$ violating, in low-
energy effective field theories at one-loop level. The general effective
Lagrangian of the colored fermion is given in QCD up to dimension-five
operators as

$\displaystyle{\cal L}_{\rm eff}$
$\displaystyle=\bar{q}\left[i{\ooalign{\hfil/\hfil\crcr$D$}}-\left(m_{q}^{\ast}P_{L}+m_{q}P_{R}\right)\right]q-\frac{1}{2}g_{s}{\mu}_{q}\bar{q}\left(\sigma\cdot
G\right)q-\frac{i}{2}g_{s}{d}_{q}\bar{q}\left(\sigma\cdot
G\right)\gamma_{5}q\,,$ (2.3)

where $D_{\mu}=\partial_{\mu}+ig_{s}T^{a}G_{\mu}^{a}$,
$T_{a}=\frac{1}{2}\lambda_{a}$ ($\lambda_{a}$ are the Gell-Mann matrices),
$(\sigma\cdot G)=\sigma^{\mu\nu}T^{a}G_{\mu\nu}^{a}$,
$\sigma_{\mu\nu}=\frac{i}{2}[\gamma_{\mu},\gamma_{\nu}]$, and
$P_{L/R}=(1\mp\gamma_{5})/2$. If the $CP$-violating phase of the colored
fermion mass, $m_{q}\equiv|m_{q}|\exp(i\theta_{q})$, is removed by the axial
rotation (chiral rotation) $q\to\exp(-\frac{i}{2}\theta_{q}\gamma_{5})q$, one
obtains Adler:1969gk ; Bell:1969ts ; Adler:1969er ; Fujikawa:1979ay ;
tHooft:1972tcz

$\displaystyle\mathcal{L}_{\rm
eff}=\bar{q}\left(i{\ooalign{\hfil/\hfil\crcr$D$}}-|m_{q}|\right)q-\frac{1}{2}g_{s}\tilde{\mu}_{q}\bar{q}\left(\sigma\cdot
G\right)q-\frac{i}{2}g_{s}\tilde{d}_{q}\bar{q}\left(\sigma\cdot
G\right)\gamma_{5}q-\theta_{q}\frac{\alpha_{s}}{8\pi}G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}\,,$
(2.6)

with

$\displaystyle\tilde{\mu}_{q}=\cos(\theta_{q})\mu_{q}+\sin(\theta_{q})d_{q}\,,\quad\tilde{d}_{q}=-\sin(\theta_{q})\mu_{q}+\cos(\theta_{q})d_{q}\,,$
(2.7)

where $\tilde{\mu}_{q}$ and $\tilde{d}_{q}$ correspond to the (anomalous)
chromo-magnetic dipole moment (CMDM) and chromo-electric dipole moment
(CEDM).#2#2#2Note that a nonrelativistic term of the CMDM is included in the
Dirac operator. Here, we ignore the dimension-six four-Fermi
operators.#3#3#3In this paper, we also discarded the QED interactions. They
are mixed with the dipole moment operators at one-loop level renormalization-
group evolution if they are present Hisano:2012cc . The renormalization-group
evolution effect will be taken into account in the next section.

We integrate out the colored fermion with the operator Schwinger method to
evaluate the QCD $\theta$ parameter at one-loop level. The method is
transparent compared with the Fock-Schwinger gauge method when considering the
one-loop corrections Hisano:2023izx ; Novikov:1983gd . They must be coincident
with each other when we calculate the gauge invariant quantities.

First, we consider the contribution of the colored fermion mass phase to the
QCD $\theta$ parameter. It is shown in Ref. Hisano:2023izx  that in the
operator-Schwinger method, the one-loop correction to the effective action
$\Delta S$ after integrating out the colored fermion field is given as

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

where $P_{\mu}\equiv iD_{\mu}$. The derivative of $\Delta S$ over $m_{q}$ for
$P_{R}$ is Hisano:2023izx ; Novikov:1983gd

$\displaystyle\frac{d}{dm_{q}}\Delta S$
$\displaystyle=i\operatorname{Tr}\left[\frac{1}{P^{2}-|m_{q}|^{2}-\frac{1}{2}g_{s}(\sigma\cdot
G)}m_{q}^{\ast}P_{R}\right]\,,$ (2.11)

where $\left[P_{\mu},P_{\nu}\right]=-ig_{s}G_{\mu\nu}$ is used. Since the
$G\tilde{G}$ term includes the Levi-Civita tensor, the QCD $\theta$ term comes
from the second order of the perturbation of $\frac{1}{2}g_{s}(\sigma\cdot G)$
in the denominator in the right-handed side of Eq. (2.11),

$\displaystyle\frac{d}{dm_{q}}\Delta S$ $\displaystyle\supset
i\operatorname{Tr}\left[\frac{1}{P^{2}-|m_{q}|^{2}}\left(\frac{1}{2}g_{s}(\sigma\cdot
G)\right)\frac{1}{P^{2}-|m_{q}|^{2}}\left(\frac{1}{2}g_{s}(\sigma\cdot
G)\right)m_{q}^{\ast}P_{R}\right]$ $\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^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}+\cdots$
$\displaystyle=i\int
d^{4}x\frac{g_{s}^{2}}{64\pi^{2}}\frac{1}{m_{q}}G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}+\cdots\,,$
(2.12)

where ${\rm Tr}(T^{a}T^{b})=(1/2)\delta^{ab}$ is used. Here, $P^{2}$ is
replaced with $(i\partial)^{2}$, and it is integrated in the 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{\cal L}$
$\displaystyle=-i\frac{\alpha_{s}}{16\pi}\log\frac{m_{q}^{\ast}}{m_{q}}G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}$
$\displaystyle=-\theta_{q}\frac{\alpha_{s}}{8\pi}G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}\,.$
(2.13)

As expected, the integration of the colored fermion gives the contribution to
QCD $\theta$ term, consistent with the Fujikawa method in Eq. (2.6).

Next, we show the contribution of the dipole moments, $\tilde{\mu}_{q}$ and
$\tilde{d}_{q}$, in Eq. (2.6). The CEDM is $T$-odd and $P$-odd, and it also
contributes to the QCD $\theta$ parameter as shown below. When the phase of
the colored fermion mass is removed by the chiral rotation, the physical CEDM
$\tilde{d}_{q}$ is given by
$\tilde{d}_{q}=\cos(\theta_{q})d_{q}-\sin(\theta_{q})\mu_{q}$ in Eq. (2.7).
The evaluation of the physical CEDM contribution to the QCD $\theta$ parameter
is more straightforward in the operator-Schwinger method. The one-loop
correction to the effective action after integrating out the colored fermion
field is

$\displaystyle\Delta S$
$\displaystyle=-i\operatorname{Tr}\log\left|\left|{\ooalign{\hfil/\hfil\crcr$P$}}-|m_{q}|-\frac{i}{2}g_{s}\tilde{d}_{q}(\sigma\cdot
G)\gamma_{5}\right|\right|\,$ (2.16)
$\displaystyle\supset-i\operatorname{Tr}\left[\frac{1}{P^{2}-|m_{q}|^{2}-\frac{1}{2}g_{s}(\sigma\cdot
G)}|m_{q}|\left(-\frac{i}{2}g_{s}\tilde{d}_{q}(\sigma\cdot
G)\gamma_{5}\right)\right]\,$
$\displaystyle=-|m_{q}|\tilde{d}_{q}\operatorname{Tr}\left[\frac{1}{P^{2}-|m_{q}|^{2}}\left(\frac{1}{2}g_{s}(\sigma\cdot
G)\right)\frac{1}{P^{2}-|m_{q}|^{2}}\left(\frac{1}{2}g_{s}(\sigma\cdot
G)\right)\gamma_{5}\right]\cdots\,$
$\displaystyle=2|m_{q}|\tilde{d}_{q}\times(16\pi^{2}\mu^{2\epsilon})\int\frac{d^{d}p}{(2\pi)^{d}}\frac{1}{(p^{2}-|m_{q}|^{2})^{2}}\times\int
d^{4}x\frac{\alpha_{s}}{8\pi}G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}$
$\displaystyle=2|m_{q}|\tilde{d}_{q}\times\left(\frac{1}{\bar{\epsilon}}-\log\frac{|m_{q}|^{2}}{\mu^{2}}\right)\times\int
d^{4}x\frac{\alpha_{s}}{8\pi}G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}\,.$ (2.17)

In the second line, we keep the leading term of the CEDM in the expansion. In
the fourth line, $P^{2}$ is replaced by $(i\partial)^{2}$, and it is
integrated out in the momentum space again. Since the integral includes UV
divergence, we introduce the dimensional regularization $d=4-2\epsilon$ with
the renormalization scale $\mu$. If the $\overline{\text{MS}}$ scheme is
adopted, $1/\bar{\epsilon}\,(\equiv 1/\epsilon-\gamma_{E}+\log 4\pi$ and
$\gamma_{E}$ is the Euler’s constant) is subtracted by the counter term. It is
found that the QCD $\theta$ parameter depends on the renormalization scale as
Jenkins:2017dyc

$\displaystyle\mu\frac{d\theta}{d\mu}$
$\displaystyle=4|m_{q}|\tilde{d}_{q}\,,$ (2.18)

and the threshold correction at the colored fermion mass ($\mu=|m_{q}|$) due
to the CEDM vanishes.

It is known that the CEDMs of the colored fermions generate the $T$-odd and
$P$-odd higher-dimensional gluon operators, such as the Weinberg operator with
dimension six, $G^{2}\tilde{G}$, in addition to the QCD $\theta$ parameter
when the colored fermions are integrated out Boyd:1990bx ; Braaten:1990gq ;
Chang:1990jv ; Dine:1990pf ; Abe:2017sam . On the other hand, the colored
fermion mass phases induce only the QCD $\theta$ parameter, but not the
higher-dimensional gluon operators. The axial rotations of the colored
fermions turn off the colored fermion mass phases and simultaneously generate
the QCD $\theta$ term in Eq. (2.6). This can be derived with the Fujikawa
method Fujikawa:1979ay  or the explicit calculation of the above one-loop
colored fermion diagram Hisano:2023izx . Now, we consider whether $T$-odd and
$P$-odd higher-dimensional gluon operators, such as the Weinberg operator, are
generated or not, when the quark mass phases are non vanishing.

In the Fujikawa method, the higher-dimensional gluon operators are not
generated from the Jacobian of the chiral transformation Fujikawa:1979ay .
However, an UV cut-off parameter has to be introduced in order to define the
anomalous Jacobian, which includes only the QCD $\theta$ term in the limit of
the infinite UV cut-off parameter. Therefore, it might be unclear whether the
higher-dimensional terms are generated when the UV cut-off is finite in the
Fujikawa method, and we would like to derive the absence of the higher-
dimensional gluon operators via the explicit one-loop calculation without the
introduction of the UV cut-off parameter. Here, by using the operator-
Schwinger method we show a theorem that the $T$-odd and $P$-odd higher-
dimensional gluon operators are not generated at one-loop level even if the
quark mass phases are not vanishing, whose diagram is shown in Fig. 1.

Figure 1:  The one-loop babble diagram contributing to the radiative
corrections to the $T$-odd and $P$-odd higher-dimensional gluon operators. It
is shown that they are not induced from the quark mass phases.

In the operator-Schwinger method, we start from the derivative of $\Delta S$
over $m_{q}$ for $P_{R}$ in Eq. (2.11). Then, by using the covariant
derivative expansion (CDE) Henning:2014wua ; Henning:2016lyp , the trace over
the functional space can be evaluated by

$\displaystyle\frac{d}{dm_{q}}\Delta S=i\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\operatorname{tr}\left[\frac{1}{p^{2}-|m_{q}|^{2}-(2ip^{\mu}D_{\mu}+D^{2}+\frac{1}{2}g_{s}(\sigma\cdot
G))}m_{q}^{\ast}P_{R}\right]$ (2.19) $\displaystyle=i\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\operatorname{tr}\left\\{\left[\frac{1}{p^{2}-|m_{q}|^{2}}\left(\frac{1}{2}g_{s}(\sigma\cdot
G)\right)\frac{1}{p^{2}-|m_{q}|^{2}}\left(\frac{1}{2}g_{s}(\sigma\cdot
G)\right)\frac{1}{p^{2}-|m_{q}|^{2}}m_{q}^{*}P_{R}\right]\right.$
$\displaystyle\quad+\operatorname{tr}\left[\frac{1}{p^{2}-|m_{q}|^{2}}\left(\frac{1}{2}g_{s}(\sigma\cdot
G)\right)\frac{1}{p^{2}-|m_{q}|^{2}}\left(\frac{1}{2}g_{s}(\sigma\cdot
G)\right)\frac{1}{p^{2}-|m_{q}|^{2}}\left(\frac{1}{2}g_{s}(\sigma\cdot
G)\right)\frac{1}{p^{2}-|m_{q}|^{2}}m_{q}^{*}P_{R}\right]$
$\displaystyle\quad+\mbox{(terms including two $G$s and two $D$s in the
numerator)}\Biggr{\\}}$ $\displaystyle\quad+\cdots\,,$ (2.20)

where the first line in Eq. (2.20) induces $G\tilde{G}$ term, while the
remaining ones we obtain

$\displaystyle\frac{d}{dm_{q}}\Delta S$ $\displaystyle=\int
d^{4}x\left(\frac{i}{64\pi^{2}}\frac{g_{s}^{2}}{m_{q}}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}-\frac{i}{192\pi^{2}}r_{GGG}\frac{g_{s}^{3}}{m_{q}|m_{q}|^{2}}f^{abc}G_{\mu\nu}^{a}G^{b\nu}_{\rho}\tilde{G}^{c\rho\mu}\right)\,,$
(2.21)

where $r_{GGG}$ is a $\mathcal{O}(1)$ real number (for example, $r_{GGG}=1$
comes from the second line of Eq. (2.20)). In the above equation, $G\tilde{G}$
term comes from the second order of $\frac{1}{2}g_{s}(\sigma\cdot G)$ in the
denominator and its coefficient is proportional to $1/m_{q}$. On the other
hand, $G^{2}\tilde{G}$ term comes from the expansions of the covariant
derivatives and/or $\frac{1}{2}g_{s}(\sigma\cdot G)$ and its coefficient is
proportional to $1/(m_{q}|m_{q}|^{2})$. The coefficients are determined with
the dimensions of the operators up to the $O(1)$ coefficients. Similarly, the
derivative of $\Delta S$ over $m_{q}^{\ast}$ for $P_{L}$ is

$\displaystyle\frac{d}{dm_{q}^{\ast}}\Delta S$
$\displaystyle=i\operatorname{Tr}\left[\frac{1}{P^{2}-|m_{q}|^{2}-\frac{1}{2}g_{s}\sigma
G}m_{q}P_{L}\right]$ $\displaystyle=-\int
d^{4}x\left(\frac{i}{64\pi^{2}}\frac{g_{s}^{2}}{m_{q}^{\ast}}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}-\frac{i}{192\pi^{2}}r_{GGG}\frac{g_{s}^{3}}{m_{q}^{\ast}|m_{q}|^{2}}f^{abc}G_{\mu\nu}^{a}G^{b\nu}_{\rho}\tilde{G}^{c\rho\mu}\right)\,.$
(2.22)

Therefore, the overall factor $1/m_{q}$ in Eq. (2.21) is replaced by
$-1/m_{q}^{\ast}$, where the relative phase comes from the difference of the
projection operators $P_{L/R}$.

The one-loop corrections to the effective action can be determined by adding
the integrals $d\Delta S/dm_{q}$ with $m_{q}$ and $d\Delta S/dm_{q}^{\ast}$
with $m_{q}^{\ast}$,

$\displaystyle\Delta S$ $\displaystyle\supset\int
d^{4}x\left(\frac{ig_{s}^{2}}{64\pi^{2}}\log
m_{q}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}+\frac{ig_{s}^{3}}{192\pi^{2}}r_{GGG}\frac{1}{|m_{q}|^{2}}f^{abc}G_{\mu\nu}^{a}G^{b\nu}_{\rho}\tilde{G}^{c\rho\mu}\right)$
$\displaystyle\quad+\int d^{4}x\left(-\frac{ig_{s}^{2}}{64\pi^{2}}\log
m_{q}^{\ast}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}-\frac{ig_{s}^{3}}{192\pi^{2}}r_{GGG}\frac{1}{|m_{q}|^{2}}f^{abc}G_{\mu\nu}^{a}G^{b\nu}_{\rho}\tilde{G}^{c\rho\mu}\right)$
$\displaystyle=\int
d^{4}x\frac{ig_{s}^{2}}{64\pi^{2}}\log\frac{m_{q}}{m_{q}^{\ast}}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}$
$\displaystyle=-\int
d^{4}x\theta_{q}\frac{\alpha_{s}}{8\pi}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}\,.$
(2.23)

Therefore, we find that only the QCD $\theta$ term is generated with the
coefficient proportional to the quark mass phase, and the Weinberg operator is
not induced from Fig. 1.

We conclude that such a cancellation occurs for all the $T$-odd and $P$-odd
higher-dimensional gluon operators with dimensions ($D$) larger than four. In
general, $d\Delta S/dm_{q}$ includes all the $T$-odd and $P$-odd dimension-$D$
gluon operators from the expansions of the covariant derivatives and/or
$\frac{1}{2}g_{s}(\sigma\cdot G)$ in the denominator, and their coefficients
are proportional to $m_{q}^{\ast}/|m_{q}|^{(D-2)}$. On the other hand,
$d\Delta S/dm_{q}^{\ast}$ has the gluon operators with the coefficients
proportional to $-m_{q}/|m_{q}|^{(D-2)}$. By integrating $d\Delta S/dm_{q}$
with $m_{q}$ and integrating $d\Delta S/dm_{q}^{\ast}$ with $m_{q}^{\ast}$, it
is found that the coefficients of the $T$-odd and $P$-odd higher-dimensional
gluon operators in Fig. 1 are canceled, while only $D=4$ term survives. Thus,
we showed that the quark mass phases do not generate the $T$-odd and $P$-odd
higher-dimensional gluon operators at the one-loop level.

## 3 Two-loop contribution to QCD $\theta$ parameter in simplified models

In this section, using the Fock-Schwinger gauge method, we evaluate the two-
loop corrections to the QCD $\theta$ parameter in simple models, in which
fermion(s) and a real scalar have $CP$-violating Yukawa interaction(s). We
find that when the fermions and the scalar have masses comparable to each
other, the Fujikawa method does not work in the evaluation of the QCD $\theta$
parameter even if taking into account the radiative correction to the
imaginary parts of the fermion masses at one-loop level.

### 3.1 One flavor case

First, we evaluate the two-loop contributions to the QCD $\theta$ parameter in
a simplified model in which a $CP$-violating Yukawa coupling $y_{q}$ with
colored fermion $q$ is contained, where the number of the flavor is one for
simplicity. Let us consider the following Lagrangian,

$\displaystyle-\mathcal{L}=\frac{1}{2}m_{\phi}^{2}\phi^{2}+\bar{q}\left[\text{Re}\left(m_{q}\right)+\text{Im}\left(m_{q}\right)i\gamma_{5}\right]q+y_{q}\bar{q}P_{R}q\phi+\text{h.c.}\,,$
(3.1)

where $\phi$ is a neutral scalar. Note that the bare QCD $\theta$ parameter
$\theta_{G}$ is set to be zero, cf. Eq. (1.5). In the following, we assume
$\text{Im}(m_{q})\ll\text{Re}(m_{q})$ at the tree level in order to
investigate the structure of the two-loop corrections using the low-energy
effective theories. On the other hand, we keep the one-loop correction formula
since the one-loop diagrams include the counter terms to cancel divergence of
subdiagrams in the two-loop diagrams and the combined results of one- and two-
loop corrections are renormalization-scale independent.

Figure 2:  The one-loop and two-loop babble diagrams which contribute to the
radiative corrections to the QCD $\theta$ parameter. The one-loop diagram
gives a simple result $\delta\theta|_{\rm
1L}=-\text{Im}(m_{q})/\text{Re}(m_{q})$. For the two-loop diagrams, the first
one generates $I_{(2;2)}$ loop function in Eq. (3.4), while the others
$2\bar{I}_{(3;1)}$ in total.

When one sets $\theta_{G}=0$, the physical QCD $\theta$ parameter can be
perturbatively calculated by the Fock-Schwinger gauge method Hisano:2023izx ;
Novikov:1983gd . The radiative corrections $\delta\theta$ are derived from
one-loop (1L) and two-loop diagrams (2L) given in Fig. 2,
$\delta\theta=\delta\theta|_{\rm 1L}+\delta\theta|_{\rm 2L}$, as
Hisano:2023izx

$\displaystyle\bar{\theta}$ $\displaystyle=\delta\theta=\delta\theta|_{\rm
1L}+\delta\theta|_{\rm 2L}+\cdots\,,$ (3.2) $\displaystyle\delta\theta|_{\rm
1L}$ $\displaystyle=-\frac{\text{Im}(m_{q})}{\text{Re}(m_{q})}\,,$ (3.3)
$\displaystyle\delta\theta|_{\rm 2L}$ $\displaystyle=\frac{{\rm
Im}\left(y_{q}^{2}\right)}{16\pi^{2}}{\left[\text{Re}(m_{q})\right]^{2}}\left[I_{(2;2)}(m_{q}^{2};m_{q}^{2};m_{\phi}^{2})+2\bar{I}_{(3;1)}(m_{q}^{2};m_{q}^{2};m_{\phi}^{2})\right]\,,$
(3.4)

with#4#4#4An accurate numerator of the two-loop correction is
$\text{Im}(y_{q}m_{q}y_{q}m_{q})$, and we expand it with
$\text{Im}(m_{q})\ll\text{Re}(m_{q})$.

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

where $Q^{2}\equiv 4\pi\mu^{2}{\rm e}^{-\gamma_{E}}$ and the dimensional
regularization is used with the renormalization scale. Here, the mass squared
for the colored-quark masses in the loop functions is understood as
$m_{i}^{2}={\text{Re}(m_{i})}^{2}+{\text{Im}(m_{i})}^{2}$. The two-loop
functions $I_{(2;2)}$ and $I_{(3;1)}$ are given in Appendix A.

The loop-function $I_{(3;1)}$ corresponds to the last two two-loop diagrams in
Fig. 2, whose subdiagram gives a one-loop UV divergence $\propto
y_{q}^{2}m_{q}/\bar{\epsilon}$, and this divergence must be canceled by a
counter term of the mass renormalization. We adopt the $\overline{\rm MS}$
scheme for the mass renormalization. The second term in
$\bar{I}_{(3;1)}(x_{1};x_{2};x_{3})$, subtracting a $1/\bar{\epsilon}$ pole,
corresponds to the counter term contribution, and the factor $(1-\epsilon\log
x_{1}/Q^{2})$ comes from the higher-order term of $\epsilon$ in the
dimensional regularization (see Eq. (2.15) in Ref. Hisano:2023izx ).#5#5#5The
second term in $\bar{I}_{(3;1)}(x_{1};x_{2};x_{3})$ corresponds to the
divergence contribution $\bar{I}_{\epsilon(3;1)}$ defined in Ref.
Hisano:2023izx  when one adopts the ${\rm MS}$ scheme for the mass
renormalization. As a result, $Q^{2}$ in ${F}_{0}(p^{2},x_{1},x_{2})$ in Ref.
Hisano:2023izx  is replaced with $\mu^{2}$ in $F_{0}(p^{2},x_{1},x_{2})$ in
Eq. (3.9).  Then, $\bar{I}_{(3;1)}(x_{1},x_{2},x_{3})$ is UV and IR finite
function Hisano:2023izx .

The above formula for the radiative corrections to the QCD $\theta$ parameter
at two-loop level becomes simpler when the colored fermion is much lighter (or
heavier) than the scalar $\phi$, since the effective theory description works.
Let us check such hierarchical situations.

First, when the scalar particle mass is much heavier than the fermion mass
$m_{\phi}\gg m_{q}$, the two-loop contributions in Eq. (3.4) can be simplified
to

$\displaystyle{I}_{(2;2)}(m_{q}^{2};m_{q}^{2};m_{\phi}^{2})+2\bar{I}_{(3;1)}(m_{q}^{2};m_{q}^{2};m_{\phi}^{2})$
$\displaystyle\quad\quad\simeq\frac{1}{m_{q}^{2}}\left(1-\log\frac{m_{\phi}^{2}}{\mu^{2}}\right)+\frac{\pi^{2}+12\log\frac{m_{q}^{2}}{m_{\phi}^{2}}+3\log^{2}\frac{m_{q}^{2}}{m_{\phi}^{2}}}{3m_{\phi}^{2}}\,.$
(3.6)

From the viewpoint of the effective field theory, first of all, the scalar
$\phi$ is integrated out at a scale of $\mu=m_{\phi}$, then the following
$CP$-violating effective interactions may be obtained up to the dimension-six
operators Hisano:2012cc ,

$\displaystyle-\mathcal{L}_{\rm eff}$
$\displaystyle=\bar{q}\left[\text{Im}\left(m_{q}+\Delta
m_{q}\right)i\gamma_{5}\right]q$
$\displaystyle\quad+\frac{i}{2}g_{s}\tilde{d}_{q}\bar{q}\sigma^{\mu\nu}\gamma_{5}T^{a}qG_{\mu\nu}^{a}-C_{4}^{q}(\bar{q}q)(\bar{q}i\gamma_{5}q)-C_{5}^{q}(\bar{q}\sigma^{\mu\nu}q)(\bar{q}i\sigma_{\mu\nu}\gamma_{5}q)$
$\displaystyle\quad+\Delta\theta|_{\rm
th}\frac{\alpha_{s}}{8\pi}G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}+\frac{1}{3}\omega
f^{abc}G^{a}_{\mu\nu}G^{b\nu}_{\rho}\tilde{G}^{c\rho\mu}\,,$ (3.7)

where $f^{abc}$ is the structure constant.

The radiative correction to the imaginary part of the fermion mass at
$\mu\simeq m_{\phi}$ is given as

$\displaystyle\text{Im}\left[\Delta m_{q}(\mu)\right]$
$\displaystyle=\frac{{\rm
Im}\left(y_{q}^{2}\right)}{16\pi^{2}}{\text{Re}(m_{q})}{F}_{0}(0,m_{q}^{2},m^{2}_{\phi})\,,$
(3.8)

where

$\displaystyle{F}_{0}(p^{2},x_{1},x_{2})$
$\displaystyle=\int_{0}^{1}dz\log\frac{-z(1-z)p^{2}+zx_{1}+(1-z)x_{2}}{\mu^{2}}\,.$
(3.9)

The above radiative correction to the imaginary part of the fermion mass
corresponds to that in a limit of the zero external momentum, since the
external momentum is negligible compared with $m_{\phi}$. In the case, the
radiative correction is explicitly given as

$\displaystyle\text{Im}\left[\Delta m_{q}(\mu)\right]$
$\displaystyle=-\frac{{\rm
Im}\left(y_{q}^{2}\right)}{16\pi^{2}}{\text{Re}(m_{q})}\left[1-\frac{1}{m_{q}^{2}-m_{\phi}^{2}}\left(m_{q}^{2}\log\frac{m_{q}^{2}}{\mu^{2}}-m_{\phi}^{2}\log\frac{m_{\phi}^{2}}{\mu^{2}}\right)\right]$
(3.10) $\displaystyle\simeq\left\\{\begin{array}[]{ll}\frac{{\rm
Im}\left(y_{q}^{2}\right)}{16\pi^{2}}{\text{Re}(m_{q})}\log\frac{m_{\phi}^{2}}{\mu^{2}}&\quad\text{for\leavevmode\nobreak\
}m_{\phi}\approx m_{q}\,,\\\ -\frac{{\rm
Im}\left(y_{q}^{2}\right)}{16\pi^{2}}{\text{Re}(m_{q})}\left(1-\log\frac{m_{\phi}^{2}}{\mu^{2}}+\frac{m_{q}^{2}}{m_{\phi}^{2}}\log\frac{m_{q}^{2}}{m_{\phi}^{2}}\right)&\quad\text{for\leavevmode\nobreak\
}m_{\phi}\gg m_{q}\,.\end{array}\right.$ (3.13)

It is found by explicit calculation that the matching conditions for the
higher-dimensional operators in Eq. (3.7) at $\mu=m_{\phi}$ are#6#6#6 The
chiral rotation of $q$ to remove ${\rm Im}[\Delta m_{q}]$ leads to the
contribution to the CEDM from the CMDM in Eq. (2.7). However, since ${\rm
Im}[\Delta m_{q}]$ is at one-loop level, the contribution to the CEDM is at
two-loop level, and then we ignore it.

$\displaystyle\tilde{d}_{q}$ $\displaystyle=-\frac{{\rm
Im}\left(y_{q}^{2}\right)}{16\pi^{2}}\frac{3{\text{Re}(m_{q})}}{2m_{\phi}^{2}}\,,$
(3.14) $\displaystyle C_{4}^{q}$
$\displaystyle=\frac{\text{Im}\left(y_{q}^{2}\right)}{2m_{\phi}^{2}}\,,$
(3.15) $\displaystyle C_{5}^{q}$ $\displaystyle=\omega=0\,.$ (3.16)

The term $\Delta\theta|_{\rm th}$ in Eq. (3.7) is a threshold correction to
the QCD $\theta$ parameter, which will be identified by comparing with the
full result and the contributions from the effective operators in the
effective theory.

In the effective Lagrangian (3.7), the CEDM operators are mixed with the four-
Fermi operators via the renormalization-group equation Hisano:2012cc , and
then the $\tilde{d}_{q}$ depends on $\mu$ when $m_{q}<\mu<m_{\phi}$. Thus,
$\tilde{d}_{q}$ is approximately given as

$\displaystyle\tilde{d}_{q}(\mu)$ $\displaystyle=$ $\displaystyle-\frac{{\rm
Im}\left(y_{q}^{2}\right)}{16\pi^{2}}\frac{{\text{Re}(m_{q})}}{m_{\phi}^{2}}\left(\frac{3}{2}+\log\frac{\mu^{2}}{m_{\phi}^{2}}\right)\,.$
(3.17)

This CEDM operator induces the QCD $\theta$ parameter via the renormalization-
group equation Jenkins:2017dyc , as reviewed in Sec. 2.

Using these results in the effective theory, the radiatively generated QCD
$\theta$ parameter at the low-energy scale is evaluated as

$\displaystyle\delta\theta$ $\displaystyle=-\frac{\text{Im}(m_{q}+\Delta
m_{q})}{\text{Re}(m_{q})}-2\int^{\log m^{2}_{q}}_{\log
m^{2}_{\phi}}\tilde{d}_{q}(\mu){\text{Re}(m_{q})}d(\log\mu^{2})+\left.\Delta\theta\right|_{\rm
th}$ (3.18)
$\displaystyle=-\frac{\text{Im}(m_{q})}{\text{Re}(m_{q})}+\frac{{\rm
Im}\left(y_{q}^{2}\right)}{16\pi^{2}}\left(1-\log\frac{m_{\phi}^{2}}{\mu^{2}}+\frac{m_{q}^{2}}{m_{\phi}^{2}}\log\frac{m_{q}^{2}}{m_{\phi}^{2}}\right)$
$\displaystyle\quad+\frac{{\rm
Im}\left(y_{q}^{2}\right)}{16\pi^{2}}\frac{m_{q}^{2}}{m_{\phi}^{2}}\left({3}\log\frac{m_{q}^{2}}{m_{\phi}^{2}}+\log^{2}\frac{m_{q}^{2}}{m_{\phi}^{2}}\right)+\left.\Delta\theta\right|_{\rm
th}\,.$ (3.19)

This result is consistent with the QCD $\theta$ parameter directly calculated
from the diagrams in Fig. 2 (see Eq. (3.6)),

$\displaystyle\delta\theta$
$\displaystyle=\left.\delta\theta\right|_{1L}+\left.\delta\theta\right|_{2L}$
$\displaystyle=-\frac{\text{Im}(m_{q})}{\text{Re}(m_{q})}+\frac{{\rm
Im}\left(y_{q}^{2}\right)}{16\pi^{2}}\left[1-\log\frac{m_{\phi}^{2}}{\mu^{2}}+\frac{m_{q}^{2}}{3m_{\phi}^{2}}\left(\pi^{2}+12\log\frac{m_{q}^{2}}{m_{\phi}^{2}}+3\log^{2}\frac{m_{q}^{2}}{m_{\phi}^{2}}\right)\right]\,,$
(3.20)

where $\mu$ should be taken for the heaviest particle mass of the loop
diagrams, $\mu\simeq m_{\phi}$. In addition, by matching the above results one
obtains the threshold correction to the QCD $\theta$ parameter as

$\displaystyle\left.\Delta\theta\right|_{\rm th}=\frac{{\rm
Im}\left(y_{q}^{2}\right)}{48}\frac{{m_{q}^{2}}}{m_{\phi}^{2}}\,.$ (3.21)

Note that these effective theory matchings are valid only when $m_{\phi}\gg
2m_{q}$ Davydychev:1992mt .

In the above calculation, we implicitly assume that the light fermion is
heavier than the $\Lambda_{\rm QCD}$ scale. On the other hand, the light
quarks in the SM cannot be integrated out. The neutron EDM is evaluated in the
low-energy effective theory, which includes higher-dimensional operators of
light quarks in addition to those of gluons. Note that when fermions with mass
larger than $\Lambda_{\rm QCD}$, such as heavier quarks in the SM, are
integrated out, the fermion CEDMs induce the Weinberg operator as a threshold
correction and it also contributes to the neutron EDM Boyd:1990bx ;
Braaten:1990gq ; Chang:1990jv ; Dine:1990pf ; Abe:2017sam , though it does not
contribute to the QCD $\theta$ term. The contribution to the neutron EDM from
the induced Weinberg operator would be smaller than the contribution from
$\delta\theta$, suppressed by $\mathcal{O}(\Lambda_{\rm
QCD}^{2}/m_{\phi}^{2})$.

(a) $m_{\phi}\mathop{}_{\textstyle\sim}^{\textstyle>}m_{q}$

(b) $m_{\phi}\mathop{}_{\textstyle\sim}^{\textstyle<}m_{q}$

Figure 3:  The radiative corrections to the QCD $\theta$ parameter,
$\delta\theta$ which is normalized by $\textrm{Im}(y_{q}^{2})$, evaluated from
the full two-loop diagrams (Eq. (3.4)) are presented by the orange lines. For
comparison, the corrections evaluated from the one-loop correction to the
colored fermion mass in a limit of zero external momentum (Eq. (3.8)) with the
Fujikawa method are also shown by the blue lines. In the left panel (a), a
parameter region of $m_{\phi}\gtrsim m_{q}$ is shown with $\mu=m_{\phi}$. On
the other hand, in the right panel (b) $m_{\phi}\lesssim m_{q}$ with
$\mu=m_{q}$.

We show the numerical examination in Fig. (3(a)) in order to stress the
importance of the full two-loop diagram calculation in the case of
$m_{\phi}\gtrsim m_{q}$. Here, $\mu=m_{\phi}$, since the effective theory of
the colored fermion works after integrating out $\phi$, as mentioned above.
The blue line represents the contribution to the QCD $\theta$ parameter from
the imaginary part of the colored fermion mass in Eq. (3.8) with the Fujikawa
method in Eq. (1.11), while the orange one comes from the direct two-loop
diagram calculation in Eq. (3.4). The CEDM of the colored fermion contributes
to the QCD $\theta$ parameter, though it is suppressed by
$m_{q}^{2}/m^{2}_{\phi}$. The threshold correction in Eq. (3.21) is also
suppressed by $m_{q}^{2}/m^{2}_{\phi}$. As a result, the correction to the QCD
$\theta$ parameter is dominated by the contribution from the correction to the
imaginary part of the colored fermion mass when $m_{\phi}\gg m_{q}$, as
expected. However, it is also shown that when $m_{\phi}/m_{q}\lesssim 10$, the
estimation of the Fujikawa method does not work well.

On the other hand, when the scalar particle mass is quite light compared to
the fermion mass $m_{\phi}\ll m_{q}$, the two-loop functions can be simplified
to Davydychev:1992mt

$\displaystyle{I}_{(2;2)}(m_{q}^{2};m_{q}^{2};m_{\phi}^{2})+2\bar{I}_{(3;1)}(m_{q}^{2};m_{q}^{2};m_{\phi}^{2})$
$\displaystyle\quad\quad\simeq\frac{1}{m_{q}^{2}}\left(1-\log\frac{m_{q}^{2}}{\mu^{2}}\right)\,,$
(3.22)

where all two-loop diagrams are contributing.

From the viewpoint of the effective theory, this contribution corresponds to a
threshold correction to the QCD $\theta$ parameter when the fermion is
integrated out,

$\displaystyle\left.\Delta\theta\right|_{\rm th}=-\frac{{\rm Im}(m_{q})}{{\rm
Re}(m_{q})}+\frac{\text{Im}\left(y_{q}^{2}\right)}{16\pi^{2}}\left(1-\log\frac{m_{q}^{2}}{\mu^{2}}\right)\,,$
(3.23)

and hence $\mu\simeq m_{q}$. The $\mu$ dependence is canceled between the
first and the second terms. It is found that the second term coincides with
the one-loop correction to the imaginary part of the fermion mass (${\rm
Im}[\Delta m_{q}(\mu)]$) in a limit of the zero external momentum, normalized
by the fermion mass,

$\displaystyle-\frac{{\rm Im}[\Delta m_{q}(\mu)]}{{\rm Re}(m_{q})}$
$\displaystyle=-\frac{{\rm
Im}(y_{q}^{2})}{16\pi^{2}}F_{0}(0,m_{q}^{2},0)=\frac{\text{Im}\left(y_{q}^{2}\right)}{16\pi^{2}}\left(1-\log\frac{m_{q}^{2}}{\mu^{2}}\right)\,.$
(3.24)

In Fig. (3(b)), we take $m_{\phi}\lesssim m_{q}$. Here, $\mu=m_{q}$. Again,
the blue line represents the contribution to the QCD $\theta$ parameter from
the imaginary part of the colored fermion mass in Eq. (3.8) with the Fujikawa
method in Eq. (1.11), while the orange one comes from the direct two-loop
diagram calculation in Eq. (3.4). The blue and orange lines coincide when
$m_{\phi}\ll m_{q}$, while all three diagrams at two-loop level contribute to
the QCD $\theta$ parameter, as mentioned above.

Thus, when $m_{\phi}\gg m_{q}$ or $m_{\phi}\ll m_{q}$, the radiative
correction to the QCD $\theta$ parameter at two-loop level can be sufficiently
evaluated by the Fujikawa method taking into account the radiative correction
to the fermion masses in a limit of the zero external momentum. We consider
that the coincidence in $m_{\phi}\ll m_{q}$ is accidental.

### 3.2 Two hierarchical flavor case

Next, we consider a more complicated situation: there are two kinds of
fermions $q_{l}$ (light fermion, $m_{l}\ll m_{\phi}$) and $q_{h}$ (heavy one,
$m_{h}\approx m_{\phi}$), and the interactions are

$\displaystyle-\mathcal{L}=\bar{q}_{i}\left[\text{Re}\left(m_{i}\right)+\text{Im}\left(m_{i}\right)i\gamma_{5}\right]q_{i}+\left(\bar{q}_{i}y_{ij}P_{R}q_{j}\phi+\text{h.c.}\right)\,,$
(3.25)

where $i,j$ run $l$ (light) and $h$ (heavy), and we take a flavor diagonal
mass basis. We suppose that the Yukawa interactions $y_{lh}$ and $y_{hl}$ are
$CP$-violating complex couplings but $y_{ll}$ and $y_{hh}$ are real ones, for
simplicity.

From the direct loop calculations, the radiative corrections to the QCD
$\theta$ parameter are

$\displaystyle\delta\theta|_{\rm 1L}$
$\displaystyle=-\sum_{i}\frac{\text{Im}(m_{i})}{\text{Re}(m_{i})}$
$\displaystyle=-\frac{\text{Im}(m_{l})}{\text{Re}(m_{l})}-\frac{\text{Im}(m_{h})}{\text{Re}(m_{h})}\,,$
(3.26) $\displaystyle\delta\theta|_{\rm 2L}$
$\displaystyle=\sum_{i,j}\frac{{\rm
Im}\left(y_{ij}y_{ji}\right)}{16\pi^{2}}{\text{Re}(m_{i})}{\text{Re}(m_{j})}\left[I_{(2;2)}(m_{i}^{2};m_{j}^{2};m_{\phi}^{2})+2\bar{I}_{(3;1)}(m_{i}^{2};m_{j}^{2};m_{\phi}^{2})\right]$
$\displaystyle=\frac{{\rm
Im}\left(y_{lh}y_{hl}\right)}{8\pi^{2}}{\text{Re}(m_{l})}{\text{Re}(m_{h})}$
$\displaystyle\times\left[I_{(2;2)}(m_{l}^{2};m_{h}^{2};m_{\phi}^{2})+\bar{I}_{(3;1)}(m_{l}^{2};m_{h}^{2};m_{\phi}^{2})+\bar{I}_{(3;1)}(m_{h}^{2};m_{l}^{2};m_{\phi}^{2})\right]\,.$
(3.27)

Here, $I_{(2;2)}(x_{1};x_{2};x_{3})=I_{(2;2)}(x_{2};x_{1};x_{3})$ is used.
When $m_{l}\ll m_{h},\,m_{\phi}$, these loop functions can be simplified to

$\displaystyle{I}_{(2;2)}(m_{l}^{2};m_{h}^{2};m_{\phi}^{2})+\bar{I}_{(3;1)}(m_{l}^{2};m_{h}^{2};m_{\phi}^{2})+\bar{I}_{(3;1)}(m_{h}^{2};m_{l}^{2};m_{\phi}^{2})$
$\displaystyle\quad\quad=\frac{1}{2m_{l}^{2}}\left[1-\frac{1}{m_{h}^{2}-m_{\phi}^{2}}\left(m_{h}^{2}\log\frac{m_{h}^{2}}{{\mu^{2}}}-m_{\phi}^{2}\log\frac{m_{\phi}^{2}}{\mu^{2}}\right)\right]$
$\displaystyle\quad\qquad+\frac{1}{2m_{h}^{2}}\left[1-\frac{1}{m_{h}^{2}-m_{\phi}^{2}}\left(m_{h}^{2}\log\frac{m_{h}^{2}}{\mu^{2}}-m_{\phi}^{2}\log\frac{m_{\phi}^{2}}{\mu^{2}}\right)\right]$
$\displaystyle\quad\qquad-\frac{m_{h}^{4}-4m_{\phi}^{2}m_{h}^{2}+3m_{\phi}^{4}+2m_{\phi}^{4}\log\frac{m_{h}^{2}}{m_{\phi}^{2}}}{2(m_{h}^{2}-m_{\phi}^{2})^{3}}\log\frac{m_{l}^{2}}{m_{\phi}^{2}}$
$\displaystyle\quad\qquad+\frac{(m_{h}^{4}-3m_{h}^{2}m_{\phi}^{2}-2m_{\phi}^{4})\log\frac{m_{h}^{2}}{m_{\phi}^{2}}-4m_{\phi}^{4}\leavevmode\nobreak\
{\rm
Li}_{2}\left(1-\frac{m_{h}^{2}}{m_{\phi}^{2}}\right)}{2(m_{h}^{2}-m_{\phi}^{2})^{3}}\,.$
(3.28)

Let us compare the radiative corrections to the QCD $\theta$ parameter in the
effective theory approach. When the heavy fermion and the scalar boson are
integrated out from the full theory, the following effective interactions are
obtained,

$\displaystyle-\mathcal{L}_{\rm eff}$
$\displaystyle=\bar{q}_{l}\left[\text{Re}\left(m_{l}+\Delta
m_{l}\right)+\text{Im}\left(m_{l}+\Delta m_{l}\right)i\gamma_{5}\right]q_{l}$
$\displaystyle\quad+\frac{i}{2}g_{s}\tilde{d}_{q_{l}}\bar{q}_{l}\sigma^{\mu\nu}\gamma_{5}T^{a}q_{l}G_{\mu\nu}^{a}+\Delta\theta|_{\rm
th}\frac{\alpha_{s}}{8\pi}G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}\,.$ (3.29)

In addition to the Weinberg operator, the $CP$-violating four-Fermi operators
are not generated since it is assumed that only $y_{lh}$ and $y_{hl}$ have
complex phases. The one-loop correction $\text{Im}(\Delta m_{l})$, which is
the correction in a limit of the zero external momentum, is

$\displaystyle\text{Im}[\Delta m_{l}(\mu)]$ $\displaystyle=\frac{{\rm
Im}\left(y_{lh}y_{hl}\right)}{16\pi^{2}}{\text{Re}(m_{h})}F_{0}(0,m_{h}^{2},m_{\phi}^{2})$
$\displaystyle=-\frac{{\rm
Im}\left(y_{lh}y_{hl}\right)}{16\pi^{2}}{\text{Re}(m_{h})}\left[1-\frac{1}{m_{h}^{2}-m_{\phi}^{2}}\left(m_{h}^{2}\log\frac{m_{h}^{2}}{\mu^{2}}-m_{\phi}^{2}\log\frac{m_{\phi}^{2}}{\mu^{2}}\right)\right]\,.$
(3.30)

Note that the one-loop corrected mass $(m_{l}+\Delta m_{l})$ is
$\mu$-independent because $m_{l}$ is the $\overline{\rm MS}$ masses. In
addition, the CEDM for $q_{l}$ at one-loop level is given as

$\displaystyle\tilde{d}_{q_{l}}$ $\displaystyle=$ $\displaystyle-\frac{{\rm
Im}\left(y_{lh}y_{hl}\right)}{16\pi^{2}}\frac{{\text{Re}(m_{h})}}{m_{\phi}^{2}}f\left(\frac{m_{h}^{2}}{m_{\phi}^{2}}\right)\,,$
(3.31)

where

$\displaystyle f(x)$ $\displaystyle=$
$\displaystyle\frac{1}{2(1-x)^{3}}(3-4x+x^{2}+2\log x)\,.$ (3.32)

Using the above results, the radiative correction to the QCD $\theta$
parameter at two-loop level in Eq. (3.28) is reduced to

$\displaystyle\delta\theta$
$\displaystyle=\left.\delta\theta\right|_{1L}+\left.\delta\theta\right|_{2L}$
$\displaystyle=-\frac{\text{Im}\left[m_{l}+\Delta
m_{l}(\mu)\right]}{\text{Re}(m_{l})}-2\tilde{d}_{q_{l}}{\text{Re}(m_{l})}\log\frac{m_{l}^{2}}{m_{\phi}^{2}}+\left.\Delta\theta\right|_{\rm
th}\,.$ (3.33)

The first two terms come from the integration of $q_{l}$ in the effective
theory. The CEDM contribution is proportional to a large log,
$\log(m_{l}/m_{\phi})$. It is expected since the renormalization-group
equation of $\theta$ involves a term proportional to the CEDM in Eq. (2.18).
These results are easily improved by the renormalization-group equations. The
threshold correction $\Delta\theta|_{\rm th}$ represents the contributions
that come from the integration of $q_{h}$ and $m_{\phi}$ in the full theory,
and we obtain it from the above matching as,

$\displaystyle\left.\Delta\theta\right|_{\rm th}$
$\displaystyle=-\frac{\text{Im}\left[m_{h}+\Delta
m_{h}(\mu)\right]}{\text{Re}(m_{h})}$ $\displaystyle\quad+\frac{{\rm
Im}\left(y_{lh}y_{hl}\right)}{16\pi^{2}}{\text{Re}(m_{l})}{\text{Re}(m_{h})}\frac{\left(m^{2}_{h}m_{\phi}^{2}+3m^{4}_{\phi}\right)\log\frac{m^{2}_{h}}{m^{2}_{\phi}}+4m^{4}_{\phi}\leavevmode\nobreak\
{\rm
Li}_{2}\left(1-\frac{m^{2}_{h}}{m_{\phi}^{2}}\right)}{(m^{2}_{\phi}-m^{2}_{h})^{3}}\,.$
(3.34)

Here, $\text{Im}[\Delta m_{h}(\mu)]$ is the radiative correction to the
imaginary part of heavy fermion mass in a limit of the zero external momentum,

$\displaystyle\text{Im}[\Delta m_{h}(\mu)]$ $\displaystyle=-\frac{{\rm
Im}\left(y_{lh}y_{hl}\right)}{16\pi^{2}}{\text{Re}(m_{l})}\left(1-\log\frac{m_{\phi}^{2}}{\mu^{2}}\right)\,.$
(3.35)

The threshold correction to the QCD $\theta$ parameter is suppressed by
$\text{Re}(m_{l})/\text{Re}(m_{h})$, since we assume that $y_{lh}$ and
$y_{hl}$ are $CP$-violating complex couplings while $y_{ll}$ and $y_{hh}$ are
real ones. The CEDM contribution in the effective theory is also suppressed by
$\text{Re}(m_{l})\text{Re}(m_{h})/m_{\phi}^{2}$. Thus, the correction to the
QCD $\theta$ parameter can be evaluated by the Fujikawa method taking into
account the radiative correction to the imaginary part of the light quark mass
in Eq. (3.30) in this case.

## 4 Nelson-Barr Model

In this section, we show an example in which the QCD $\theta$ parameter is
induced at two-loop level while it vanishes at tree level. The $CP$ symmetry
is spontaneously broken in the Nelson-Barr models Nelson:1983zb ; Barr:1984qx
. The QCD $\theta$ term vanishes at the tree level by assuming quark mass
matrices ${\cal M}_{q}(q=u,d)$ with ${\rm arg\,det}\,({\cal M}_{q})=0$ even
after the spontaneous $CP$-violating sector.

The minimal models were derived by Bento, Branco, and Parada Bento:1991ez .
They are classified into two types of models, depending on whether up- or
down-type quarks are coupled with the $CP$-violating sectors. In this paper,
we consider the $d$-type model,#7#7#7 The $u$-type model is given by
$\displaystyle-{\cal L}^{u}_{Y}$ $\displaystyle=$ $\displaystyle
y_{u}^{ij}\tilde{H}\bar{Q}_{i}u_{j}+y_{d}^{ij}{H}\bar{Q}_{i}d_{j}+g^{aj}\Sigma_{a}\bar{\psi}u_{j}+m\bar{\psi}\psi^{c}+\textrm{h.c.}\,.$
given by

$\displaystyle-{\cal L}^{d}_{Y}$ $\displaystyle=$ $\displaystyle
y_{u}^{ij}{\tilde{H}}\bar{Q}_{i}u_{j}+y_{d}^{ij}{H}\bar{Q}_{i}d_{j}+g^{aj}\Sigma_{a}\bar{\psi}d_{j}+m\bar{\psi}\psi^{c}+\textrm{h.c.}\,,$
(4.1)

where the down-type quarks are coupled with the $CP$-violating sector. Here,
$Q$, $u$, and $d$ are for $SU(2)_{L}$ doublet left-handed quarks, singlet
right-handed up and down quarks, respectively. The subscripts $i,j,k$ are used
for the generation in the SM ($i,j,k=1,2,3$). We introduce an $SU(2)_{L}$
singlet vector-like quark, whose left and right-handed components are
represented by $\psi$ and $\psi^{c}$, respectively. The vacuum expectation
value of the SM Higgs doublet field $H$ is $\langle
H\rangle=(0,v/\sqrt{2})^{T}$. The conjugate is $\tilde{H}=\epsilon H^{\ast}$.
The complex singlet scalar fields $\Sigma_{a}$ ($a=1,\cdots,N_{\Sigma}$,
$N_{\Sigma}>1$) are assumed to get the complex vacuum expectation values so
that the $CP$ symmetry is spontaneously broken. The gauge charges of the
matter contents are summarised in Table 1. In this model, the down-type quark
mass matrix is given by

$\displaystyle\left(\begin{array}[]{cc}\bar{Q}&\bar{\psi}\end{array}\right){\cal
M}_{d}\left(\begin{array}[]{c}d\\\ \psi^{c}\end{array}\right)$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\bar{Q}&\bar{\psi}\end{array}\right)\left(\begin{array}[]{cc}y_{d}\langle{H}\rangle&0\\\
\xi^{\dagger}&m\end{array}\right)\left(\begin{array}[]{c}d\\\
\psi^{c}\end{array}\right)\,,$ (4.10)

where $\xi^{\ast}_{i}=g^{ai}\langle\Sigma_{a}\rangle$. We have to forbid terms
such as $\Sigma_{a}\bar{\psi}\psi^{c}$ and ${H}\bar{Q}\psi^{c}$ since they
lead to ${\rm arg\,det}({\cal M}_{d})\neq 0$. Those terms are forbidden by the
following discrete $\mathbb{Z}_{N}$ symmetry imposed as,

$\displaystyle\Sigma_{a}$ $\displaystyle\rightarrow$ $\displaystyle{\rm
e}^{\frac{2\pi ik}{N}}\Sigma_{a}\,,$ $\displaystyle\psi$
$\displaystyle\rightarrow$ $\displaystyle{\rm e}^{\frac{2\pi ik}{N}}\psi\,,$
(4.11) $\displaystyle\psi^{c}$ $\displaystyle\rightarrow$ $\displaystyle{\rm
e}^{\frac{2\pi ik}{N}}\psi^{c}\,.$

---
| chirality | $SU(3)_{C}$ | $SU(2)_{L}$ | $U(1)_{Y}$ | $\mathbb{Z}_{N}$
$Q$ | L | $\bf 3$ | $\bf 2$ | $1/6$ | $0$
$u$ | R | $\bf 3$ | $\bf 1$ | $2/3$ | $0$
$d$ | R | $\bf 3$ | $\bf 1$ | $-1/3$ | $0$
$H$ | – | $\bf 1$ | $\bf 2$ | $1/2$ | $0$
$\psi$ | L | $\bf 3$ | $\bf 1$ | $-1/3$ | $k$
$\psi^{c}$ | R | $\bf 3$ | $\bf 1$ | $-1/3$ | $k$
$\Sigma_{a}$ | – | $\bf 1$ | $\bf 1$ | $0$ | $k$
$S$ | – | $\bf 1$ | $\bf 1$ | $0$ | $0$
Table 1: The matter contents and their gauge charges in a Nelson-Barr type
(spontaneous $CP$ violating) model.

It is pointed out by Ref. Dine:2015jga  that the four-point scalar
interactions of $H$ and $\Sigma_{a}$ lead to ${\rm arg\,det}({\cal M}_{d})\neq
0$ at one-loop level and the QCD $\theta$ parameter is radiatively
generated.#8#8#8 When the four-point scalar interactions are absent, the QCD
$\theta$ parameter is still generated at three-loop level. It is evaluated
using the dimensional analysis in Refs. Nelson:1984hg ; Valenti:2021rdu . It
corresponds to the two-loop correction to the QCD $\theta$ parameter, shown in
the left diagram of Fig. 4.

Figure 4:  The two-loop diagrams that contribute to the radiative QCD theta
parameter in the Nelson-Barr model. The external gluon fields are omitted. The
cross ($\times$) represents the chirality flip. Subscripts $A$ and $B$ are for
mass eigenstates of fermions.

Furthermore, we introduce a real scalar field $S$ whose real vacuum
expectation value leads to the singlet vector-like quark mass, $m$, as
$m=f\langle S\rangle$ with a Yukawa coupling $-{\cal
L}=fS\bar{\psi}\psi^{c}+{\rm h.c.}$. From a viewpoint of model-building,
however, the origin of $m$ should be related to the vacuum expectation values
of $\Sigma_{a}$ since the observed CKM $CP$ phase of $O(1)$ is realized only
when $m\sim|\xi_{i}|=|g^{ai}\langle\Sigma_{a}\rangle|$. It is found that the
four-point scalar coupling of $S$ and $\Sigma_{a}$ also generates the QCD
$\theta$ parameter at two-loop level, denoted as a right diagram of Fig. 4.
Then, we consider the following four-point scalar interactions,

$\displaystyle V(\Sigma_{a},H)$ $\displaystyle=$
$\displaystyle\gamma_{ab}\Sigma_{a}^{\ast}\Sigma_{b}|H|^{2}+\frac{1}{2}\tilde{\gamma}_{ab}\Sigma_{a}^{\ast}\Sigma_{b}S^{2}\cdots\,.$
(4.12)

We evaluate the QCD $\theta$ parameter at two-loop level assuming
$\gamma_{ab}$ and $\tilde{\gamma}_{ab}$ are non-vanishing.

The down-type quark mass matrix is diagonalized as

$\displaystyle{\cal M}_{d}$ $\displaystyle=$ $\displaystyle O_{d}\bar{{\cal
M}}_{d}U^{\dagger}_{d}\,,$ (4.13)

where $\bar{{\cal M}}_{d}$ is the diagonal mass eigenvalue matrix $\bar{{\cal
M}}_{d}={\rm diag}(\bar{m}_{d_{A}})$ ($A=1,\cdots,4)$. The inverse of ${\cal
M}_{d}$ is given by the unitary matrices as

$\displaystyle{\cal M}^{-1}_{d}$ $\displaystyle=$ $\displaystyle
U_{d}\bar{{\cal M}}^{-1}_{d}O_{d}^{\dagger}$ (4.16) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}y_{d}^{-1}\langle H\rangle^{-1}&0\\\
-\frac{1}{m}\xi^{\dagger}y_{d}^{-1}\langle
H\rangle^{-1}&\frac{1}{m}\end{array}\right)\,.$

We assume that $\langle H\rangle\ll m$, $|\xi_{k}|$. In this case, the unitary
matrices $U_{d}$ and $O_{d}$ are given as#9#9#9 This diagonalization comes
from a mathematical identity Valenti:2021rdu ,  $\displaystyle{\cal
M}_{d}\tilde{U}_{d}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}Y_{d}\langle{H}\rangle&Y\langle{H}\rangle\\\
0&M\end{array}\right)\,.$ (4.19)

$\displaystyle U_{d}=\tilde{U}_{d}\left(\begin{array}[]{cc}u_{d}&0\\\
0&1\end{array}\right)\,,$ $\displaystyle
O_{d}=\tilde{O}_{d}\left(\begin{array}[]{cc}o_{d}&0\\\
0&1\end{array}\right)\,.$ (4.24)

Here,

$\displaystyle\tilde{U}_{d}=\left(\begin{array}[]{cc}1-\frac{\xi\xi^{\dagger}}{|\xi|^{2}}(1-\frac{m}{M})&\frac{\xi}{M}\\\
-\frac{\xi^{\dagger}}{M}&\frac{m}{M}\end{array}\right)\,,$
$\displaystyle\tilde{O}_{d}=\left(\begin{array}[]{cc}1&\frac{Y\langle
H\rangle}{M}\\\ -\frac{Y^{\dagger}\langle H\rangle}{M}&1\end{array}\right)\,,$
(4.29)

where $M^{2}=m^{2}+\xi^{\dagger}\xi$ and
$Y^{i}=y_{d}^{ij}\frac{\xi_{j}}{M}=Y_{d}^{ij}\frac{\xi_{j}}{m}$. The mass
eigenvalues of the ${\cal M}_{d}$ are given by $\bar{m}_{d_{4}}=M$ and
eigenvalues of $Y_{d}^{ij}\langle H\rangle$ for three SM down quarks,

$\displaystyle Y_{d}^{ij}$ $\displaystyle=$ $\displaystyle
y_{d}^{ik}\left[1-\frac{\xi\xi^{\dagger}}{|\xi|^{2}}x\left(1-\frac{m}{M}\right)\right]_{kj}.\,$
(4.30)

Two $(3\times 3)$ unitary matrices in Eq. (4.24), $u_{d}$ and $o_{d}$, are for
diagonalization of $Y_{d}$.

Now we evaluate the corrections to the QCD $\theta$ parameter at two-loop
level, as shown in Fig. 4. They come from mixings between $H$ and $\Sigma_{a}$
and also between $S$ and $\Sigma_{a}$. Assuming $\gamma_{ab}$ and
$\tilde{\gamma}_{ab}$ much less than one, we keep only the leading terms in
the perturbation. In the case of diagrams in Fig. 4, they give the corrections
to the QCD $\theta$ parameter as

$\displaystyle\Delta\theta_{\Sigma}$ $\displaystyle=$
$\displaystyle\frac{1}{8\pi^{2}}y_{d}^{ij}g^{ak}\gamma_{ab}{\rm
Im}\left[[U_{d}]_{kA}\bar{m}_{d_{A}}[O^{\dagger}_{d}]_{Ai}[U_{d}]_{jB}\bar{m}_{d_{B}}[O^{\dagger}_{d}]_{B4}\right]$
(4.31) $\displaystyle\times\frac{\langle
H\rangle\langle\Sigma_{b}\rangle}{m_{h}^{2}-m_{\Sigma_{a}}^{2}}\left[\left(I_{(2;2)}+\bar{I}_{(1;3)}+\bar{I}_{(3;1)}\right)(\bar{m}_{d_{A}}^{2};\bar{m}_{d_{B}}^{2};m_{h}^{2})-(m_{h}^{2}\rightarrow
m_{\Sigma_{a}}^{2})\right]\,,$ $\displaystyle\Delta\theta_{S}$
$\displaystyle=$
$\displaystyle\frac{1}{8\pi^{2}}fg^{ak}\tilde{\gamma}_{ab}{\rm
Im}\left[[U_{d}]_{kA}\bar{m}_{d_{A}}[O^{\dagger}_{d}]_{A4}[U_{d}]_{4B}\bar{m}_{d_{B}}[O^{\dagger}_{d}]_{B4}\right]$
(4.32) $\displaystyle\times\frac{\langle
S\rangle\langle\Sigma_{b}\rangle}{m_{S}^{2}-m_{\Sigma_{a}}^{2}}\left[\left(I_{(2;2)}+\bar{I}_{(1;3)}+\bar{I}_{(3;1)}\right)(\bar{m}_{d_{A}}^{2};\bar{m}_{d_{B}}^{2};m_{S}^{2})-(m_{S}^{2}\rightarrow
m_{\Sigma_{a}}^{2})\right]\,,$

where $m_{h}$, $m_{S}$, and $m_{\Sigma_{a}}$ are masses of $H$, $S$, and
$\Sigma_{a}$, respectively. These results still include fake IR bad behavior
due to light quarks such as $1/\bar{m}_{d_{k}}$. However, using
$\bar{I}_{(3;1)}(x_{1},x_{2},x_{3})\simeq-1/(2x_{1})F_{0}(0,x_{2},x_{3})$ for
$x_{1}\ll x_{3}$ and inverse matrix of ${\cal M}_{d}$ in Eq. (4.16), the IR
behavior can be removed. Then, when $\langle H\rangle\ll m$, $|\xi_{k}|$, the
above results are reduced as

$\displaystyle\Delta\theta_{\Sigma}$ $\displaystyle=$
$\displaystyle-\frac{1}{16\pi^{2}}\gamma_{ab}g^{ak}g^{ck}{\rm
Im}[\langle\Sigma_{b}\rangle\langle\Sigma_{c}\rangle^{\ast}]$ (4.33)
$\displaystyle\times\frac{1}{m_{h}^{2}-m_{\Sigma_{a}}^{2}}\left[F_{0}(0;M^{2};m_{h}^{2})-(m_{h}^{2}\rightarrow
m_{\Sigma_{a}}^{2})\right]\,,$ $\displaystyle\Delta\theta_{S}$
$\displaystyle=$
$\displaystyle\frac{1}{8\pi^{2}}\tilde{\gamma}_{ab}g^{ak}g^{ck}{\rm
Im}[\langle\Sigma_{b}\rangle\langle\Sigma_{c}\rangle^{\ast}]\textbf{}$ (4.34)
$\displaystyle\times\frac{Mf\langle
S\rangle}{m_{S}^{2}-m_{\Sigma_{a}}^{2}}\left[\left(I_{(2;2)}+2\bar{I}_{(3;1)}\right)(M^{2};M^{2};m_{S}^{2})-(m_{S}^{2}\rightarrow
m_{\Sigma_{a}}^{2})\right]\,.$

Since the $CP$ symmetry is spontaneously broken, the corrections to the QCD
$\theta$ parameter are UV finite. It is found that $\Delta\theta_{\Sigma}$
comes mainly from the correction to light quark masses since it is given by
the function $F_{0}$. On the other hand, $\Delta\theta_{S}$ is generated by
the vector-like quark loop. When $M$ is comparable to $m_{S}$ and/or
$m_{\Sigma_{a}}$, the correction is not dominated by the correction to the
vector-like quark mass, and the full evaluation of the two-loop diagram is
required to evaluate them.

## 5 Conclusion and Discussion

The QCD $\theta$ parameter is generated radiatively in spontaneous $CP$ or $P$
symmetry-breaking models, which solve the strong $CP$ problem. We scrutinized
the QCD $\theta$ parameter at the two-loop level analysis. In the simplified
models with $CP$-violating Yukawa interactions, we observed that the two-loop
calculation of the radiative QCD $\theta$ parameter using the Fock-Schwinger
gauge method is consistent with the effective field theory approach at the
low-energy scale. Furthermore, we clarified the application scope of the
Fujikawa method. When there is a scale hierarchy in the particle masses in
$CP$-violating sector, the Fujikawa method is sufficient for evaluating the
QCD $\theta$ parameter. On the other hand, in the case of a small hierarchy,
the Fujikawa method does not work well. The Nelson-Barr model is an example
that the Fujikawa method cannot evaluate the radiatively generated QCD
$\theta$ parameter correctly. If the vector-like quark and additional scalars
have comparable masses, the Fock-Schwinger gauge method should be used to
evaluate the radiative QCD $\theta$ parameter.

It is an important subject to evaluate the three-loop contributions to the QCD
$\theta$ parameter in some well-motivated models for the strong $CP$ problem,
such as the left-right models, and to compare the predictions with the
experimental bounds on the QCD $\theta$ parameter. In the evaluation, the
effective theory approach would be useful if the CP-violating interaction can
be integrated out.

###### Acknowledgements.

This work is supported by the JSPS Grant-in-Aid for Scientific Research Grant
No. 20H01895 (J.H.) and No. 21K03572 (J.H.). 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 (N.O.) 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}^{2}-2x_{1}x_{2}-2x_{2}x_{3}-2x_{3}x_{1}}$
and $Q^{2}\equiv 4\pi\mu^{2}{\rm e}^{-\gamma_{E}}$ Ford:1992pn ;
Espinosa:2000df ; Martin:2001vx . Note that although the last terms of
$\bar{I}_{\epsilon}$ are UV finite, they do not affect any physical quantity
Martin:2018emo , whose terms are suppressed by $\mathcal{O}(\epsilon)$ at one-
loop level, while they are uplifted to $\mathcal{O}(\epsilon^{0})$ at two-loop
level.

The radiative correction to the QCD $\theta$ parameter at two-loop level is
proportional to
${I}_{(2;2)}(x_{1};x_{2};x_{3})+\bar{I}_{(3;1)}(x_{1};x_{2};x_{3})+\bar{I}_{(3;1)}(x_{2};x_{1};x_{3})$
($\bar{I}_{(3;1)}$ is defined as Eq. (3.5)). We provide it with some limits
here. When $x_{1}\ll x_{2},\leavevmode\nobreak\ x_{3}$, we have

$\displaystyle{I}_{(2;2)}(x_{1};x_{2};x_{3})+\bar{I}_{(3;1)}(x_{1};x_{2};x_{3})+\bar{I}_{(3;1)}(x_{2};x_{1};x_{3})$
$\displaystyle\quad\quad=\frac{1}{2x_{1}}\left[1-\frac{1}{x_{3}-x_{2}}\left(x_{3}\log\frac{x_{3}}{\mu^{2}}-x_{2}\log\frac{x_{2}}{{\mu^{2}}}\right)\right]$
$\displaystyle\quad\qquad+\frac{1}{2x_{2}}\left[1-\frac{1}{x_{3}-x_{2}}\left(x_{3}\log\frac{x_{3}}{\mu^{2}}-x_{2}\log\frac{x_{2}}{\mu^{2}}\right)\right]$
$\displaystyle\quad\qquad+\frac{3x_{3}^{2}-4x_{3}x_{2}+x_{2}^{2}+2x_{3}^{2}\log\frac{x_{2}}{x_{3}}}{2(x_{3}-x_{2})^{3}}\log\frac{x_{1}}{x_{3}}$
$\displaystyle\quad\qquad-\frac{({x_{2}^{2}-3x_{2}x_{3}-2x_{3}^{2}})\log\frac{x_{2}}{x_{3}}-4x_{3}^{2}\leavevmode\nobreak\
{\rm Li}_{2}\left(1-\frac{x_{2}}{x_{3}}\right)}{2(x_{3}-x_{2})^{3}}\,.$ (A.9)

To obtain this analytic formula, we expand $R$ by $x_{1}$,

$\displaystyle
R\simeq\pm\left[(x_{3}-x_{2})-\frac{x_{2}+x_{3}}{x_{3}-x_{2}}x_{1}-\frac{2x_{2}x_{3}}{(x_{3}-x_{2})^{3}}x_{1}^{2}-\frac{2x_{2}x_{3}(x_{2}+x_{3})}{(x_{3}-x_{2})^{5}}x_{1}^{3}\right]\,,$
(A.10)

where the overall $+$ ($-$) sign is for $x_{3}>x_{2}$ ($x_{2}>x_{3}$).

When $x_{1},\leavevmode\nobreak\ x_{2}\ll x_{3}$, we obtain

$\displaystyle{I}_{(2;2)}(x_{1};x_{2};x_{3})+\bar{I}_{(3;1)}(x_{1};x_{2};x_{3})+\bar{I}_{(3;1)}(x_{2};x_{1};x_{3})$
$\displaystyle\quad\quad=\frac{1}{2x_{1}}\left(1-\log\frac{x_{3}}{\mu^{2}}{+\frac{x_{2}}{x_{3}}\log\frac{x_{2}}{x_{3}}}\right)+\frac{1}{2x_{2}}\left(1-\log\frac{x_{3}}{\mu^{2}}{+\frac{x_{1}}{x_{3}}\log\frac{x_{1}}{x_{3}}}\right)$
$\displaystyle\quad\qquad+\frac{2\pi^{2}+9(\log\frac{x_{1}}{x_{3}}+\log\frac{x_{2}}{x_{3}})+6\log\frac{x_{1}}{x_{3}}\log\frac{x_{2}}{x_{3}}}{6x_{3}}\,,$
(A.11)

and

$\displaystyle{I}_{(2;2)}(x_{1};x_{1};x_{3})+2\bar{I}_{(3;1)}(x_{1};x_{1};x_{3})$
$\displaystyle\quad\quad=\frac{1}{x_{1}}\left(1-\log\frac{x_{3}}{\mu^{2}}\right)+\frac{\pi^{2}+{12}\log\frac{x_{1}}{x_{3}}+3\log^{2}\frac{x_{1}}{x_{3}}}{3x_{3}}\,.$
(A.12)

## References

  * (1) C. Abel et al., “Measurement of the Permanent Electric Dipole Moment of the Neutron,” Phys. Rev. Lett. 124 (2020) 081803 [arXiv:2001.11966].
  * (2) J. Liang, et al., “Nucleon Electric Dipole Moment from the $\theta$ Term with Lattice Chiral Fermions.” arXiv:2301.04331.
  * (3) ACME Collaboration, “Improved limit on the electric dipole moment of the electron,” Nature 562 (2018) 355–360.
  * (4) T. S. Roussy et al., “An improved bound on the electron’s electric dipole moment,” Science 381 (2023) adg4084 [arXiv:2212.11841].
  * (5) V. V. Flambaum, M. Pospelov, A. Ritz, and Y. V. Stadnik, “Sensitivity of EDM experiments in paramagnetic atoms and molecules to hadronic CP violation,” Phys. Rev. D 102 (2020) 035001 [arXiv:1912.13129].
  * (6) R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38 (1977) 1440–1443.
  * (7) S. Weinberg, “A New Light Boson?” Phys. Rev. Lett. 40 (1978) 223–226.
  * (8) F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978) 279–282.
  * (9) M. Kamionkowski and J. March-Russell, “Planck scale physics and the Peccei-Quinn mechanism,” Phys. Lett. B 282 (1992) 137–141 [hep-th/9202003].
  * (10) 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].
  * (11) S. M. Barr and D. Seckel, “Planck scale corrections to axion models,” Phys. Rev. D 46 (1992) 539–549.
  * (12) A. E. Nelson, “Naturally Weak CP Violation,” Phys. Lett. B 136 (1984) 387–391.
  * (13) S. M. Barr, “Solving the Strong CP Problem Without the Peccei-Quinn Symmetry,” Phys. Rev. Lett. 53 (1984) 329.
  * (14) S. M. Barr, “A Natural Class of Nonpeccei-quinn Models,” Phys. Rev. D 30 (1984) 1805.
  * (15) M. A. B. Beg and H. S. Tsao, “Strong P, T Noninvariances in a Superweak Theory,” Phys. Rev. Lett. 41 (1978) 278.
  * (16) R. N. Mohapatra and G. Senjanovic, “Natural Suppression of Strong p and t Noninvariance,” Phys. Lett. B 79 (1978) 283–286.
  * (17) K. S. Babu and R. N. Mohapatra, “A Solution to the Strong CP Problem Without an Axion,” Phys. Rev. D 41 (1990) 1286.
  * (18) S. M. Barr, D. Chang, and G. Senjanovic, “Strong CP problem and parity,” Phys. Rev. Lett. 67 (1991) 2765–2768.
  * (19) N. Craig, I. Garcia Garcia, G. Koszegi, and A. McCune, “P not PQ,” JHEP 09 (2021) 130 [arXiv:2012.13416].
  * (20) L. J. Hall and K. Harigaya, “Implications of Higgs Discovery for the Strong CP Problem and Unification,” JHEP 10 (2018) 130 [arXiv:1803.08119].
  * (21) J. de Vries, P. Draper, and H. H. Patel, “Do Minimal Parity Solutions to the Strong $CP$ Problem Work?” arXiv:2109.01630.
  * (22) J. R. Ellis and M. K. Gaillard, “Strong and Weak CP Violation,” Nucl. Phys. B 150 (1979) 141–162.
  * (23) K. Fujikawa, “Path Integral Measure for Gauge Invariant Fermion Theories,” Phys. Rev. Lett. 42 (1979) 1195–1198.
  * (24) S. L. Adler, “Axial vector vertex in spinor electrodynamics,” Phys. Rev. 177 (1969) 2426–2438.
  * (25) 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.
  * (26) 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.
  * (27) G. ’t Hooft and M. J. G. Veltman, “Regularization and Renormalization of Gauge Fields,” Nucl. Phys. B 44 (1972) 189–213.
  * (28) J. Hisano, T. Kitahara, N. Osamura, and A. Yamada, “Novel loop-diagrammatic approach to QCD $\theta$ parameter and application to the left-right model,” JHEP 03 (2023) 150 [arXiv:2301.13405].
  * (29) V. Fock, “Proper time in classical and quantum mechanics,” Phys. Z. Sowjetunion 12 (1937) 404–425.
  * (30) J. S. Schwinger, “On gauge invariance and vacuum polarization,” Phys. Rev. 82 (1951) 664–679.
  * (31) A. Schwarz, V. Fateev, and Y. Tyupkin, “On the particle-like solutions in the presence of fermions,” 155, Lebedev Institute, 1976.
  * (32) C. Cronstrom, “A SIMPLE AND COMPLETE LORENTZ COVARIANT GAUGE CONDITION,” Phys. Lett. B 90 (1980) 267–269.
  * (33) M. A. Shifman, “Wilson Loop in Vacuum Fields,” Nucl. Phys. B 173 (1980) 13–31.
  * (34) M. S. Dubovikov and A. V. Smilga, “Analytical Properties of the Quark Polarization Operator in an External Selfdual Field,” Nucl. Phys. B 185 (1981) 109–132.
  * (35) 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.
  * (36) 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)].
  * (37) S. N. Nikolaev and A. V. Radyushkin, “Method for Computing Higher Gluonic Power Corrections to QCD Charmonium Sum Rules,” Phys. Lett. B 110 (1982) 476. [Erratum: Phys.Lett.B 116, 469 (1982)].
  * (38) S. N. Nikolaev and A. V. Radyushkin, “Vacuum Corrections to QCD Charmonium Sum Rules: Basic Formalism and O ($G^{3}$) Results,” Nucl. Phys. B 213 (1983) 285–304.
  * (39) J. Hisano, K. Tsumura, and M. J. S. Yang, “QCD Corrections to Neutron Electric Dipole Moment from Dimension-six Four-Quark Operators,” Phys. Lett. B 713 (2012) 473–480 [arXiv:1205.2212].
  * (40) E. E. Jenkins, A. V. Manohar, and P. Stoffer, “Low-Energy Effective Field Theory below the Electroweak Scale: Anomalous Dimensions,” JHEP 01 (2018) 084 [arXiv:1711.05270].
  * (41) G. Boyd, A. K. Gupta, S. P. Trivedi, and M. B. Wise, “Effective Hamiltonian for the Electric Dipole Moment of the Neutron,” Phys. Lett. B 241 (1990) 584–588.
  * (42) E. Braaten, C.-S. Li, and T.-C. Yuan, “The Evolution of Weinberg’s Gluonic CP Violation Operator,” Phys. Rev. Lett. 64 (1990) 1709.
  * (43) D. Chang, W.-Y. Keung, C. S. Li, and T. C. Yuan, “QCD Corrections to CP Violation From Color Electric Dipole Moment of $b$ Quark,” Phys. Lett. B 241 (1990) 589–592.
  * (44) M. Dine and W. Fischler, “Constraints on New Physics From Weinberg’s Analysis of the Neutron Electric Dipole Moment,” Phys. Lett. B 242 (1990) 239–244.
  * (45) B. Henning, X. Lu, and H. Murayama, “How to use the Standard Model effective field theory,” JHEP 01 (2016) 023 [arXiv:1412.1837].
  * (46) B. Henning, X. Lu, and H. Murayama, “One-loop Matching and Running with Covariant Derivative Expansion,” JHEP 01 (2018) 123 [arXiv:1604.01019].
  * (47) A. I. Davydychev and J. B. Tausk, “Two loop selfenergy diagrams with different masses and the momentum expansion,” Nucl. Phys. B 397 (1993) 123–142.
  * (48) L. Bento, G. C. Branco, and P. A. Parada, “A Minimal model with natural suppression of strong CP violation,” Phys. Lett. B 267 (1991) 95–99.
  * (49) M. Dine and P. Draper, “Challenges for the Nelson-Barr Mechanism,” JHEP 08 (2015) 132 [arXiv:1506.05433].
  * (50) A. E. Nelson, “Calculation of $\theta$ Barr,” Phys. Lett. B 143 (1984) 165–170.
  * (51) A. Valenti and L. Vecchi, “The CKM phase and $\overline{\theta}$ in Nelson-Barr models,” JHEP 07 (2021) 203 [arXiv:2105.09122].
  * (52) 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)].
  * (53) 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].
  * (54) 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].
  * (55) S. P. Martin and H. H. Patel, “Two-loop effective potential for generalized gauge fixing,” Phys. Rev. D 98 (2018) 076008 [arXiv:1808.07615].