Model independent evaluation of the Wilson coefficient of the Weinberg operator in QCD
Tomohiro Abe Junji Hisano Ryo Nagai
hep-ph





  Institute for Advanced Research, Nagoya University,
  Furo-cho Chikusa-ku, Nagoya, Aichi, 464-8602 Japan


  Kobayashi-Maskawa Institute for the Origin of Particles and the
  Universe, Nagoya University,
  Furo-cho Chikusa-ku, Nagoya, Aichi, 464-8602 Japan




  Kobayashi-Maskawa Institute for the Origin of Particles and the
  Universe, Nagoya University,
  Furo-cho Chikusa-ku, Nagoya, Aichi, 464-8602 Japan


Department of Physics, Nagoya University,
Furo-cho Chikusa-ku, Nagoya, Aichi, 464-8602 Japan


  Kavli IPMU (WPI), UTIAS, University of Tokyo, Kashiwa, Chiba 277-8584, Japan



Department of Physics, Tohoku University, Sendai, Miyagi 980-8578, Japan


We derive a Wilson coefficient of a CP-violating purely gluonic dimension-6 operator 
called the Weinberg operator (GGG̃) 
generated by a scalar and two fermions at the two-loop level. 
We do not specify the representation of SU(3)_c for the scalar and the fermions,
and thus our result can be applied to a variety of models beyond the standard model. 
We estimate the nucleon EDMs induced by the Weinberg operator in some examples and discuss the importance of measuring EDMs. It is found that future measurements of the EDMs can probe physics at higher energy scale beyond the reach of collider experiments.




Model independent evaluation of the Wilson coefficient 
 of the Weinberg operator in QCD
    Ryo Nagai
    December 30, 2023
========================================================================================







§ INTRODUCTION

Measurements of electric dipole moments (EDMs) are very powerful for exploring physics beyond the standard model (SM).
The SM predicts small values of EDMs.
Its prediction for the neutron EDM is 
d_n ≃ 10^-32 e cm <cit.>,
and for the electron EDM is d_e ≤ 10^-38 e cm <cit.>. 
It is much smaller than the current upper bounds,
|d_n| < 2.9 × 10^-26 e cm (90% CL) <cit.>,

and |d_e| < 8.7 × 10^-29 e cm (90% CL) <cit.>.



On the other hand, models beyond the SM often have new CP violation sources,
and they can predict larger values of EDMs compared to the SM.
Therefore observation of EDMs is equivalent to a discovery of physics beyond the SM.
Moreover, EDMs can probe physics at higher energy scale beyond the reach of collider experiments.
The current data of the Large Hadron Collider (LHC) experiment imply that 
the scale of new physics is higher than O(1) TeV. 
Thus the importance of measuring EDMs is increasing.
This situation motivates us to evaluate EDMs in a variety of models with small theoretical uncertainties.




Weinberg pointed out that 
a pure gluonic dimension-6 operator, GGG̃, is a source of CP violation in QCD, 
and it affects to the EDM of neutrons <cit.>.
We call the operator as the Weinberg operator.
In renormalizable models, two-loop diagrams are at the leading order for the Weinberg operator.
Full two-loop calculations have done for some models:
two-Higgs doublet models with CP violation <cit.>,
the quark-squark-gluino diagram in MSSM <cit.>,
and
Left-Right symmetric models <cit.>.
Besides full two-loop calculations, it has been shown that the Weinberg operator is induced 
as a result of integrating out heavy colored particles which carry chromo-EDMs (cEDMs) <cit.>. 
The calculation is based on the effective theory, and the calculation is simpler than the full two-loop calculation.
However, it is justified only if the following two conditions are satisfied;
(i) One of the particles in the loop is lighter than the others and
(ii) the lightest particle is color non-singlet. 
It is useful to derive more general formulae to evaluate contributions to the Weinberg operator
from new physics in general setup.


In this paper, we calculate two-loop diagrams that contain one scalar and two fermions
without specifying their color representations,
and derive a general formula for the Wilson coefficient of the Weinberg operator in the system.
Our result can be applied to a variety models.
We also compare our full two-loop calculation to the result derived in the effective theory.
It is shown that the calculation based on the effective theory can overestimate the contribution
if the mass difference among the particles in the loop is small.
For the application of our general formula, we estimate the nucleon EDMs in some examples. 
Comparing our results with the latest upper bound on nucleon EDMs, 
we survey the parameter region which is consistent with current experimental results. 
We also discuss the impact of measuring EDMs by showing the parameter region which can be covered by the future observations of the EDMs. 




The rest of this paper is organized as follows.
At first, we introduce our setup and fix the notation for calculating the Wilson coefficient of the Weinberg operator in section <ref>. In section <ref>, we perform full two-loop calculation and derive a general formula for the Wilson coefficient of the Weinberg operator. 

The result of the two-loop calculation is summarized in section <ref>. In section <ref>, we discuss the relationship between our result and the result based on the effective theoretical approach. We estimate the nucleon EDMs in some examples and discuss the importance of measuring EDMs in section <ref>. Finally, we summarize our discussion in section <ref>. 




§ SETUP


We consider two fermions (A and B) and a scalar field (S) 
which have the following interactions,

    L⊃    
    - ψ̅_B g_B̅ASψ_A S
    - ψ̅_A g_A̅BS̅ψ_B S^*,




where




    g_B̅AS =    X_B̅AS (s + γ_5 a), 
    
     g_A̅BS̅ =    X^†_A̅BS̅ (s^* - γ_5 a^*).

All the indices of SU(3)_c representations are described by X.
Explicit expressions for some examples are shown in Table. <ref>.
s and a are complex numbers.
The Feynman rules are shown in Fig. <ref>.













In this setup, we will calculated the Wilson coefficient of the Weinberg operator, C_G, defined

    L⊃ g_s^3 C_G  O,




where




    O 
    =   
     -1/3 f^abc G^a_μν G^b^ν_ ρG̃^c ρμ 
    ,




and g_s is the SU(3)_c gauge coupling defined in the covariant derivative by




    D_μ = ∂_μ + i g_s T^a A^a_μ.




In Eq. (<ref>), 
f^abc is the structure constant of SU(3)_c, 
and G̃^c ρμ = 1/2ϵ^ρμαβ G^c_αβ 
with ϵ^0123 = +1.







§ EVALUATION OF DIAGRAMS


The leading order contributions to the Weinberg operator are given by two-loop diagrams.
We use the Fock-Schwinger gauge for the external gluon fields.
The technical details of the Fock-Schwinger gauge is reviewed, for example, in Refs. <cit.>.


In addition to the Weinberg operator (GGG̃), 
we also need to calculate operators which consist of G_μν, G̃_μν,
and two covariant derivatives, such as DGDG̃.

This is because 
the external gluon fields are described by covariant derivatives and gluon field strength in
the Fock-Schwinger gauge,
and commutation of two covariant derivatives is equivalent to gluon field strength.






 §.§ GGG̃ from DGDG̃ terms


We calculate GGG̃ from DGDG̃ terms. We evaluate diagrams with two external gluon fields.
Diagrams should have at least one gluon field from fermion line, otherwise terms with γ_5 all vanish.
In addition, diagrams with one gluon field from the scalar field also vanish.
See appendix <ref> for the detail.
Figure <ref> shows the diagrams we have to evaluate.
We call each diagram from left to right as (2,0), (1,1), and (0,2), respectively.






In the calculation, we obtain terms that are combination of X and SU(3)_c generator, such as
(X_B̅AS (T^a)_AA' (T^b)_A'A” X^†_A̅”BS̅).
Here we denote the SU(3)_c indices of T^a and T^b symbolically. 
Since X consists of ϵ tensors and Kronecker deltas, 
(X_B̅AS (T^a)_AA' (T^b)_A'A” X^†_A̅”BS̅) is proportional to (T^a T^b) ∝δ^ab.
Other similar terms are also proportional to δ^ab.
Thus we can introduce the following notations.

    (X_B̅AS (T^a)_AA' (T^b)_A'A” X^†_A̅”BS̅)
    ≡   (X T_A T_A X^†) δ^ab
    ,
    (X_B̅AS (T^a)_AA' X^†_A̅'B'S̅ (T^b)_B'B)
    ≡   (X T_A X^† T_B ) δ^ab
    ,
    (X_B̅AS X^†_A̅B'S̅ (T^a)_B'B” (T^b)_B”B)
    ≡   (X X^† T_B T_B ) δ^ab
    .

where the subscript F=A,B denotes that T_F is the generator for the fermion field F.
The values of (X T_A T_A X^†) etc. depend on the representations of A, B, and S.
We show some examples in Table <ref>.
We use  <cit.> to evaluate traces with γ-matrices.
We also use the following identities that can be proved by using Jacobi identities and Bianchi identities.




    D_α_1 G^a_αμ D_β_1 G^a_βν(
     g^αβϵ^α_1 β_1 μν
    +g^αβ_1ϵ^α_1 βμν
    +g^α_1 βϵ^αβ_1 μν
    +g^α_1 β_1ϵ^α_1 β_1 μν)
    =   
     27 g_s  O
    ,
    
    G^a_αμ D_β_2  D_β_1 G^a_βν(
     g^ββ_1ϵ^αβ_2 μν
    +g^ββ_2ϵ^αβ_1 μν
    +g^β_1 β_2ϵ^αβμν)
    =   
    - 18 g_s  O
    ,
    
    (D_α_2 D_α_1 G^a_αμ) G^a_βν(
     g^αα_1ϵ^α_2 βμν
    +g^αα_2ϵ^α_1 βμν
    +g^α_1 α_2ϵ^αβμν)
    =   
    - 18 g_s  O
    .




Since the subdiagrams of the (2,0) and (0,2) contain UV-divergences,
we need counterterms to renormalize them during our calculation.
However, as we will see in later, 
the final result after summing all of diagrams is independent from the UV-divergence.


We find the following DGDG̃ terms.

    (2,0)-diagram
    =   
    i 6/(4π)^4 g_s^3 (X T_A T_A X^†)  O
       ×[
    i Re(sa^*)/2 
    h_1(m_A^2, m_B^2, m_S^2)
    + Im(sa^*) 
    h_2(m_A^2, m_B^2, m_S^2)
    ]
    ,
    
    
     (0,2)-diagram
    =   
    i 6/(4π)^4 g_s^3 (X X^† T_B T_B )  O
       ×[
    i Re(sa^*)/2 
    h_1(m_B^2, m_A^2, m_S^2)
    + Im(sa^*) 
    h_2(m_B^2, m_A^2, m_S^2)
    ]
    ,
    
    
     (1,1)-diagram
    =   
    i 6/(4π)^4 g_s^3 (X T_A X^† T_B )  O
       ×[
    i Re(sa^*)/2 
    h_3(m_A^2, m_B^2, m_S^2)
    + Im(sa^*)
    h_4(m_A^2, m_B^2, m_S^2)
    + (A ↔ B)
    
    
    
    ]
    ,




where




    h_1(m_A^2, m_B^2, m_S^2)
    =   ( i/(4π)^2)^-1∫_ℓℓ^2/(ℓ^2 - m_A^2)^4
    (B̃_0 + B̃_1)_(ℓ^2, m_S^2, m_B^2)
    ,
    
    h_2(m_A^2, m_B^2, m_S^2)
    =   ( i/(4π)^2)^-1∫_ℓm_A m_B/(ℓ^2 - m_A^2)^4B̃_0_(ℓ^2, m_S^2, m_B^2) 
    ,
    
    h_3(m_A^2, m_B^2, m_S^2)
    =   ( i/(4π)^2)^-1∫_ℓℓ^2/(ℓ^2 - m_B^2)^3∂ (B_0 + B_1)_(ℓ^2, m_S^2, m_A^2)/∂ m_A^2
    ,
    
    h_4(m_A^2, m_B^2, m_S^2)
    =   ( i/(4π)^2)^-1∫_ℓm_A m_B/(ℓ^2 - m_B^2)^3∂ B_0_(ℓ^2, m_S^2, m_A^2)/∂ m_A^2
    ,
    
     B_0_(ℓ^2, m_S^2, m_A^2)
    =   ( i/(4π)^2)^-1∫_q 1/[q^2 - m_S^2][(q+ℓ)^2 - m_A^2]
     ,
    ℓ^μ B_1_(ℓ^2, m_S^2, m_A^2)
     =   ( i/(4π)^2)^-1∫_q q^μ/[q^2 - m_S^2][(q+ℓ)^2 - m_A^2]
    ,




where ∫_ℓ = ∫d^4 ℓ/(2π)^4,
and 
B̃_i = B_i +(counter terms).
For example, 
B̃_i_(ℓ^2, m_S^2, m_B^2) = B_i_(ℓ^2, m_S^2, m_B^2) - B_i_(m_A^2, m_S^2, m_B^2)
with the on-shell renormalization condition.
However, the contributions from counterterms vanish after we add other diagrams discussed in the next section.
In the calculation for Eq. (<ref>), we have used relations given in Eqs. (<ref>) and (<ref>).






 §.§ Diagrams with three external gluon fields

We evaluate diagrams with three external gluon fields to obtain GGG̃ terms.
Because of the same reason as in DGDG̃ terms, 
diagrams should have at least two gluon field from fermion line. 
In addition, diagrams with one gluon field from the scalar field vanish
thanks to a virtue of the Fock-Schwinger gauge.


See appendix <ref> for the detail.
The diagrams we have to evaluate are shown in Fig. <ref>.
We call each diagram from left to right as (3,0), (2,1), (1,2), and (0,3), respectively.





There are some terms that consist of X and SU(3)_c generators.
Direct calculations shows that they are proportional to the structure constant of SU(3)_c.
We can introduce the following notations.

    X_B̅AS (T^a)_AA' (T^b)_A'A” (T^c)_A”A”' X^†_A̅”̅'̅BS̅≡    
    (X T_A T_A T_A X^†) i f^abc
    ,
    
    X_B̅AS (T^a)_AA' (T^b)_A'A” X^†_A̅”̅B'S̅ (T^c)_B'B≡    
    (X T_A T_A X^† T_B) i f^abc
    ,
    
    X_B̅AS (T^a)_AA' X^†_A̅'̅B'S̅ (T^b)_B'B” (T^c)_B”B≡    
    (X T_A X^† T_B T_B) i f^abc
    ,
    
    X_B̅AS X^†_A̅B'S̅ (T^b)_B'B” (T^b)_B”B”' (T^c)_B”'B≡    
    (X X^† T_B T_B T_B) i f^abc
    .




We find relations

    ( X T_A T_A T_A X^†)  O
    =   1/2( X T_A T_A X^†)  O
    ,
    
     ( X T_A T_A X^† T_B )  O
    =   1/2( X T_A X^† T_B )  O
    ,
    
     ( X T_A X^† T_B T_B )  O
    =   1/2( X T_A X^† T_B )  O
    ,
    
     ( X X^† T_B T_B T_B )  O
    =   1/2( X X^† T_B T_B )  O
    .



It is easy to prove these relations. For example,



    ( X T_A T_A T_A X^†) (-6 O)
    =   
     -i ( X T^a T^b T^c X^†) G^a_μν G^b^ν_ ρ G^c_αβϵ^ρμαβ
    
    =   
    -i/2( X [T^a, T^b] T^c X^†) G^a_μν G^b^ν_ ρ G^c_αβϵ^ρμαβ
    
    =   1/2 f^abd( X T^d T^c X^†) G^a_μν G^b^ν_ ρ G^c_αβϵ^ρμαβ
    
    =   1/2( X T_A T_A X^†) (-6  O)
    .

Here we used a fact that 
G^a_μν G^b^ν_ ρϵ^ρμαβ is anti-symmetric under exchange of a and b.
Using the relations, we find 




    (3,0)-diagram
    =   
    -i 6/(4π)^4 g_s^3 ( X T_A T_A X^†)  O
         ×[
    iRe(sa^*)/2 
    h_1(m_A^2, m_B^2, m_S^2)
    
             
    + Im(sa^*)
    (
    h_2(m_A^2, m_B^2, m_S^2)
    +
    m_A m_B
    f_1(m_A^2, m_B^2, m_S^2)
    )
    ]
    ,
    
    
     (0,3)-diagram
    =   
    -i 6/(4π)^4 g_s^3 ( X X^† T_B T_B )  O
         ×[
    i Re(sa^*)/2
    h_1(m_B^2, m_A^2, m_S^2)
    
             
    + Im(sa^*)
    (
    h_2(m_B^2, m_A^2, m_S^2)
    +
    m_A m_B
    f_1(m_B^2, m_A^2, m_S^2)
    )
    ]
    ,
    
    
    (2,1)-diagram
    +
     (1,2)-diagram
    =   
    -i 6/(4π)^4 g_s^3 ( X T_A X^† T_B )  O
         ×[
     iRe(sa^*)/2
    h_3(m_A^2, m_B^2, m_S^2)
    
             
    + Im(sa^*)
    (
    h_4(m_A^2, m_B^2, m_S^2)
    +
    m_A m_B
    f_2(m_A^2, m_B^2, m_S^2)
    )
    
             
    +
    (A ↔ B)
    ]
    ,




where h_i are defined in Eqs. (<ref>)–(<ref>), and 




    f_1(m_A^2, m_B^2, m_S^2)
    ≡   ( i/(4π)^2)^-1∫_ℓℓ^2/(ℓ^2 - m_A^2)^4∂ B_0_(ℓ^2, m_S^2, m_B^2)/∂ℓ^2
    
    =   ∫_0^∞ dℓ_E^2
    ∫_0^1 dz
    -ℓ_E^4 z (1-z)/( m_S^2 z + m_B^2 (1-z) + ℓ_E^2 z (1-z) ) ( ℓ_E^2 + m_A^2 )^4
    ,
    
    f_2(m_A^2, m_B^2, m_S^2)
    ≡   ( i/(4π)^2)^-1∫_ℓℓ^2/(ℓ^2 - m_A^2)^4∂ B_0_(ℓ^2, m_S^2, m_B^2)/∂ m_B^2
    
    =   ∫_0^∞ dℓ_E^2
    ∫_0^1 dz
    ℓ_E^4 (1-z)/( m_S^2 z + m_B^2 (1-z) + ℓ_E^2 z (1-z) ) ( ℓ_E^2 + m_A^2 )^4
    .




Similar to (2,0) and (0,2)-diagrams, (3,0) and (0,3)-diagrams contain UV-divergences,
and thus h_1 and h_2 contain counterterms.
In the derivation of f_1, we have used a relation given in Eq. (<ref>).




§ RESULTS


We sum up all of the diagrams calculated in Sec. <ref>,
and find the following coefficient of the Weinberg operator.




    C_G
    =   
    -1/(4π)^4 6 
    Im(sa^*) m_A m_B 
    
       ×{(X T_A T_A X^†) 
    f_1(m_A^2, m_B^2, m_S^2)
    +
    (X X^† T_B T_B ) 
    f_1(m_B^2, m_A^2, m_S^2)
    
    
            
    +
    (X T_A X^† T_B )
    [
    f_2(m_A^2, m_B^2, m_S^2) + f_2(m_B^2, m_A^2, m_S^2)
    ]
    }
    ,




where f_1 and f_2 are defined in Eqs. (<ref>) and (<ref>).
From the definitions, we see that f_1 < 0 and f_2 >0.

If we use the same interaction term twice, 
we have to multiply 1/2 by Eq. (<ref>).
For instance, if ψ_A=ψ_B and S^* = S, then Eq. (<ref>)
is modified with extra 1/2 as follows.

    C_G
    =   
    -1/21/(4π)^4 6 
    Im(sa^*) m^2_A  
    
       ×{(X T_A T_A X^†) 
    f_1(m_A^2, m_A^2, m_S^2)
    +
    (X X^† T_A T_A ) 
    f_1(m_A^2, m_A^2, m_S^2)
    
            
    +2
    (X T_A X^† T_A )
    [
    f_2(m_A^2, m_A^2, m_S^2)
    ]
    }
    .





We make two comments. First one is that the cancellation of Re(sa^*) terms. 
This cancellation is expected because these terms break hermicity of the Lagrangian. 
The second comment is the cancellation of the contributions from the counterterms,
which are in h_1 and h_2. Thus the result is independent from the renormalization
conditions for the subdiagrams in (2,0), (0,2), (3,0), and (0,3)-diagrams.


For numerical evaluations, we need to calculate group factors such as X T_A T_A X^†.

Table <ref> shows some examples.
Note that 
if there is an SU(3)_c singlet field, then
the non-zero factors are N(r) that is defined through (T^a T^b) = N(r) δ^ab.





Our result is model independent in the sense that it can be applied to any models that generate the Weinberg operator from the diagrams shown in Fig. <ref>. 


For example, we can apply Eq. (<ref>) to models that predict CP-violating quark-squark-gluino coupling by taking A =q, B = g̃, and S = q̃. 
Another example is a CP-violating tth-coupling, where we can use Eq. (<ref>) by taking A = B = t and S = h. 
We have checked that our result is consistent with the results reported in Refs. <cit.>.






§ A VIEW POINT FROM EFFECTIVE THEORY


In this section, we discuss our result from a view point of effective theory.
For m_A ≪ m_B, m_S,
the Weinberg operator can be 
generated from the chromo-electric dipole moment (cEDM) of the particle A 
that is generated by the heavy particles (B and S).
The authors in Refs. <cit.>
have calculated the Weinberg operator in the framework of effective theory assuming the existence of the cEDM 
at the UV cutoff scale for a fermion in the SU(3)_c fundamental representation, and they have found




    C_g^eff. =   d^c/32π^2 m_A,

where d^c is defined through

    L⊃   Q̅( -i/2 g_s d^c γ_5 σ^μν G_μν) Q.




Here Q is in the SU(3)_c fundamental representation.
In the following, we show that our full two-loop result agree with C_g^eff. in m_A ≪ m_B, m_S regime.


In order to compare our result given in Eq. (<ref>)
with Eq. (<ref>),
we need to calculate the cEDM of A, which is generated by the diagrams shown in Fig. <ref>.[
We calculate the diagrams shown in Fig. <ref> with the 't Hooft Feynman gauge, not with the Fock-Schwinger gauge.]







We introduce the following notations,




    ( X^†_A̅BS (T^a)_S'S X_B̅A'S')
    ≡   ( X^† T_S X) (T^a)_AA'
    ,
    ( X^†_A̅BS (T^a)_BB' X_B̅'̅A'S)
    ≡   ( X^† T_B X) (T^a)_AA'
    .




We denote the contribution to the cEDM from the diagram in the left (right) panel of Fig. <ref> as
d^c_S (d^c_B).
We find




    d^c_S =   
    - 1/(4π)^2 2Im(s a^*) m_B
    ( X^† T_S X)  
    ∫_0^1 dz z(1-z)/m_S^2 z + m_B^2 (1-z) - m_A^2 z(1-z)
    ,
    
     d^c_B =   
    - 1/(4π)^2 2Im(s a^*) m_B
    ( X^† T_B X)  
    ∫_0^1 dz (1-z)^2/m_S^2 z + m_B^2 (1-z) - m_A^2 z(1-z)
    .




The result depends on T_S as we can see from Eq. (<ref>).

However, T_A, T_B, and T_S must be related to each other by the gauge invariance,
and thus T_S can be written by T_A, T_B, and X. 
We find 

    X^† T_S X =   
    X^† T_B X - 1/N(r_A) (X T_A T_A X^†)
    
    ,
    
     X^† T_B X =   1/N(r_A) (X T_A X^† T_B)
    
    .




We take N(r_A) = 1/2
because Eq. (<ref>) is derived with a fermion in the SU(3)_c fundamental representation
in Refs. <cit.>.


A derivation of Eq. (<ref>) is given in Appendix <ref>.
Eq. (<ref>) is derived by multiplying (T^b)_A'A and Eq. (<ref>).
Using these relations, we can eliminate T_S from Eq. (<ref>), 
and results can be expressed only by X T_A T_A X^† and X T_A X^† T_B.
We find




    C_g^eff. =    
    -
    1/(4π)^4 6 Im(s a^*) m_A m_B
    
    
       ×(
    - (X T_A T_A X^†)
    1/3 m_A^2∫_0^1 dz z(1-z)/m_S^2 z + m_B^2 (1-z) - m_A^2 z(1-z)
           
    +
     (X T_A X^† T_B)
    1/3 m_A^2∫_0^1 dz (1-z)/m_S^2 z + m_B^2 (1-z) - m_A^2 z(1-z))
    .




We compare our full two-loop result, C_g,
with the result calculated based on the effective theory, C_g^eff..
We consider three cases,
(A,B,S) ∼ (3, 1, 3̅), (3, 3, 1), (3, 3, 8).
We calculate the ratio of C_g^eff. and C_g
by using Eqs. (<ref>), (<ref>), and Table. <ref>. 
Numerical results are shown in Fig. <ref>.
Here we take m_A = 100 GeV, and m_S = m_B.
We find that our calculation agrees with the result in Refs. <cit.> very well for m_A ≪ m_B, m_S.
This result means the validity of the calculation based on the effective theory for the large mass hierarchy regime.
On the other hand, we also find discrepancies between C_g^eff. and C_g for small mass hierarchy regime.
In particular, the discrepancy is not negligible for (A,B,S) ∼ (3, 1, 3̅).
Therefore, it is better to use C_g instead of C_g^eff. for the purpose of reducing theoretical uncertainty.








Even for the large mass hierarchy regime, C_g has an advantage over C_g^eff..
Eq. (<ref>) is applicable 
only if
a lighter particle is a color non-singlet fermion.
On the other hand, Eq. (<ref>) can be applied for any color representation with arbitrary mass spectra.









§ NUMERICAL ANALYSIS





Let us next estimate the nucleon EDMs induced by the Weinberg operator and discuss the importance of measuring EDMs. 
The value of EDMs for proton and neutron induced by the Weinberg operator (d_p(w) and d_n(w)) can be estimated by applying the QCD sum rules <cit.>:

    d_N(w)
    =
    ± e Λ_nEDM w(1 ),
        (N=n,p),

where Λ_nEDM=10-30 MeV and w(1 ) denotes w=g^3_sC_G at the scale of 1GeV. 
It should be noted that the factor Λ_nEDM can also be estimated by naive dimensional analysis (NDA) <cit.> and the NDA predicts a factor of two larger value than the estimation using QCD sum rules <cit.>. Therefore, the following phenomenological analysis may be under-estimating the hadronic uncertainty.
We also note that we cannot determine the sign of the contribution from the Weinberg operator solely from the analysis based on QCD sum rules or the NDA.

In order to estimate the values of w(1 ), we need to take into account the effect of the renormalization group equation (RGE) evolution between the matching energy scale (μ_match) and the scale of 1 GeV. The RGE for w at the leading order is given as <cit.> 

    d/d lnμw(μ)
    =
    g^2_s(μ)/16π^2(N_C+2N_f)w(μ),

where N_C=3 and N_f are the number of colors and quark flavors, respectively. Solving the RGE and plugging the solution into Eq. (<ref>), we estimate the nucleon EDMs induced by the Weinberg operator.

The Weinberg operator induces quark EDMs and cEDMs through the operator mixing effects and they also affect the nucleon EDMs. However, the contributions are subdominant compared to that from the Weinberg operator because the operator mixing effects are suppressed by the one-loop factor. For the models discussed below, we have checked |d^(c)EDM_N/d_N(w)|≲𝒪(0.01) numerically where d^(c)EDM_N is the nucleon EDMs induced by quark (c)EDMs and d_N(w) is the nucleon EDMs induced by Weinberg operator. Therefore, we neglect the operator mixing effects in the following analysis.

We are now ready to estimate the value of d_N(w). Let us consider three cases, (A,B,S)∼(3,1,3), (3,3,1), and (3,3,8) as concrete examples.
In the numerical evaluation of the nucleon EDMs, we use 
m_Z=91.1876 , 
m_t=173.1 , 
m_b=4.18 , 
m_c=1.28 ,
and
α_s(m_Z)=0.1182
as input parameters <cit.> and set the matching energy scale, μ_match, as the minimum value among m_A, m_B, and m_S (μ_match=[m_A,m_B,m_S]).



Figure <ref> and figure <ref> show the absolute value of d_N(w) (|d_N(w)|) with taking m_S=m_B and (sa^∗)=0.25. In figure <ref>, blue, red, and green bands correspond to |d_N(w)/e| in the case where m_A=100 GeV, m_A=1 TeV, and m_A=10 TeV, respectively. In these colored bands, the upper lines, middle lines, and lower lines correspond to |d_N(w)/e| estimated with Λ_nEDM=30 MeV, Λ_nEDM=20 MeV, and Λ_nEDM=10 MeV, respectively. In figure <ref>, we show the contours of |d_N(w)/e| in (m_A, m_B) plane. In the left panels, we take Λ_nEDM=10 MeV. In the right panels, we take Λ_nEDM=30 MeV.











In m_A≫ m_B=m_S regime, we find that the mass dependence of |d_N(w)| in (A,B,S)∼(3,1,3) model is different from that in the other models. In (A,B,S)∼(3,1,3) model, |d_N(w)| behaves as |d_N(w)|∝ m_B/m^3_A in m_A≫ m_B=m_S regime. On the other hand, in the (A,B,S)∼(3,3,1) and (3,3,8) models, |d_N(w)| behaves as |d_N(w)| ∝1/(m_Am_B) in m_A≫ m_B=m_S regime. The difference of the mass dependence can be understood by taking into account the mass dependence of loop functions f_1 and f_2, which are summarized in Appendix <ref>.  


The region filled with magenta color in the figures shows the excluded region by the latest upper bound on neutron EDM, 

|d_n|≤ 2.9× 10^-26 e cm <cit.>. We find that the models with m_A≃ m_B≃ m_S≃𝒪(100) GeV and (sa^∗)=0.25 can be severely constrained by the current upper limit on the neutron EDM. 
For example, 
if we take (sa^∗)=0.25 and m_B=m_S=100 GeV in (3,3,1) model, 
m_A should be larger than
245-787 GeV
to avoid the latest upper bound on neutron EDM.
The uncertainty of the lower mass bound comes from the uncertainty of the scale of Λ_nEDM in Eq. (<ref>). 

It should be noted that the existence of 𝒪(100) GeV colored particles might be also severely constrained by the hadron collider experiments. However, in order to estimate the collider bounds on our models explicitly, we need to specify the production and decay properties of the colored particles. The properties of colored particles highly depend on the detail of the models which is not specified in this paper. Thus, 

we defer the estimation as future work.


Magenta dashed lines show the future prospects for the observation of nucleon EDMs. Here we take |d_n|=5×10^-28 <cit.> and |d_p|=2.5×10^-29 <cit.> as reference values of the future sensitivities. 
It is found that the future measurements for the EDMs may have sensitivities to the mass scale of 𝒪(1-10) TeV in each scenario. 
For example, 
if we take (sa^∗)=0.25 and m_B=m_S=1 TeV in (3,3,1) model, the future measurement of proton EDM proposed by Ref. <cit.> can cover m_A≲15.6-47.1 TeV mass region.
It is implied that future measurements of EDMs enable us to investigate physics at higher energy scale beyond the reach of collider experiments.



Before closing this section, we would like to give a comment on the uncertainty of our numerical calculation.
The main uncertainty of our numerical analysis comes from 
the QCD sum rules which we use in order to estimate the nucleon EDMs induced by the Weinberg operator. If the value of Λ_nEDM is determined more precisely, we can reduce the uncertainty of our numerical calculation. It means that the improvement of the determination Λ_nEDM gives us more detailed information of physics beyond the SM.









§ SUMMARY




We have investigated the contribution to EDMs from a CP-violating purely gluonic dimension-6 operator
called the Weinberg operator, GGG̃. 
We have calculated its Wilson coefficient at the leading order in a setup 
that contains CP-violating interactions constructed by one scalar and two fermion fields.
The SU(3)_c gauge invariance is imposed in the interaction terms.
The SU(3)_c representations of the scalar and fermion fields are kept arbitrary in our analysis,
and thus our result can be applied to a variety of models. 
Our main result is given in Eq. (<ref>).


We have compared our full two-loop result with the result calculated based on the effective theory approach <cit.>.
The effective theory is applicable if the lightest particle in the loop is a colored fermion.
We have numerically shown that our result agrees with the effective theory in the region of parameter space where
the lightest particle is much lighter than the other particles.
However, if the mass difference is small, the discrepancy between two calculations are sizable 
as can be seen from Fig. <ref>. In such a regime, we have to use the result of the full two-loop calculation.
It is worth noting that the full two-loop result given in Eq. (<ref>)
is as simple as the result derived from the effective theory Eq. (<ref>).




Finally, 
we have estimated the nucleon EDMs generated by the Weinberg operator, 
and have compared it with the current upper bound and future sensitivities on the neutron EDM.
We have found that all the newly introduced particles have to be heavier than O(100) GeV
due to the current upper bound on the EDM. 
We also have found that future experiments can search new particles whose masses are heavier than O(1) TeV.
For example, in (A,B,S) ∼ (3,3,1) model with m_B = m_S =1 TeV, it is possible to observe EDM 
for m_A < 16–47 TeV. 
Thus the future measurements of EDMs are powerful tools for probing physics at higher energy scale beyond the reach of collider experiments.
On the other hand, the uncertainty originates from the estimation of the EDM based on the naive dimensional analysis.
It is necessary to improve the evaluation of the contribution to EDMs from the Weinberg operator
for the determination of the new physics scale by the EDM measurements.[
For a current status of a lattice calculation of the neutron EDM from the Weinberg operator, see Ref. <cit.>.
]










§ ACKNOWLEDGMENTS

This work was supported by JSPS KAKENHI Grant Number 16K17715 [T.A.],
16H06490 [R.N.]
and by Grant-in-Aid for Scientific research from the Ministry of Education, Science, Sports, and Culture (MEXT), Japan, No. 16H06492 [J.H.].
The work of J.H. is also supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. 





Appendix



§ ONE GLUON FIELD FROM THE SCALAR FIELD IN FOCK-SCHWINGER GAUGE


We briefly review the Fock-Schwinger gauge and explain 
the reason why diagrams with gluon emission from the scalar field are absent in our calculation. 




We treat gluon fields as background fields obeying the Fock-Schwinger gauge,

    (x-x_0)^μ A_μ(x).

This gauge fixing condition allows us to expand the gauge field by its field strength.
In practice, we need to know the Feynman rules with this background field.
Vertices with fermion propagators is modified as




    i/p + k - m 
    (-ig_s γ^μ T^a ϵ^a_μ(k))
    i/p - m→   i/p + k - m 
    (-ig_s γ^μ T^a Ã^a_μ(k))
    i/p - m 
    ,




where




    Ã_μ(k) 
    =   (
    1/2 G_αμ(x_0) 1/i∂/∂ k_α
       
    +
    1/3 D_α_1 G_α_2 μ(x_0) 
    1/i∂/∂ k_α_11/i∂/∂ k_α_2
       
    +
    1/41/2! D_α_1 D_α_2 G_α_3 μ(x_0) 
    1/i∂/∂ k_α_11/i∂/∂ k_α_21/i∂/∂ k_α_3
       
    + ⋯) (2π)^4 δ(k)
    .




We also have to perform integration with respect to k.



Vertices with scalar propagators are similarly modified as follows.




    i/(p+k)^2 - m^2 
    ((2p + k)^μ g_s T^a ϵ^a_μ(k))
    i/p^2 - m^2 →    i/(p+k)^2 - m^2 
    ((2p + k)^μ g_s T^a Ã^a_μ(k))
    i/p^2 - m^2
    ,
    
    i/(p+k)^2 - m^2 
    (g_s^2 T^a T^b ϵ^a μ(k) ϵ^b_μ(k))
    i/p^2 - m^2 →    i/(p+k)^2 - m^2 
    (g_s^2 T^a T^b Ã^a μ(k) Ã^b_μ(k))
    i/p^2 - m^2 
    .




Using these rules, we can show the absence of diagrams with gluon emission from the scalar field.



We start by considering diagrams that generate DGDG̃. 
In this case, we need two gluon fields, Ã^a_μ(k) and Ã^b_ν(k').
If two of them are emitted from scalar fields,
we only have two gamma matrices from the fermion propagator, and thus
cannot obtain ϵ tensor. 
Moreover, we can show that 
diagrams with a gluon emitted from a fermion field and the other gluon from the scalar field
also vanish. Let us consider a diagrams with a gluon field emitted from fermion field A. 
The diagram is proportional to




    ∫_p,q,k(
    γ^5
    p+k + m_A/(p+k)^2 - m_A^2γ^μp + m_A/p^2 - m_A^2p+k+q + m_B/(p+k+q)^2 - m_B^2)
    × (⋯)
    
    
    =   
    -4i
    ∫_p,q,kϵ^αμβσ
    k_α p_β q_σ/
    ((p+k)^2 - m_A^2)
    (p^2 - m_A^2)
    ((p+k+q)^2 - m_B^2)
    × 
    (⋯)
    ,




where (⋯) is a term independent from p. 
The integral in Eq. (<ref>) is zero 
because of the symmetry in the p-integral and ϵ tensor.
Therefore, two gluon fields have to be emitted from fermion fields.
This is the reason for the absence of diagrams with gluon emission from the scalar field in Fig. <ref>.


We move to discuss diagrams that generate GGG̃. 
There are three gluon fields in this case.
We can show that at least two gluon fields have to be emitted from fermion fields 
because of the same reason in the previous paragraph.
The remains are diagrams with zero gluon fields from the scalar field and one gluon field from the scalar field.
The latter diagrams vanish thanks to a virtue of the Fock-Schwinger gauge.
The proof is the following.
From Eqs. (<ref>) and (<ref>), 
it is easy to find the scalar propagator with one gluon field strength,

    g_s/2 T^a G^a_ρμi/(p+k)^2 - m^2(
    (2p + k)^μ∂/∂ k_ρδ^(4)(k)
    )
    i/p^2 - m^2 
    .
 
Under the exchange of ρ and μ, G_ρμ is anti-symmetric
but δ^(4)(k) ∂/∂ k_ρ( [(p+k)^2 - m^2]^-1 (2p + k)^μ) 
is symmetric.
Thus we can simplify the expression as

    g_s/2 T^a G^a_ρμi/p^2 - m^2(
    2p^μ∂/∂ k_ρδ^(4)(k)
    )
    i/p^2 - m^2 
    .
 
This propagator vanishes if fermion propagators are independent from k.
The fermion propagators are independent from k if four-momentum is not conserved at one of the vertices.
The violation of the momentum conservation at a vertex happens in the calculation with the Fock-Schwinger gauge
because the gauge fixing condition breaks translational invariance. 
As a result, the fermion propagators are independent from k, and thus the diagrams with one gluon
emission from the scalar boson vanish.
This is the reason for the absence of diagrams with gluon emission from the scalar field in Fig. <ref>.








§ SOME FORMULAE OF INTEGRALS

In section <ref>, we have used the following relations.

    -4/3∫_ℓℓ^4/(ℓ^2 - m_A^2)^3∂^2 (B_11 + 2 B_1 + B_0)_(ℓ^2, m_S^2, m_B^2)/∂ (m_B^2)^2
    =
    ∫_ℓ2 ℓ^2 m_A^2/(ℓ^2 - m_A^2)^4∂ (B_0 + B_1)_(ℓ^2, m_S^2, m_B^2)/∂ m_B^2
    +
    (A ↔ B)
    
    ,
    
       ∫_ℓ2 ℓ^2/(ℓ^2 - m_A^2)^3∂^2 (B_0 + B_1)_(ℓ^2, m_S^2, m_B^2)/∂ (m_B^2)^2
    =
    - 
    (
    1/(ℓ^2 - m_A^2)^3 
    +
    3 m_A^2/(ℓ^2 - m_A^2)^4)
     ∂ B_0_(ℓ^2, m_S^2, m_B^2)/∂ m_B^2
     +
     (A ↔ B)
    
    ,
    
       ∫_ℓ( 2/(ℓ^2 - m_A^2 )^4 + 4 m_A^2/(ℓ^2 - m_A^2 )^5)
    B̃_0_(ℓ^2, m_S^2, m_B^2)
    =
    ∫_ℓℓ^2/(ℓ^2 - m_A^2 )^4∂/∂ℓ^2
    B_0(ℓ^2, m_S^2, m_B^2)
    .




To derive Eq. (<ref>), we use integration by parts with
the following relation.

    2/(ℓ^2 - m_A^2 )^4 + 4 m_A^2/(ℓ^2 - m_A^2 )^5
    =   
    - 1/2∂/∂ℓ^μℓ^μ/(ℓ^2 - m_A^2 )^4 
    .
 













§ APPROXIMATE FORMULAE


In this appendix, we derive approximate formulae for loop functions f_1 and f_2 which can be applied in m_A≪ m_B=m_S and m_A≫ m_B=m_S regime. The formulae help us to understand the results given in section <ref> qualitatively.

To begin with, we focus on the case where m_A≪ m_B=m_S. In this regime, we find

    f_1(m_A,m_B,m_B)   ≃ -1/18m^2_Am^2_B,
    
    f_1(m_B,m_A,m_B)   ≃ -1/18m^4_B,
    
    f_2(m_A,m_B,m_B)   +f_2(m_B,m_A,m_B)≃1/6m^2_Am^2_B.

where we ignore 𝒪(m^2_A/m^2_B) corrections. 
Applying above expressions, we obtain

    C_G
       ≃1/6(4π)^4Im(s a^*)1/m_Am_B,
        (A,B,S)∼(3,1,3),
    
    C_G
       ≃
    -1/3(4π)^4Im(s a^*)1/m_Am_B,
        (A,B,S)∼(3,3,1),
    
    C_G
       ≃11/36(4π)^4Im(s a^*)1/m_Am_B,
        (A,B,S)∼(3,3,8),
 
in m_A≪ m_B=m_S regime. On the other hand, in m_A≫ m_B=m_S regime, we find

    f_1(m_A,m_B,m_B)   ≃ -1/6m^4_A,
    
    f_1(m_B,m_A,m_B)   ≃ -1/6m^2_Am^2_B,
    
    f_2(m_A,m_B,m_B)   +f_2(m_B,m_A,m_B)≃1/3m^2_Am^2_B,

where we ignore 𝒪(m^2_B/m^2_A) corrections. Thus, we find that the Wilson coefficient of the Weinberg operator behaves as

    C_G
       ≃1/2(4π)^4Im(s a^*)m_B/m_A^3,
        (A,B,S)∼(3,1,3),
    
    C_G
       ≃
    -1/2(4π)^4Im(s a^*)1/m_Am_B,
        (A,B,S)∼(3,3,1),
    
    C_G
       ≃5/6(4π)^4Im(s a^*)1/m_Am_B,
        (A,B,S)∼(3,3,8),
 
in m_A≫ m_B=m_S regime, respectively. Note that the mass dependence of C_G in (3,1,3) is different from the other models.





§ A RELATION AMONG X^† T_S X, X^† T_B X, AND X T_A T_A X^†



The SU(3)_c gauge invariance impose the following relations.

    (U^†)_B'B (U)_AA' (U)_SS' X_B̅ A S =    X_B̅' A' S',




where U_AA', U_SS', and U^†_B'B are the unitary matrices 
under the SU(3)_c gauge transformation
for the fields A, S, and B, respectively.
Since X consists of ϵ tensors and Kronecker deltas, X is an invariant tensor,
and thus 

    (T^a)_AA' X_B̅'̅ A S' +  (T^a)_SS' X_B̅'̅ A' S - (T^a)_B'B X_B̅ A' S'
    =0.




Using this relation, we find




    (X^† T_S X) (T^a)_AA'
    =   
     (X^† T_B X) (T^a)_AA'
     - X^†_A̅ B S̅ (T^a) _A” A X_B̅A”S
    .




It is easy to show that 




    X^†_A̅ B S̅ (T^a) _A” A' X_B̅A”S
    =   1/N(r_A) (X T_A T_A X^†) (T^a)_AA',
 



where N(r_A) is defined through (T^a T^b) = N(r_A) δ^ab.
Finally, we find Eq. (<ref>).


























99

PhysLett.B109.215 
  M. B. Gavela, A. Le Yaouanc, L. Oliver, O. Pene, J. C. Raynal and T. N. Pham,
  
  Phys. Lett.  109B, 215 (1982).
  doi:10.1016/0370-2693(82)90756-0
  
  



PhysLett.109B.490 
  I. B. Khriplovich and A. R. Zhitnitsky,
  
  Phys. Lett.  109B, 490 (1982).
  doi:10.1016/0370-2693(82)91121-2
  
  



PhysLett.B197.556 
  B. H. J. McKellar, S. R. Choudhury, X. G. He and S. Pakvasa,
  
  Phys. Lett. B 197, 556 (1987).
  doi:10.1016/0370-2693(87)91055-0
  
  



NUPHA.B244.77 
  J. O. Eeg and I. Picek,
  
  Nucl. Phys. B 244, 77 (1984).
  doi:10.1016/0550-3213(84)90182-2
  
  



hep-ph/0504231 
  M. Pospelov and A. Ritz,
  
  Annals Phys.  318, 119 (2005)
  doi:10.1016/j.aop.2005.04.002
  [hep-ph/0504231].
  
  



Pospelov:1991zt 
  M. E. Pospelov and I. B. Khriplovich,
  
  Sov. J. Nucl. Phys.  53, 638 (1991)
  [Yad. Fiz.  53, 1030 (1991)].
  
  



hep-ex/0602020 
  C. A. Baker et al.,
  
  Phys. Rev. Lett.  97, 131801 (2006)
  doi:10.1103/PhysRevLett.97.131801
  [hep-ex/0602020].
  
  



1310.7534 
  J. Baron et al. [ACME Collaboration],
  
  Science 343, 269 (2014)
  doi:10.1126/science.1248213
  [arXiv:1310.7534 [physics.atom-ph]].
  
  



PRLTA.63.2333 
  S. Weinberg,
  
  Phys. Rev. Lett.  63, 2333 (1989).
  doi:10.1103/PhysRevLett.63.2333
  
  



PHRVA.D41.999 
  D. A. Dicus,
  
  Phys. Rev. D 41, 999 (1990).
  doi:10.1103/PhysRevD.41.999
  
  



PHLTA.B237.216 
  J. Dai, H. Dykstra, R. G. Leigh, S. Paban and D. Dicus,
  
  Phys. Lett. B 237, 216 (1990)
  Erratum: [Phys. Lett. B 242, 547 (1990)].
  doi:10.1016/0370-2693(90)91432-B, 10.1016/0370-2693(90)91810-X
  
  



PHRVA.D42.867 
  D. Chang, C. S. Li and T. C. Yuan,
  
  Phys. Rev. D 42, 867 (1990).
  doi:10.1103/PhysRevD.42.867
  
  



PHLTA.B249.467 
  I. Z. Rothstein,
  
  Phys. Lett. B 249, 467 (1990).
  doi:10.1016/0370-2693(90)91018-7
  
  



PRLTA.64.1709 
  E. Braaten, C. S. Li and T. C. Yuan,
  
  Phys. Rev. Lett.  64, 1709 (1990).
  doi:10.1103/PhysRevLett.64.1709
  
  



PHLTA.B241.584 
  G. Boyd, A. K. Gupta, S. P. Trivedi and M. B. Wise,
  
  Phys. Lett. B 241, 584 (1990).
  doi:10.1016/0370-2693(90)91874-B
  
  



PHLTA.B241.589 
  D. Chang, W. Y. Keung, C. S. Li and T. C. Yuan,
  
  Phys. Lett. B 241, 589 (1990).
  doi:10.1016/0370-2693(90)91875-C
  
  



PHLTA.B242.239 
  M. Dine and W. Fischler,
  
  Phys. Lett. B 242, 239 (1990).
  doi:10.1016/0370-2693(90)91464-M
  
  



Novikov:1983gd 
  V. A. Novikov, M. A. Shifman, A. I. Vainshtein and V. I. Zakharov,
  
  Fortsch. Phys.  32, 585 (1984).
  
  



1007.2601 
  J. Hisano, K. Ishiwata and N. Nagata,
  
  Phys. Rev. D 82, 115007 (2010)
  doi:10.1103/PhysRevD.82.115007
  [arXiv:1007.2601 [hep-ph]].
  
  



Mertig:1990an 
  R. Mertig, M. Bohm and A. Denner,
  
  Comput. Phys. Commun.  64, 345 (1991).
  doi:10.1016/0010-4655(91)90130-D
  
  



1601.01167 
  V. Shtabovenko, R. Mertig and F. Orellana,
  
  Comput. Phys. Commun.  207, 432 (2016)
  doi:10.1016/j.cpc.2016.06.008
  [arXiv:1601.01167 [hep-ph]].
  
  



Ellis:2008zy 
  J. R. Ellis, J. S. Lee and A. Pilaftsis,
  
  JHEP 0810, 049 (2008)
  doi:10.1088/1126-6708/2008/10/049
  [arXiv:0808.1819 [hep-ph]].
  
  



Jung:2013hka 
  M. Jung and A. Pich,
  
  JHEP 1404, 076 (2014)
  doi:10.1007/JHEP04(2014)076
  [arXiv:1308.6283 [hep-ph]].
  
  



Chang:1991ry 
  D. Chang, T. W. Kephart, W. Y. Keung and T. C. Yuan,
  
  Phys. Rev. Lett.  68, 439 (1992).
  doi:10.1103/PhysRevLett.68.439
  
  



Chang:1992vs 
  D. Chang, T. W. Kephart, W. Y. Keung and T. C. Yuan,
  
  Nucl. Phys. B 384, 147 (1992).
  doi:10.1016/0550-3213(92)90465-N
  
  



Demir:2002gg 
  D. A. Demir, M. Pospelov and A. Ritz,
  
  Phys. Rev. D 67, 015007 (2003)
  doi:10.1103/PhysRevD.67.015007
  [hep-ph/0208257].
  
  



Manohar:1983md 
  A. Manohar and H. Georgi,
  
  Nucl. Phys. B 234, 189 (1984).
  doi:10.1016/0550-3213(84)90231-1
  
  



Chang:1992bg 
  D. Chang, X. G. He, W. Y. Keung, B. H. J. McKellar and D. Wyler,
  
  Phys. Rev. D 46, 3876 (1992)
  doi:10.1103/PhysRevD.46.3876
  [hep-ph/9209284].
  
  



Chang:2000ck 
  D. Chang, X. G. He and B. H. J. McKellar,
  
  Phys. Rev. D 63, 096005 (2001)
  doi:10.1103/PhysRevD.63.096005
  [hep-ph/9909357].
  
  



Degrassi:2005zd 
  G. Degrassi, E. Franco, S. Marchetti and L. Silvestrini,
  
  JHEP 0511, 044 (2005)
  doi:10.1088/1126-6708/2005/11/044
  [hep-ph/0510137].
  
  



Patrignani:2016xqp 
  C. Patrignani et al. [Particle Data Group],
  
  Chin. Phys. C 40, no. 10, 100001 (2016).
  doi:10.1088/1674-1137/40/10/100001
  
  



Altarev:2009zz 
  I. Altarev et al.,
  
  Nucl. Instrum. Meth. A 611, 133 (2009).
  doi:10.1016/j.nima.2009.07.046
  
  



Lehrach:2012eg 
  A. Lehrach, B. Lorentz, W. Morse, N. Nikolaev and F. Rathmann,
  
  arXiv:1201.5773 [hep-ex].
  
  



1711.04730 
  J. Dragos, T. Luu, A. Shindler and J. de Vries,
  
  EPJ Web Conf.  175, 06018 (2018)
  doi:10.1051/epjconf/201817506018
  [arXiv:1711.04730 [hep-lat]].
  
  




 



