format: text

text:
# Contents

## 1 Parameters for a precision test

The running fine structure constant or QED coupling $`\alpha `$ at the scale of the $`Z^0`$ mass and the anomalous magnetic moment $`a_\mu `$ of the muon  are two parameters which have gained great deal of interest during the last years. The reason for this interest is that these parameters can be determined experimentally with high precision. Theoretical predictions, therefore, can serve as precision tests of the Standard Model (SM) of elementary particle physics because deviations to the measurement opens windows for the appearence of so-called “New Physics”.

### The QED coupling

For the QED coupling, an accurate knowledge of $`\alpha (M_Z)`$ is instrumental in narrowing down the mass window for the last particle of the Standard Model, the Higgs particle, which might have been already detected at the LEP collider at CERN at $`115\mathrm{GeV}`$ but which needs to be verified. It is worth to consider this parameter at the pole of the $`Z^0`$ because a huge amount of data have been taken there during a rather long period at the LEP I run. One could expect a rather precise experimental measurement at this point.

### The anomalous magnetic moment

The above mentioned window seems to open already for the anomalous magnetic moment of the muon. In February 2001, the Brookhaven National Laboratory (BNL) reported about a precision measurement of the anomalous magnetic moment of the positive muon ,

$$a_\mu ^{\mathrm{exp}}=(\mathrm{116\hspace{0.17em}592\hspace{0.17em}023}\pm 151)\times 10^{11}$$
(1)
which had to be contrasted with the theory prediction. At the first moment, a deviation of roughly $`2.6`$ standard deviations gave rise to suggestions for new-physics effects published during the last few months, including concepts of supersymmetry, leptoquarks, lepton number violating models, technicolor models, string theory concepts, extra dimensions and so on (cf. for instance Ref. ). However, the situation changed suddenly when an error in the sign was found in the calculation of the light-by-light contribution, obtained by Tom Kinoshita and collaborators  (see e.g. Refs. ). This sign was later on corrected by the authors themselves . The discrepancy was reduced again. New data sets taken in Novosibirsk again churn the water, but the physics community is split into two parts at least since these events.

### The hadronic contribution

From the theoretical point of view the uncertainty in the determination of $`\alpha (M_Z)`$ and $`a_\mu `$ is dominated by the uncertainty of the hadronic contribution $`\alpha ^{\mathrm{had}}`$ and $`a_\mu ^{\mathrm{had}}`$. This talk, therefore, will deal with these contributions only.

## 2 The integral representation

In this section I want to explain where the hadronic contribution comes from. Actually, the main ingredient to this contribution is the hadronic two-point correlator of the photon current, and the details and reasons of their origin are detailed in the following subsections.

### Hadronic contribution to the QED coupling

The corrections $`\mathrm{\Delta }\alpha (q^2)`$ to the QED coupling constant are given by the corrections to the photon propagator due to the calculation of a chain of inserted vacuum-polarization terms ,

$`{\displaystyle \frac{\alpha (q^2)}{q^2}}`$ $`=`$ $`\alpha \left({\displaystyle \frac{1}{q^2}}+{\displaystyle \frac{1}{q^2}}\mathrm{\Pi }_\gamma (q^2){\displaystyle \frac{1}{q^2}}+\mathrm{}\right)=`$ (2)
$`=`$ $`{\displaystyle \frac{\alpha }{q^2(1+\mathrm{\Pi }_\gamma (q^2)/q^2)}}={\displaystyle \frac{\alpha }{q^2(1\mathrm{\Delta }\alpha (q^2))}}`$
where $`\mathrm{\Pi }_\gamma (q^2)=e^2q^2\mathrm{\Pi }(q^2)=4\pi \alpha q^2\mathrm{\Pi }(q^2)`$ and $`\mathrm{\Pi }(q^2)`$ is the two-point correlator. This function is an analytical function except for a cut along the real axis for $`q^2<4m_\pi ^2`$, and it vanishes for $`|q^2|\mathrm{}`$. In order to understand what is meant by this, I take a simple example, namely

$$\mathrm{\Pi }_M(q^2)=\sqrt{4m_\pi ^2q^2}.$$
(3)
This function takes an imaginary value for $`q^2>4m_\pi ^2`$. The sign of the imaginary value depends on whether we approach the real axis from the upper or lower half plane, as it is shown in the figure on the right hand side. This should be indicated by $`q^2=se^{\pm i0}`$, $`s>4m_\pi ^2`$. So we obtain

$`\sqrt{4m_\pi ^2se^{i0}}`$ $`=`$ $`\sqrt{4m_\pi ^2+se^{i\pi }}=\sqrt{(s4m_\pi ^2)e^{i\pi }}=`$
$`=`$ $`e^{i\pi /2}\sqrt{s4m_\pi ^2}=i\sqrt{s4m_\pi ^2},`$
$`\sqrt{4m_\pi ^2se^{i0}}`$ $`=`$ $`e^{i\pi /2}\sqrt{s4m_\pi ^2}=i\sqrt{s4m_\pi ^2}.`$ (4)
We define the discontinuity by

$$\text{Disc }\mathrm{\Pi }_M(s)=\mathrm{\Pi }_M(se^{i0})\mathrm{\Pi }_M(se^{i0})=2i\sqrt{s4m_\pi ^2}.$$
(5)
The spectral density is finally given by

$$\rho _M(s)=\frac{1}{2\pi i}\text{Disc }\mathrm{\Pi }_M(s)=\frac{1}{\pi }\sqrt{s4m_\pi ^2}.$$
(6)
After having investigated the cut by introducing the discontinuity, we can use Cauchy’s theorem in order to express the correlator function by the corresponding spectral density. For this we take the circle about the specific point $`z=q^2`$ and expand this circle to a circle with infinite radius. If the circle reaches the cut, it will circumvent this cut along the real axis, resulting in two line integrals. In assuming that in the limit of an infinite radius the circle part of the integral vanishes,

$`\mathrm{\Pi }_M(q^2)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi i}}{\displaystyle \frac{\mathrm{\Pi }_M(z)dz}{zq^2}}=`$ (7)
$`=`$ $`{\displaystyle \frac{1}{2\pi i}}{\displaystyle _{\mathrm{}e^{i0}}^{4m_\pi ^2e^{i0}}}{\displaystyle \frac{\mathrm{\Pi }_M(z)dz}{zq^2}}+{\displaystyle \frac{1}{2\pi i}}{\displaystyle _{4m_\pi ^2e^{i0}}^{\mathrm{}e^{i0}}}{\displaystyle \frac{\mathrm{\Pi }_M(z)dz}{zq^2}}=`$
$`=`$ $`{\displaystyle \frac{1}{2\pi i}}{\displaystyle _{4m_\pi ^2}^{\mathrm{}}}{\displaystyle \frac{\mathrm{\Pi }_M(s^{i0})\mathrm{\Pi }_M(s^{i0})}{sq^2}}𝑑s={\displaystyle _{4m_\pi ^2}^{\mathrm{}}}{\displaystyle \frac{\rho _M(s)}{sq^2}}.`$
This identity, known as dispersion relation, is valid not only for the used example but generally if $`\mathrm{\Pi }_M(q^2)`$ falls off sufficiently fast for $`|q^2|0`$. However, this is not the case for the correlator function we have to deal with in our real application,

$$\mathrm{\Pi }(q^2)=_{4m_\pi ^2}^{\mathrm{}}\frac{\rho (s)ds}{s+q^2},\rho (s)=\frac{1}{2\pi i}\text{Disc }\mathrm{\Pi }(s),\text{Disc }\mathrm{\Pi }(s)=\mathrm{\Pi }(se^{i\pi })\mathrm{\Pi }(se^{i\pi })$$
(8)
(it is convenient to write the correlator function in Euclidean space). In this case $`\mathrm{\Pi }(q^2)`$ is singular. But we can redefine $`\mathrm{\Pi }(q^2)`$ by a subtracted quantity,

$$\mathrm{\Pi }(q^2)\mathrm{\Pi }(q^2)\mathrm{\Pi }(0)=_{4m_\pi ^2}^{\mathrm{}}\left(\frac{1}{s+q^2}\frac{1}{q^2}\right)\rho (s)𝑑s=_{4m_\pi ^2}^{\mathrm{}}\frac{q^2\rho (s)ds}{s(s+q^2)}.$$
(9)
This subtraction method is known as momentum subtraction at a specified momentum square, in this case at the point $`q^2=0`$. Taking only the hadronic contribution to the correction, the spectral density is related to the relative hadronic cross section in $`e^+e^{}`$ annihilations

$$R=\frac{\sigma (e^+e^{}\mathrm{hadrons})}{\sigma (e^+e^{}\mu ^+\mu ^{})}$$
(10)
by $`R(s)=12\pi ^2\rho (s)`$. For low energies this relative cross section is determined by the pion form factor $`F_\pi `$,

$$R(s)=\frac{v_\pi ^3}{4}|F_\pi (s)|^2,v_\pi =\sqrt{1\frac{4m_\pi ^2}{s}}.$$
(11)
So the expression we deal with for the hadronic contribution to the fine structure correction is given by

$$\mathrm{\Delta }\alpha (M_Z^2)=4\pi \alpha _{4m_\pi ^2}^{\mathrm{}}H(s)\rho (s)𝑑s,H(s)=\frac{M_Z^2}{s(M_Z^2s)}$$
(12)
The integrand consists of the spectral density and a factor which is called weight function.

### Hadronic contribution to the anomalous magnetic moment

The anomalous magnetic moment can be determined by looking at the first order QED correction to the muon-photon vertex with an self energy insertion for the loop photon. Using again the dispersion relation, we obtain

$$a_\mu =4\pi ^2\left(\frac{\alpha }{\pi }\right)^2_{4m_\pi ^2}^{\mathrm{}}H(s)\rho (s)𝑑s$$
(13)
where the weight function is $`H(s)=K(s)/s`$ with

$$K(s)=_0^1\frac{z^2(1z)dz}{z^2+(1z)s/m^2}.$$
(14)
Details about how to obtain this function are shown in Appendix A. The integral can be calculated and gives rise to

$`K(s)`$ $`=`$ $`{\displaystyle \frac{(1+x^2)(1+x)^2}{x^2}}\mathrm{ln}(1+x)+{\displaystyle \frac{x^2(1+x)}{1x}}\mathrm{ln}(x)+`$ (15)
$`+{\displaystyle \frac{x^2}{2}}(2x^2)+{\displaystyle \frac{(1+x^2)(1+x)^2}{x^2}}\left(x+{\displaystyle \frac{x^2}{2}}\right)`$
where

$$x=\frac{1v}{1+v},v=\sqrt{1\frac{4m_\mu ^2}{s}}.$$
(16)
In this case $`K(s)`$ has a cut along the real axis so that the total singularity of the weight function $`H(s)=K(s)/s`$ at the origin is stronger than in the case of the QED coupling. If we again only take the hadronic contribution, $`\rho (s)`$ is meant to be the spectral density corresponding to the hadronic correlator function.

### Local and global duality

The question may arise why we cannot exclude the experimental values from our considerations at all and take the relation in Eq. (8) in order to calculate the spectral density from a correlator function calculated perturbatively. The reason is that there is actually an obstacle in using this relation. As depending on methods of functional analysis, the inverse of the dispersion relation is only valid if there are no poles encircled by the path in the complex plane, a condition which is necessary to obtain this relation. These poles can have their origin from weight functions in combination with the spectral density. This means that if there is such a weight function included in the integration of the spectral density, the inverse relation shown above is only valid globally and not locally. We call this global resp. local duality. Two paths out of this situation are shown in the following.

## 3 The polynomial adjustment method

A singularity in the encirled part of the complex plane occurs for the two parameters we want to look at. Therefore we cannot fully keep out the experimental measurements from our considerations. There is, however, a way to include them in an optimal way. This is done by our method  which will be presented in the following section. In the meantime this method has also been used by other authors (see e.g. Ref. ). The method is based on the fact that we can use global duality when the weight function is non-singular. This is the case for a polynomial function. Therefore, we might mimic the weight function by a polynomial function obeying different conditions which I will explain later. By adding and subtracting this polynomial function $`P_N(s)`$ of given order $`N`$ to the weight function $`H(s)`$, without any restrictions we obtain

$$_{s_a}^{s_b}\rho (s)H(s)𝑑s=_{s_a}^{s_b}\rho (s)\left(H(s)P_N(s)\right)𝑑s+_{s_a}^{s_b}\rho (s)P_N(s)𝑑s$$
(17)
where $`[s_a,s_b]`$ is any interval out of the total integration range. But because the second term has now a polynomial weight, we can use global duality to write

$`{\displaystyle _{s_a}^{s_b}}\rho (s)P_N(s)𝑑s={\displaystyle \frac{1}{2\pi i}}{\displaystyle _{s_a}^{s_b}}\text{Disc }\mathrm{\Pi }(s)P_N(s)𝑑s=`$ (18)
$`=`$ $`{\displaystyle \frac{1}{2\pi i}}{\displaystyle _{|s|=s_a}}\mathrm{\Pi }(s)P_N(s)𝑑s+{\displaystyle \frac{1}{2\pi i}}{\displaystyle _{|s|=s_b}}\mathrm{\Pi }(s)P_N(s)𝑑s.`$
Thus this part can be represented by a difference of two circle integrals in the complex plane. On the other hand, the difference $`H(s)P_N(s)`$ suppresses the contribution of the first part. Our method consists thus of the following steps:

* replacing $`\rho (s)`$ in the first part of Eq. (17) by the value of the experimentally measured total cross section $`R(s)`$ (see e.g. Ref. ),
* replacing the circle integral contribution to flavours at their threshold by zero,
* in all other cases inserting the QCD perturbative and non-perturbative parts of $`\mathrm{\Pi }(s)`$ on the circle.

These replacements can be seen as a concept within QCD sum rules. To obtain the best efficience of our method, we have to restrict the polynomial function by the following contraints:

* The method of least squares should be used to mimic the weight
* However, the degree $`N`$ should not be higher than the order of the highest perturbative resp. non-perturbative contribution increased by one (this is a consequence of the Cauchy’s theorem which is involved in the analytical integration of the circle integrals).
* Especially for the low energy region, the polynomial function should vanish on the real axis to avoid instanton effects.
* In regions where resonances occur, the polynomial function should fit the weight function to suppress those contributions which constitute the highest uncertainty of the experimental data.

What are the ingredients for this mixture then? We have to look at what experimental data we can use and how the perturbative expansions for the correlator function are given.

### The experiment side

Data from experiments are crucial in the low energy range and in threshold regions where perturbative QCD cannot be applied. Here combined data sets from various electron positron annihilation experiments are used  which are accomplished by recent BES measurements .

In addition the use of precise $`\tau `$-decay data from Ref.  by isospin rotation promises to be a rewarding step in the low energy region. The vector $`\tau `$ spectral functions are related to the isovector $`e^+e^{}`$ cross sections for the corresponding hadronic states $`X`$ by

$`\sigma ^{I=1}(e^+e^{}X^0)={\displaystyle \frac{4\pi \alpha ^2}{s}}v_{J=1}(\tau ^{}X^{}\nu _\tau ).`$ (19)
$`v_{J=1}(\tau ^{}X^{}\nu _\tau )`$ is obtained by dividing the normalized invariant mass-squared distribution $`dN_X^{}/N_X^{}ds`$ for a given hadronic mass $`\sqrt{s}`$ by the appropriate kinematic factor,

$`v_{J=1}(\tau ^{}X^{}\nu _\tau )`$ $`=`$ $`{\displaystyle \frac{m_\tau ^2}{6|V_{ud}|^2S_{\mathrm{EW}}}}{\displaystyle \frac{B(\tau ^{}X^{}\nu _\tau )}{B(\tau ^{}e^{}\overline{\nu }_e\nu _\tau )}}\times `$ (20)
$`\times {\displaystyle \frac{dN_X^{}}{N_X^{}ds}}\left[\left(1{\displaystyle \frac{s}{m_\tau ^2}}\right)^2\left(1+{\displaystyle \frac{2s}{m_\tau ^2}}\right)\right]^1`$
where $`|V_{ud}|=0.9752\pm 0.0007`$ denotes the CKM weak mixing matrix element  and $`S_{EW}=1+\delta _{EW}=1.0194`$ accounts for electroweak second order corrections . The spectral functions are normalized by the ratio of the respective vector branching fraction $`B(\tau ^{}X^{}\nu _\tau )`$ to the branching fraction of the electron channel $`B(\tau ^{}e^{}\overline{\nu }_e\nu _\tau )=17.79\pm 0.04`$ .

### The theory side

The two-point correlator  is given by

$$i0|j_\alpha ^{\mathrm{em}}(x)j_\beta ^{\mathrm{em}}(0)|0e^{iqx}d^4x=(g_{\alpha \beta }q^2+q_\alpha q_\beta )\mathrm{\Pi }(q^2)$$
(21)
where we only included the isospin contribution $`I=1`$, in contrast to corresponding considerations for the $`\tau `$ decay. The scalar correlator function $`\mathrm{\Pi }(q^2)`$ consists of perturbative and non-perturbative contributions which we include to the extend we need them to keep the accuracy. For the perturbative contribution to the correlator we use a result given in Ref. . I only write down the first few terms,

$$\mathrm{\Pi }^\mathrm{P}(q^2)=\frac{3}{16\pi ^2}\underset{i=1}{\overset{n_f}{}}Q_i^2\left[\frac{20}{9}+\frac{4}{3}L+C_F\left(\frac{55}{12}4\zeta (3)+L\right)\frac{\alpha _s}{\pi }+O(\alpha _s^2,m_q^2/q^2)\right]$$
(22)
with $`L=\mathrm{ln}(\mu ^2/q^2)`$ while in Ref.  the expression is given up to the order $`O(\alpha _s^2,m_q^{12}/q^{12})`$. The number of active flavours is denoted by $`n_f`$. For the zeroth order term in the $`m_q^2/q^2`$ expansion we have added higher order terms in $`\alpha _s`$,

$$\frac{3}{16\pi ^2}\underset{i=1}{\overset{n_f}{}}Q_i^2\left[\left(c_3+k_2L+\frac{1}{2}(k_0\beta _1+2k_1\beta _0)L^2+\frac{1}{3}k_0\beta _0^2L^3\right)\left(\frac{\alpha _s}{\pi }\right)^3+O(\alpha _s^4)\right]$$
(23)
with $`k_0=1`$, $`k_1=1.63982`$ and $`k_2=6.37101`$. We have denoted the yet unknown constant term in the four-loop contribution by $`c_3`$. Remark, however, that the constant non-logarithmic terms will not contribute to our calculations. The non-perturbative contributions are given in ref. ,

$`\mathrm{\Pi }^{\mathrm{NP}}(q^2)={\displaystyle \frac{1}{18q^4}}\left(1+{\displaystyle \frac{7\alpha _s}{6\pi }}\right){\displaystyle \frac{\alpha _s}{\pi }}G^2+`$
$`+{\displaystyle \frac{8}{9q^4}}\left(1+{\displaystyle \frac{\alpha _s}{4\pi }}C_F+\mathrm{}\right)m_u\overline{u}u+{\displaystyle \frac{2}{9q^4}}\left(1+{\displaystyle \frac{\alpha _s}{4\pi }}C_F+\mathrm{}\right)m_d\overline{d}d+`$
$`+{\displaystyle \frac{2}{9q^4}}\left(1+{\displaystyle \frac{\alpha _s}{4\pi }}C_F+(5.8+0.92L){\displaystyle \frac{\alpha _s^2}{\pi ^2}}\right)m_s\overline{s}s+`$
$`+{\displaystyle \frac{\alpha _s^2}{9\pi ^2q^4}}(0.6+0.333L)m_u\overline{u}u+m_d\overline{d}d+`$
$`{\displaystyle \frac{C_Am_s^4}{36\pi ^2q^4}}\left(1+2L+(0.7+7.333L+4L^2){\displaystyle \frac{\alpha _s}{\pi }}\right){\displaystyle \frac{448\pi }{243q^6}}\alpha _s|\overline{q}q|^2+O(q^8)`$
where we have included the $`m_s^4/q^4`$-contribution arising from the unit operator. In this expression we used the $`SU(3)`$ colour factors $`C_F=4/3`$, $`C_A=3`$, and $`T_F=1/2`$. For the coupling constant $`\alpha _s`$ as well as for the running quark mass we use four-loop expression given in Refs.  even though in both cases the three-loop accuracy would already have been sufficient for the present application.

### The evaluation of the QED coupling

The calculation for the running fine structure constant or QED coupling is published in Ref. . An important step in preparing for the application of our method is to select points $`s_a`$ and $`s_b`$ for the limits of the integrals resp. the radii of the circles. Except for the threshold regions there seems to be a wide range for placing these points. Indeed we have shown that the method is fairly independent of this choice but is in some cases limited by computational constraints as for instance the fact that matrices fail to be inverted in some special cases. For this reason it is not advisable to make the intervals too narrow.

As an example for a threshold case I show the first contribution,

$$_{s_0}^{s_1}R(s)H(s)𝑑s=_{s_0}^{s_1}R^{\mathrm{exp}}(s)\left(H(s)P_N(s)\right)𝑑s+6\pi i_{|s|=s_1}\mathrm{\Pi }^{\mathrm{QCD}}(s)P_N(s)𝑑s.$$
(25)
Here we selected the range from the light flavour production threshold $`s_0=4m_\pi ^2`$ going to the next threshold marked by the mass of the $`\psi `$, $`s_1=m_\psi ^2(3.1\mathrm{GeV})^2`$. In this case of course the inner circle integral did not contribute. The polynomial weight functions are compared to the weight function itself in Fig. 2, the results for different polynomial degrees are shown in Fig. 2. Further subdivisions can be found in Tab. 1. The part of the integral starting from $`s_4=(40\mathrm{GeV})^2`$ up to infinity is done using local duality, i.e. by inserting the function $`R(s)`$ obtained for perturbative QCD,

$$\rho ^{\mathrm{had}}(s)=\frac{\alpha N_c}{12\pi }\underset{f}{}Q_f^2\sqrt{1\frac{4m_f^2}{s}}\left(1+\frac{2m_f^2}{s}\right)$$
(26)
($`N_c`$ is the number of colours) into the second part of Eq. (17).

To obtain the results shown in Tab. 1, we used the condensate values

$$\frac{\alpha _s}{\pi }GG=0.04\pm 0.04\mathrm{GeV}^4,\alpha _s\overline{q}q^2=(4\pm 4)\times 10^4\mathrm{GeV}^6.$$
(27)
For the errors coming from the uncertainty of the QCD scale we take

$$\mathrm{\Lambda }_{\overline{\mathrm{MS}}}=380\pm 60\mathrm{MeV}$$
(28)
The errors resulting from the uncertainty in the QCD scale in different energy intervals are clearly correlated and will have to be added linearly in the end. We also include the error of the strange quark mass in the light quark region which is taken as

$$\overline{m}_s(1\mathrm{GeV})=200\pm 60\mathrm{MeV}$$
(29)
For the charm and bottom quark masses we use the values

$$\overline{m}_c(m_c)=1.4\pm 0.2\mathrm{GeV},\overline{m}_b(m_b)=4.8\pm 0.3\mathrm{GeV}.$$
(30)
Summing up the contributions from the five flavours $`u`$, $`d`$, $`s`$, $`c`$, and $`b`$, our result for the hadronic contribution to the dispersion integral including the systematic error due to the dependence on $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$ (column 5 in Table 1) reads

$$\mathrm{\Delta }\alpha _{\mathrm{had}}^{(5)}(M_Z)=(277.6\pm 4.1)\times 10^4.$$
(31)
In order to obtain the total result for $`\alpha (M_Z)`$, we have to add the lepton and top contributions. Since we have nothing new to add to the calculation of these contributions we simply take the values from Ref. ,

$$\mathrm{\Delta }\alpha _{\mathrm{had}}^t(M_Z)=(0.70\pm 0.05)\times 10^4,\mathrm{\Delta }\alpha _{\mathrm{lep}}(M_Z)314.97\times 10^4.$$
(32)
Writing $`\mathrm{\Delta }\alpha (M_Z)=\mathrm{\Delta }\alpha _{\mathrm{lep}}(M_Z)+\mathrm{\Delta }\alpha _{\mathrm{had}}(M_Z)`$ our final result is

$$\alpha (M_Z)^1=\alpha (0)^1(1\mathrm{\Delta }\alpha (M_Z))=128.925\pm 0.056$$
(33)
where we used $`\alpha (0)^1=137.036`$.

### The evaluation of the magnetic moment

The procedure we choose for the anomalous magnetic moment of the muon differs from the one used for the QED coupling because of two reasons. As mentioned earlier, the singularity of the weight function at the origin on the one hand is much stronger and it therefore it would need more effort to approximate the weight by a polynomial in the lower energy region. On the other hand we can and should therefore make use of the fact that the $`\tau `$ data between $`s=(0.28\mathrm{GeV})^2`$ and $`s_1=(1.4\mathrm{GeV})^2`$ were measured with a high precision, so they should influence the evaluation of the dispersion integral at full extent. The less precise $`e^+e^{}`$ data from the above energy region, however, are worth being replaced by QCD expressions as much as possible.

Our results are presented in Table 2. The dealing with error estimates is at some point different from the precious case of the QED coupling because of the different data sets we used. We only keep error estimates steming from the uncertainties for the quark masses and a systematic error over all ranges from the uncertainty of the parameter $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$. In comparison to this, all other error estimates for the theory contribution of the “mixed” regions (like errors of vacuum expectation values) can be neglected. The value obtained is comparable with the predictions in Ref.  and can therefore be seen as one of the most accurate on this field. It will be published hopefully soon .

## 4 Estimates for the next-to-leading order

At this point we come back to the point whether the experimental measurements for the hadronic $`e^+e^{}`$ channel really have to be taken into account. The question is to be taken serious because in this case it is at least strange to denote the estimate “theory estimate”. The method to use experimental data for a theory estimate is different from the (in this respect more preferable) method to use the experimental data in order to extract a few parameters like the locations of resonances and their widths. In order to see whether such an approach is possible, we use simple models to extract a single parameter from the given value of the anomalous magnetic moment of the muon. On the first sight this looks like a “zero-sum game”, but using the model we then can estimate the next-to-leading order contribution to the anomalous magnetic moment which is found to be very close to data-based results in the literature.

### Integration in the Euclidean domain

As mentioned earlier, the non-analytic structure of the complex plane does not allow for the calculation of the spectral density (given in Minkowskian domain) by using perturbation theory results (obtained in Euclidean domain). Therefore, at this point we return to the Euclidean domain. The transition from

$$a_\mu =4\pi ^2\left(\frac{\alpha }{\pi }\right)^2_{4m_\pi ^2}^{\mathrm{}}\frac{K(s)}{s}\rho (s)𝑑s,K(s)=_0^1\frac{z^2(1z)dz}{z^2+(1z)s/m^2}$$
(34)
to the integral representation

$$a_\mu =4\pi ^2\left(\frac{\alpha }{\pi }\right)^2_0^{\mathrm{}}W(t)\left(\mathrm{\Pi }^{\mathrm{had}}(t)\right)𝑑t$$
(35)
where

$$W(t)=\frac{4m^4}{\sqrt{t^2+4m^2t}\left(t+2m^2+\sqrt{t^2+4m^2t}\right)^2}$$
(36)
is calculated in Appendix B. Finally, an integration-by-parts leads to

$$a_\mu =4\pi ^2\left(\frac{\alpha }{\pi }\right)^2_0^{\mathrm{}}\left(\frac{d\mathrm{\Pi }^{\mathrm{had}}(t)}{dt}\right)F(t)𝑑t$$
(37)
where

$$F(t)=_t^{\mathrm{}}W(t^{})𝑑t^{}=\frac{1}{2}\left(\frac{t+2m^2\sqrt{t^2+4m^2t}}{t+2m^2+\sqrt{t^2+4m^2t}}\right)=\frac{2m^4}{\left(t+2m^2+\sqrt{t^2+4m^2t}\right)^2}.$$
(38)
For high values of the Euclidean squared energy $`t`$ the expression $`d\mathrm{\Pi }^{\mathrm{had}}(t)/dt`$ can be taken from perturbation theory while for the whole range different models can be used which are introduced in the following. However, we can take a short look on the weight function $`F(t)`$. This function shown in Fig. 3 for $`m=m_\mu =(769.9\pm 0.8)\mathrm{MeV}`$  supports to very high extend low energy contributions. This will become important in the following.

### The step model

The simplest model we can find for the $`e^+e^{}`$ hadronic spectrum is the model of steps which appear when a new quark flavour is opened, i.e. when quark-antiquark pair production becomes kinematically possible at this specific energy. The model would simply read

$$\rho _1(s)=N_cQ^2\theta (s4m_Q^2)$$
(39)
where $`N_c=3`$ is the number of colours and $`Q`$ is the charge (in units of the elementary charge) of the quark flavour with mass $`m_Q`$. Actually we need not separate between the three lightest flavours but instead can take them in combination,

$$\rho _1^{uds}=N_c(Q_u^2+Q_d^2+Q_s^2)\theta (s4m_{\mathrm{eff}}^2)=N_cQ_{\mathrm{eff}}^2\theta (s4m_{\mathrm{eff}}^2),Q_{\mathrm{eff}}^2=\frac{2}{3}$$
(40)
where $`m_{\mathrm{eff}}`$ is the mass parameter which has to be adjusted in this model. For the step functions of the heavy quarks we use their masses from literature as cited before. Starting with this model spectral function, we obtain for the hadronic correlator function

$$\mathrm{\Pi }_1^{uds}(t)=t_{4m_Q^2}^{\mathrm{}}\frac{\rho _1(s)}{s(s+t)}𝑑s=N_cQ_{\mathrm{eff}}^2t_{4m_{\mathrm{eff}}^2}^{\mathrm{}}\frac{ds}{s(s+t)}=N_cQ_{\mathrm{eff}}^2\mathrm{ln}\left(\frac{4m_{\mathrm{eff}}^2}{t+4m_{\mathrm{eff}}^2}\right)$$
(41)
and for it’s derivative

$$\frac{d\mathrm{\Pi }_1^{uds}(t)}{dt}=\frac{N_cQ_{\mathrm{eff}}^2}{t+4m_Q^2}.$$
(42)
This actually makes sense, since perturbative QCD for high predicts

$$\frac{d\mathrm{\Pi }^{\mathrm{had}}(t)}{dt}=\frac{N_cQ_{\mathrm{eff}}^2}{t}\left(1+\frac{\alpha _s(t)}{\pi }\right)$$
(43)
While for small values of $`t`$ this expression remains finite.

### The threshold model

One can of course use more sophisticated models. One of these is given by the fermionic spectral density

$$\rho _2(s)=N_cQ_{\mathrm{eff}}^2\sqrt{1\frac{4m_q^2}{s}}\left(1+\frac{2m_q^2}{s}\right).$$
(44)
The corresponding hadronic correlator function is given by

$$\mathrm{\Pi }_2(t)=N_cQ_{\mathrm{eff}}^2\left\{\left(\frac{1}{z}3\right)\phi (z)\frac{1}{3}\right\}$$
(45)
where

$$\phi (z)=\frac{1}{\sqrt{z}}artanh(\sqrt{z})1,z=\frac{t}{t+4m_q^2},$$
(46)
and

$$\frac{d\mathrm{\Pi }_2(t)}{dt}=N_cQ_{\mathrm{eff}}^2\left\{\frac{t6m_q^2}{t^2}+\frac{24m_q^4}{t^3}\sqrt{\frac{t}{t+4m_q^2}}artanh\left(\sqrt{\frac{t}{t+4m_q^2}}\right)\right\}.$$
(47)
It figures out that the derivatives of $`\mathrm{\Pi }_1^{uds}(t)`$ and $`\mathrm{\Pi }_2(t)`$ coincide at the origin if we chose $`m_{\mathrm{eff}}/m_q=\sqrt{5}/2`$.

### The resonance model

Finally, one can also include the lowest lying resonance, namely the $`\rho `$ meson resonance. A one-scale no-parameter model that satisfies the duality constraints from the operator product expansion is given by

$$\rho _4(s)=N_cQ_{\mathrm{eff}}^2\left(2m_\rho ^2\delta (sm_\rho ^2)+\theta (s2m_\rho ^2)\right).$$
(48)
The derivative of the corresponding correlator function reads

$$\frac{d\mathrm{\Pi }_4(t)}{dt}=N_cQ_{\mathrm{eff}}^2\left(\frac{2m_\rho ^2}{(t+m_\rho ^2)^2}+\frac{1}{t+2m_\rho ^2}\right).$$
(49)
Alternatively, the narrow width $`\rho `$ resonance can be replaced by a Breit–Wigner resonance,

$$\rho _3(s)=N_cQ_{\mathrm{eff}}^2\left(\theta (2m_\rho ^2s)\rho _{\mathrm{BW}}(s)\theta (s4m_\pi ^2)+\theta (s2m_\rho ^2)\right)$$
(50)
where

$$\rho _{\mathrm{BW}}(s)=\frac{2m_\rho ^2}{\pi }\left(\frac{\mathrm{\Gamma }_\rho m_\rho }{(sm_\rho ^2+\mathrm{\Gamma }_\rho ^2/4)^2+\mathrm{\Gamma }_\rho ^2m_\rho ^2}\right),$$
(51)
the derivative of the correlator function is given by

$$\frac{d\mathrm{\Pi }_3(t)}{dt}=N_cQ_{\mathrm{eff}}^2\left(_{4m_\pi ^2}^{2m_\rho ^2}\frac{\rho _{\mathrm{BW}}(s)ds}{(s+t)^2}+\frac{1}{t+2m_\rho ^2}\right).$$
(52)
### The big surprize

The spectral densities we have used are quite different. A glance on Fig. 4 convinces us that at least the spectral density of the third model is quite different from the two others.

However, calculating the derivative of the correlator function for the three models, we obtain the result shown on the top of Fig. 5. The two curves of model 1 and model 2 coincide at the origin because we chose the relation $`m_{\mathrm{eff}}/m_q=\sqrt{5}/2`$. We can even adjust $`m_{\mathrm{eff}}/m_q`$ (and, therefore, $`m_q`$) in order to obtain a conicidence with the third curve at the origin. The surprizing result is that despite the fact that the spectral densities in Fig. 4 look quite different, all curves for the correlator derivatives have nearly the same shape. But because the weight function $`F(t)`$ supports contributions at low values of $`t`$ this means that the method is independent of the choice for the model!

The procedure, therefore, is as follows: we take the first model as the easiest one to obtain our results. The contributions to the anomalous magnetic moment of the muon for the heavy quark flavours can be calculated to be $`a_\mu ^{\mathrm{had}}(\mathrm{LO};\mathrm{heavy})=(71\pm 18)\times 10^{11}`$. The leading order result from the last chapter is reduced by this value to result in $`a_\mu ^{\mathrm{had}}(\mathrm{LO};\mathrm{light})=(6870\pm 72)\times 10^{11}`$. This value is used to adjust $`m_{\mathrm{eff}}`$. We obtain

$$m_{\mathrm{eff}}=(201.7\pm 1.2)\mathrm{MeV}.$$
(53)
### Next-to-leading order estimate

As I mentioned earlier, up to this point the method seems to be a “zero-sum game”. However, there is really makes sense if we want to calculate the next-to-leading order correction. The correction shown in Fig. 6(b) is the convolution of the spectral density with the two-loop kernel $`K^{(2)}(s)`$,

$$a_\mu ^{\mathrm{had}}(\mathrm{NLO}(\mathrm{b}))=4\pi \left(\frac{\alpha }{\pi }\right)^3_{4\pi ^2}^{\mathrm{}}\frac{K^{(2)}(s)}{s}\rho (s)𝑑s.$$
(54)
Using model 1, we obtain  $`a_\mu ^{\mathrm{had}}(\mathrm{NLO}(\mathrm{b}))=(211\pm 5)\times 10^{11}`$. For the two bubbles in Fig. 6(c), the so-called double-bubble diagram, model 1 can again employed to give  $`a_\mu ^{\mathrm{had}}(\mathrm{NLO}(\mathrm{c}))=(106\pm 2)\times 10^{11}`$. Both results, as well as their sum are very close to the data-based results of the literature . Finally, we add also the two next-to-leading order contributions, employing the four-point correlator, namely the light-by-light contribution (Fig. 6(d)) and the two photon Green function (Fig. 6(e)) (see the discussion in Ref. ) to obtain our estimate

$$a_\mu ^{\mathrm{had}}(\mathrm{NLO})=(85\pm 20)\times 10^{11}.$$
(55)
This estimate is still in good agreement with the difference

$`a_\mu ^{\mathrm{exp}}(\mathrm{NLO})`$ $`=`$ $`\left(7165\pm 151|_{\mathrm{exp}}\pm 2.9|_{\mathrm{QED}}\pm 4|_{\mathrm{EW}}(6941\pm 70)\right)\times 10^{11}=`$ (56)
$`=`$ $`(224\pm 189)\times 10^{11}`$
between the experimental value and the theory estimate for the LO correction. The improvement of the accuracy of the experiment will show whether this can be kept or not.

## 5 Conclusion

I have presented results for the QED coupling $`\alpha `$ at the $`Z^0`$ pole and the anomalous magnetic moment $`a_\mu `$ of the muon, showing the main features of the polynomial adjustment method. In the second part I have used the LO theory estimate to adjust a mass parameter and calculate the NLO contributions.

### Acknowledgements

I thank J.G. Körner, K. Schilcher, A.A. Pivovarov, N.F. Nasrallah, and J. Maul for a fruitful collaboration on this field. I acknowledge a grant given by the Graduiertenkolleg “Eichtheorien” in Mainz, Germany.

## Appendix A How does the function $`K(s)`$ appear?

If we look at the figure at the left hand side which is the diagram we have to consider, the correlator (in Minkowskian domain) shown as shaded bubble is given by

$$\mathrm{\Pi }_{\rho \sigma }(k^2)=(k_\rho k_\sigma k^2g_{\rho \sigma })\mathrm{\Pi }(k^2).$$
(57)
For the scalar correlator function we take

$$\mathrm{\Pi }(k^2)=\frac{k^2\rho (s)ds}{s(sk^2)}.$$
(58)
The contribution for the diagram, therefore, is given by

$`ie\mathrm{\Lambda }^\mu ={\displaystyle }{\displaystyle \frac{d^Dk}{(2\pi )^D}}\overline{u}(p^{})(ie\gamma ^\alpha ){\displaystyle \frac{i}{p/^{}+/km}}(ie\gamma ^\mu ){\displaystyle \frac{i}{p/+/km}}(ie\gamma ^\beta )u(p)\times `$ (59)
$`\times {\displaystyle \frac{i}{k^2}}(g_{\alpha \rho }(1\chi ){\displaystyle \frac{k_\alpha k_\rho }{k^2}})(ik^2)(g^{\rho \sigma }{\displaystyle \frac{k^\rho k^\sigma }{k^2}})\mathrm{\Pi }(k^2){\displaystyle \frac{i}{k^2}}(g_{\sigma \beta }(1\chi ){\displaystyle \frac{k_\sigma k_\beta }{k^2}})=`$
$`=`$ $`e^3{\displaystyle \frac{d^Dk}{(2\pi )^D}\overline{u}(p^{})\frac{\gamma ^\alpha (p/^{}+/k+m)\gamma ^\mu (p/+/k+m)\gamma ^\beta }{((p^{}+k)^2m^2)((p+k)^2m^2)k^2}\left(g_{\alpha \beta }\frac{k_\alpha k_\beta }{k^2}\right)\mathrm{\Pi }(k^2)}`$
while the anomalous magnetic moment is the part of the diagram proportional to $`(p+p^{})^\mu `$ for $`q^2=(pp^{})^2=0`$. The contraction with the second part of the effective transverse photon propagator simplifies the expression remarkably. Because of

$`\overline{u}(p^{})/k(p/^{}+/k+m)`$ $`=`$ $`\overline{u}(p^{})(2p^{}k+k^2+(p/^{}+m)/k)=u(p^{})(2p^{}k+k^2),`$
$`(p/+/k+m)/k`$ $`=`$ $`(2pk+k^2+/k(p/+m))u(p)=(2pk+k^2)u(p)`$ (60)
(note that $`\overline{u}(p^{})(p/^{}m)=(p/m)u(p)=0`$), we obtain the two denominator factors

$`(p^{}+k)^2m^2`$ $`=`$ $`p^2+2p^{}k+k^2m^2=2p^{}k+k^2,`$
$`(p+k)^2m^2`$ $`=`$ $`p^2+2pk+k^2m^2=2pk+k^2`$ (61)
($`p^2=p^2=m^2`$) which cancel out. However, the contribution figures out to be proportional to $`\gamma ^\mu `$ and therefore does not contribute to the anomalous magnetic moment. Therefore, we can contract with $`g_{\alpha \beta }`$ and simplify the Dirac structure.

### Simplification of the Dirac structure

The simplification is done systematically using $`\overline{u}(p^{})(p/^{}m)=(p/m)u(p)=0`$, $`\gamma ^\alpha \gamma _\alpha =D`$ and the usual relations for the Dirac matrices,

$`\overline{u}(p^{})\gamma ^\alpha (p/^{}+/k+m)\gamma ^\mu (p/+/k+m)\gamma _\alpha u(p)=`$ (62)
$`=`$ $`\overline{u}(p^{})(2(p^{}+k)^\alpha +(p/^{}/k+m)\gamma ^\alpha )\gamma ^\mu (2(p+k)_\alpha +\gamma _\alpha (p//k+m))u(p)=`$
$`=`$ $`\overline{u}(p^{})(2(p^{}+k)^\alpha /k\gamma ^\alpha )\gamma ^\mu (2(p+k)_\alpha \gamma _\alpha /k)u(p)=`$
$`=`$ $`4(pp^{}+(p+p^{})k+k^2)\overline{u}(p^{})\gamma ^\mu u(p)+\overline{u}(p^{})/k\gamma ^\alpha \gamma ^\mu \gamma _\alpha /ku(p)+`$
$`2\overline{u}(p^{})\gamma ^\mu (p/^{}+/k)/ku(p)2\overline{u}(p^{})/k(p/+/k)\gamma ^\mu u(p)=`$
$`=`$ $`(4m^2+4(p+p^{})k+4k^2)\overline{u}(p^{})\gamma ^\mu u(p)+(2D)\overline{u}(p^{})/k\gamma ^\mu /ku(p)+`$
$`4p^\mu \overline{u}(p^{})/ku(p)+2m\overline{u}(p^{})\gamma ^\mu /ku(p)2k^2\overline{u}(p^{})\gamma ^\mu u(p)+`$
$`4p^\mu \overline{u}(p^{})/ku(p)+2m\overline{u}(p^{})/k\gamma ^\mu u(p)2k^2\overline{u}(p^{})\gamma ^\mu u(p)=\mathrm{}`$
$`=`$ $`\overline{u}(p^{})((4m^2+4(p+p^{})k+(D2)k^2)\gamma ^\mu +`$
$`(4(p+p^{})^\mu +2(D2)k^\mu )/k+4mk^\mu )u(p)`$
where $`0=q^2=(pp^{})^2=2m^22pp^{}pp^{}=m^2`$ was used. Again, the first part is proportional to $`\gamma ^\mu `$ and does not contribute to the anomalous magnetic moment.

### Scalar, vector and tensor integrals

Inserting the dispersion relation for $`\mathrm{\Pi }(k^2)`$, the denominator factor $`k^2`$ is replaced by $`(sk^2)`$. The basic integral to deal with will be calculated by using Feynman parametrization,

$`I`$ $`=`$ $`{\displaystyle \frac{d^Dk}{(2\pi )^D}\frac{1}{(k^2+2pk)(k^2+2p^{}k)(k^2s)}}=`$ (63)
$`=`$ $`{\displaystyle \frac{d^Dk}{(2\pi )^D}\mathrm{\Gamma }(3)_0^1_0^{1z_1}\frac{dz_1dz_2}{(z_1(k^2+2pk)+z_2(k^2+2p^{}k)+(1z_1z_2)(k^2s))^3}}=`$
$`=`$ $`\mathrm{\Gamma }(3){\displaystyle _0^1}𝑑z_1{\displaystyle _0^{1z_1}}𝑑z_2{\displaystyle \frac{d^Dk}{(2\pi )^D}\frac{1}{(k^2+2(z_1p+z_2p^{})k(1z_1z_2)s)^3}}=`$
$`=:`$ $`\mathrm{\Gamma }(3){\displaystyle _0^1}𝑑z_1{\displaystyle _0^{1z_1}}𝑑z_2{\displaystyle \frac{d^Dk}{(2\pi )^D}\frac{1}{(k^2+2PkM^2)^3}}=`$
$`=`$ $`\mathrm{\Gamma }(3){\displaystyle _0^1}𝑑z_1{\displaystyle _0^{1z_1}}𝑑z_2{\displaystyle \frac{d^Dk}{(2\pi )^D}\frac{1}{((k+P)^2P^2M^2)^3}}=`$
$`=`$ $`\mathrm{\Gamma }(3){\displaystyle _0^1}𝑑z_1{\displaystyle _0^{1z_1}}𝑑z_2{\displaystyle \frac{d^Dk^{}}{(2\pi )^D}\frac{1}{(k^2P^2M^2)^3}}=`$
$`=`$ $`{\displaystyle \frac{i\mathrm{\Gamma }(3D/2)}{(4\pi )^{D/2}}}{\displaystyle _0^1}𝑑z_1{\displaystyle _0^{1z_1}}𝑑z_2(P^2+M^2)^{D/23}`$
where for the last step a standard integral was used which is calculated by employing the Wick rotation. Note in addition that a shift $`kk^{}=k+P`$ was performed. If we start with a vector integral including a factor $`k^\mu `$ in the numerator, one obtains

$`I^\mu `$ $`=`$ $`\mathrm{\Gamma }(3){\displaystyle _0^1}𝑑z_1{\displaystyle _0^{1z_1}}𝑑z_2{\displaystyle \frac{d^Dk^{}}{(2\pi )^D}\frac{(k^{}P)^\mu }{(k^2P^2M^2)^3}}=`$ (64)
$`=`$ $`{\displaystyle \frac{i\mathrm{\Gamma }(3D/2)}{(4\pi )^{D/2}}}{\displaystyle _0^1}𝑑z_1{\displaystyle _0^{1z_1}}𝑑z_2P^\mu (P^2+M^2)^{D/23},`$
the integrand proportional to $`k^\mu `$ vanishes because the function is odd in $`k^\mu `$. For the tensor integral we obtain

$$\stackrel{~}{I}^{\mu \nu }=\frac{i\mathrm{\Gamma }(3D/2)}{(4\pi )^{D/2}}_0^1dz_1_0^{1z_1}dz_2P^\mu P^\nu (P^2+M^2)^{D/23}+g^{\mu \nu }J=:I^{\mu \nu }+g^{\mu \nu }J.$$
(65)
### Combining both ingredients

If we combine the Dirac structure and the integrals, we obtain (up to general factors)

$$4(p+p^{})^\mu \gamma _\lambda I^\lambda 2(D2)\gamma _\rho I^{\mu \rho }2(D2)\mathrm{\Gamma }^\mu J+4mI^\mu .$$
(66)
Again, the term proportional to $`J`$ does not contribute, we do not have to care about this integral in the following. Including all general factors, we obtain

$`\mathrm{\Lambda }^\mu `$ $`=`$ $`{\displaystyle \frac{e^2\mathrm{\Gamma }(3D/2)}{(4\pi )^{D/2}}}{\displaystyle _0^1}dz_1{\displaystyle _0^{1z_1}}dz_2\times `$ (67)
$`\times \overline{u}(p^{})(4(p+p^{})^\mu P/2(D2)P^\mu P/4mP^\mu )u(p)(P^2+M^2)^{D/23}.`$
The integral is finite. Therefore, we can replace $`D=4`$. Now we use

$$\overline{u}(p^{})P/u(p)=\overline{u}(p^{})(z_1p/+z_2p/^{})u(p)=m(z_1+z_2)\overline{u}(p^{})u(p)$$
(68)
and write $`M=(1z_1z_2)s`$ and

$$P^\mu =z_1p^\mu +z_2p^\mu =\frac{1}{2}(z_1+z_2)(p+p^{})^\mu +\frac{1}{2}(z_1z_2)(pp^{})^\mu $$
(69)
where again only the first part counts for the anomalous magnetic moment. Finally, we can use $`z=z_1+z_2`$ as a new variable and obtain

$`\mathrm{\Lambda }^\mu `$ $`=`$ $`{\displaystyle \frac{e^2}{16\pi ^2}}\overline{u}(p^{})(p+p^{})^\mu u(p){\displaystyle _0^1}𝑑z_1{\displaystyle _{z_1}^1}𝑑z{\displaystyle \frac{4mz2mz^22mz}{m^2z^2+(1z)s}}=`$ (70)
$`=`$ $`{\displaystyle \frac{e^2}{8\pi ^2m}}\overline{u}(p^{})(p+p^{})^\mu u(p){\displaystyle _0^1}𝑑z{\displaystyle _0^z}𝑑z_1{\displaystyle \frac{z(1z)}{z^2+(1z)s/m^2}}=`$
$`=`$ $`{\displaystyle \frac{\alpha }{2\pi m}}\overline{u}(p^{})(p+p^{})^\mu u(p){\displaystyle _0^1}{\displaystyle \frac{z^2(1z)dz}{z^2+(1z)s/m^2}}.`$
At this point the function $`K(s)`$ appears.

### Explicit form of $`K(s)`$

The integral for $`K(s)`$ can be calculated. First we can write

$`K(s)`$ $`=`$ $`{\displaystyle _0^1}{\displaystyle \frac{z^2(1z)dz}{z^2+(1z)s/m^2}}={\displaystyle _0^1}{\displaystyle \frac{z(z^2+(1z)s/m^2)+(1z)zs/m^2+z^2}{z^2+(1z)s/m^2}}=`$
$`\mathrm{}`$ $`=`$ $`{\displaystyle _0^1}z𝑑z+\left(1{\displaystyle \frac{s}{m^2}}\right){\displaystyle _0^1}𝑑z+{\displaystyle _0^1}{\displaystyle \frac{zs/m^2(2s/m^2)s/m^2(1s/m^2)}{z^2+(1z)s/m^2}}𝑑z.`$ (71)
The denominator in the integrand can be written as product

$$z^2+(1z)\frac{s}{m^2}=(zz_+)(zz_{})$$
(72)
where

$$z_\pm =\frac{s}{2m^2}(1\pm v)=\frac{2}{1v},v=\sqrt{1\frac{4m^2}{s}}.$$
(73)
Therefore, we make the ansatz

$$\frac{A_+}{zz_+}+\frac{A_{}}{zz_{}}=\frac{(A_++A_{})z(A_+z_{}+A_{}z_+)}{(zz_+)(zz_{})}$$
(74)
and solve the system of equations

$$A_++A_{}=\frac{s}{m^2}\left(2\frac{s}{m^2}\right),A_+z_{}+A_{}z_+=\frac{s}{m^2}\left(1\frac{s}{m^2}\right)$$
(75)
to obtain finally

$$vA_+=\frac{(1+v)^2}{(1v)^2},vA_{}=+\frac{(1v)^2}{(1+v)^2}.$$
(76)
Using

$$1z_\pm =1\frac{2}{1v}=\frac{1\pm v}{1v}\frac{1z_\pm }{z_\pm }=\frac{1\pm v}{2},$$
(77)
we obtain

$`K(s)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{4}{1v^2}}+A_+\mathrm{ln}\left({\displaystyle \frac{1z_+}{z_+}}\right)+A_{}\mathrm{ln}\left({\displaystyle \frac{1v_{}}{v_{}}}\right)=`$ (78)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{4}{1v^2}}{\displaystyle \frac{(1+v)^2}{v(1v)^2}}\mathrm{ln}\left({\displaystyle \frac{1+v}{2}}\right)+{\displaystyle \frac{(1v)^2}{v(1+v)^2}}\mathrm{ln}\left({\displaystyle \frac{1v}{2}}\right).`$
In using $`x=(1v)/(1+v)`$, we finally obtain

$`K(s)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{(1+x)^2}{x}}{\displaystyle \frac{1+x}{x^2(1x)}}\mathrm{ln}\left({\displaystyle \frac{1}{1+x}}\right)+x^2{\displaystyle \frac{1+x}{1x}}\mathrm{ln}\left({\displaystyle \frac{x}{1+x}}\right)=`$ (79)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{(1+x)^2}{x}}+{\displaystyle \frac{1+x}{1x}}\left({\displaystyle \frac{1}{x^2}}x^2\right)\mathrm{ln}(1+x)+x^2{\displaystyle \frac{1+x}{1x}}\mathrm{ln}x`$
which can be shown to be identical with the expression in the main text.

## Appendix B Calculation details for the three models

The main point here is to calculate the Euclidean analogon to the function $`K(s)`$, i.e. $`W(t)`$. We start with a reformulation of the denominator in the integrand of $`K(s)`$,

$$z^2+(1z)\frac{s}{m^2}=\frac{1z}{m^2}\left(s+\frac{z^2}{1z}m^2\right)$$
(80)
and conclude that

$$t=\frac{z^2}{1z}m^2z=\frac{1}{2m^2}\left(t+\sqrt{t^2+4m^2t}\right)$$
(81)
is the substitution we have to use. We obtain

$`dz`$ $`=`$ $`{\displaystyle \frac{1}{2m^2}}\left(1+{\displaystyle \frac{t+2m^2}{\sqrt{t^2+4m^2t}}}\right)dt={\displaystyle \frac{t+2m^2\sqrt{t^2+4m^2t}}{2m^2\sqrt{t^2+4m^2t}}}dt,`$
$`z^2`$ $`=`$ $`{\displaystyle \frac{1}{4m^4}}\left(t^2+t^2+4m^2t2t\sqrt{t^2+4m^2t}\right)={\displaystyle \frac{t}{2m^4}}\left(t+2m^2\sqrt{t^2+4m^2t}\right),`$
$`1z`$ $`=`$ $`{\displaystyle \frac{1}{2m^2}}\left(t+2m^2\sqrt{t^2+4m^2t}\right)`$ (82)
and are lucky that a general factor occurs at all places. Moreover, for this general factor we can use

$$\left(t+2m^2\sqrt{t^2+4m^2t}\right)\left(t+2m^2+\sqrt{t^2+4m^2t}\right)=4m^4$$
(83)
to switch freely from numerator to denominator and vice versa. For the denominator we obtain

$$z^2+(1z)\frac{s}{m^2}=\frac{1}{2m^4}\left(t+2m^2\sqrt{t^2+4m^2t}\right)(s+t),$$
(84)
for the numerator we get

$$z^2(1z)dz=\frac{t(t+2m^2\sqrt{t^2+4m^2t})^3}{8m^8\sqrt{t^2+4m^2t}}dt,$$
(85)
therefore, finally

$`K(s)`$ $`=`$ $`{\displaystyle _0^1}{\displaystyle \frac{z^2(1z)dz}{z^2+(1z)s/m^2}}={\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{t(t+2m^2\sqrt{t^2+4m^2t})^2}{4m^4\sqrt{t^2+4m^2t}(s+t)}}𝑑t=`$ (86)
$`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{4m^4tdt}{\sqrt{t^2+4m^2t}(t+2m^2+\sqrt{t^2+4m^2t})^2(s+t)}}={\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{tW(t)dt}{s+t}}.`$
Finally, we have

$`a_\mu `$ $`=`$ $`4\pi ^2\left({\displaystyle \frac{\alpha }{\pi }}\right)^2{\displaystyle _{4m_\pi ^2}^{\mathrm{}}}{\displaystyle \frac{K(s)}{s}}\rho (s)𝑑s=4\pi ^2\left({\displaystyle \frac{\alpha }{\pi }}\right)^2{\displaystyle _{4m_\pi ^2}^{\mathrm{}}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{tW(t)}{s(s+t)}}\rho (s)𝑑t𝑑s=`$ (87)
$`=`$ $`4\pi ^2\left({\displaystyle \frac{\alpha }{\pi }}\right)^2{\displaystyle _0^{\mathrm{}}}W(t){\displaystyle _{4m_\pi ^2}^{\mathrm{}}}{\displaystyle \frac{t\rho (s)}{s(s+t)}}𝑑s𝑑t=4\pi ^2\left({\displaystyle \frac{\alpha }{\pi }}\right)^2{\displaystyle _0^{\mathrm{}}}W(t)\left(\mathrm{\Pi }(t)\right)𝑑t.`$
### The Euclidean weight $`F(t)`$

$`W(t)`$ was not the final weight function we selected for our considerations. Instead, we performed an integration-by-parts, taking the derivative of $`\mathrm{\Pi }(t)`$ and integrating $`W(t)`$. The Euclidean weight function $`F(t)`$ should be given by

$$F(t)=_t^{\mathrm{}}W(t^{})𝑑t^{}.$$
(88)
We show this by checking that $`F(\mathrm{})=0`$ and $`F^{}(t)=W(t)`$ for the fiven

$$F(t)=\frac{1}{2}\left(\frac{t+2m^2\sqrt{t^2+4m^2t}}{t+2m^2+\sqrt{t^2+4m^2t}}\right)=\frac{2m^4}{(t+2m^2+\sqrt{t^2+4m^2t})^2}.$$
(89)
The first condition is granted. For the second one we take

$`F^{}(t)`$ $`=`$ $`{\displaystyle \frac{4m^4(1+(t+2m^2)/\sqrt{t^2+4m^2t})}{(t+2m^2+\sqrt{t^2+4m^2t})^3}}={\displaystyle \frac{4m^4(t+2m^2+\sqrt{t^2+4m^2t})}{\sqrt{t^2+4m^2t}(t+2m^2+\sqrt{t^2+4m^2t})^3}}=`$ (90)
$`=`$ $`{\displaystyle \frac{4m^4}{\sqrt{t^2+4m^2t}(t+2m^2+\sqrt{t^2+4m^2t})^2}}=W(t).`$
### The correlator for model 1

The first calculation is quite easy and short,

$`{\displaystyle \frac{\mathrm{\Pi }_1(t)}{N_cQ_{\mathrm{eff}}^2}}`$ $`=`$ $`t{\displaystyle _{4m_\pi ^2}^{\mathrm{}}}{\displaystyle \frac{\rho _1(s)ds}{s(s+t)}}=t{\displaystyle _{4m_\pi ^2}^{\mathrm{}}}{\displaystyle \frac{\theta (s4m_{\mathrm{eff}}^2)}{s(s+t)}}=t{\displaystyle _{4m_{\mathrm{eff}}^2}^{\mathrm{}}}{\displaystyle \frac{ds}{s(s+t)}}=`$ (91)
$`=`$ $`{\displaystyle _{4m_{\mathrm{eff}}^2}^{\mathrm{}}}\left({\displaystyle \frac{ds}{s}}{\displaystyle \frac{ds}{s+t}}\right)=\mathrm{ln}\left({\displaystyle \frac{s}{s+t}}\right)|_{4m_{\mathrm{eff}}^2}^{\mathrm{}}=\mathrm{ln}\left({\displaystyle \frac{4m_{\mathrm{eff}}^2}{t+4m_{\mathrm{eff}}^2}}\right),`$
$`{\displaystyle \frac{d\mathrm{\Pi }_1(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{N_cQ_{\mathrm{eff}}^2}{t+4m_{\mathrm{eff}}^2}}.`$ (92)
### The correlator for model 2

Starting with

$$\rho (s)=N_cQ_{\mathrm{eff}}^2\sqrt{1\frac{4m^2}{s}}\left(1+\frac{2m^2}{s}\right)=\frac{1}{2}N_cQ_{\mathrm{eff}}^2v(3v^2),v=\sqrt{1\frac{4m^2}{s}},s=\frac{4m^2}{1v^2}$$
(93)
we obtain

$`{\displaystyle \frac{\mathrm{\Pi }_2(t)}{N_cQ_{\mathrm{eff}}^2}}`$ $`=`$ $`t{\displaystyle _{4m^2}^{\mathrm{}}}{\displaystyle \frac{ds}{s(s+t)}}\sqrt{1{\displaystyle \frac{4m^2}{s}}}\left(1+{\displaystyle \frac{2m^2}{s}}\right)=t{\displaystyle _0^1}{\displaystyle \frac{v^2(3v^2)dv}{4m^2+t(1v^2)}}=`$ (94)
$`=`$ $`t{\displaystyle _0^1}{\displaystyle \frac{v^2(3v^2)dv}{4m^2+ttv^2}}={\displaystyle _0^1}{\displaystyle \frac{v^2(3v^2)dv}{a^2v^2}},a^2=1+{\displaystyle \frac{4m^2}{t}}.`$
We continue with

$`{\displaystyle \frac{\mathrm{\Pi }_2(t)}{N_cQ_{\mathrm{eff}}^2}}`$ $`=`$ $`{\displaystyle _0^1}{\displaystyle \frac{v^2(3v^2)dv}{a^2v^2}}={\displaystyle _0^1}v^2𝑑v(a^23){\displaystyle _0^1}𝑑v+a^2(a^23){\displaystyle _0^1}{\displaystyle \frac{dv}{a^2v^2}}=`$
$`=`$ $`{\displaystyle \frac{1}{3}}(a^23)+(a^23){\displaystyle _0^1}{\displaystyle \frac{dv}{1v^2/a^2}}={\displaystyle \frac{1}{3}}+(a^23)\left(artanh\left({\displaystyle \frac{1}{a}}\right)1\right).`$


text:
# Entanglement capability of self-inverse Hamiltonian evolution

## Abstract

We determine the entanglement capability of self-inverse Hamiltonian evolution, which reduces to the known result for Ising Hamiltonian, and identify optimal input states for yielding the maximal entanglement rate. We introduce the concept of the operator entanglement rate, and find that the maximal operator entanglement rate gives a lower bound on the entanglement capability of a general Hamiltonian.

Two-qubit unitary gates occupy a central role in quantum information science Nie00 , particularly because of the capability of these gates to generate and enhance entanglement of the state of the system. In general, entanglement capability can be enhanced by introducing ancillary states Dur01 ; Kra01 , but for the important and interesting case of the Ising Hamiltonian $`H_{\text{ISING}}=\sigma _z\sigma _z`$ it is ancilla-independent Chi02 ; Childs  (where $`\sigma _z`$=diag$`(1,1)`$ is a Pauli matrix). The independence of the entanglement capability on ancillas is a consequence of the self-inverse property $`H_{\text{ISING}}=H_{\text{ISING}}^1`$.

We generalize this result to all Hamiltonian evolution of the type

$$H=X_AX_B$$
(1)
such that $`X_i=X_i^1_i`$ for $`i\{A,B\}`$ and $`H=H^1`$, with the Ising Hamiltonian being a special case. Here, we assume that $`X_i`$ is not an identity operator $`𝟙_𝕚`$. Due the self-inverse property of the Hamiltonian, we have the evolution operator ($`\mathrm{}=1`$)

$$U(t)=e^{iHt}=\mathrm{cos}t𝟙_𝔸𝟙_𝔹𝕚\mathrm{sin}𝕥𝕏_𝔸𝕏_𝔹.$$
(2)
We employ operator entanglement Zan01 ; Nie02  and operator entanglement rate to characterize the entanglement capability of the self-inverse Hamiltonian evolution: these approaches yield elegant analyses that are simple to apply for determining entanglement capabilities. We find that the maximal operator entanglement rate always gives a lower bound on the entanglement capability of a general Hamiltonian.

The entanglement capability of a Hamiltonian Dur01 ; Chi02 ; Kra02  is defined relative to a specified entanglement measure. We use von Neumann entropy as our entanglement measure of a pure state $`|\mathrm{\Psi }_{AB}`$:

$$E(|\mathrm{\Psi })=\text{Tr}_A(\rho _A\mathrm{log}\rho _A),$$
(3)
where $`\rho _A`$=$`\text{Tr}_B(\rho _{AB})`$, $`\rho _{AB}=|\mathrm{\Psi }\mathrm{\Psi }|`$, and $`\mathrm{log}`$ is always base 2. The entanglement capability of Hamiltonian $`H`$ is defined as the maximum entanglement rate when a pure state is acted on by the associated evolution operator $`U(t)=\mathrm{exp}(iHt)`$. Mathematically, the ancilla-unassisted and ancilla-assisted entanglement capabilities are defined as Dur01 ; Chi02

$`E_H:`$ $`=\text{max}_{|\mathrm{\Psi }_{AB}}\mathrm{\Gamma }(t)|_{t0},`$ (4)
$`E_H^{\text{anc}}:`$ $`=\text{sup}_{|\mathrm{\Psi }_{A^{}ABB^{}}}\mathrm{\Gamma }(t)|_{t0},`$ (5)
respectively, where $`\mathrm{\Gamma }(t)=\text{d}E[U(t)|\mathrm{\Psi }(0)]/\text{d}t`$ is the state entanglement rate, and ancilla systems $`A^{}`$ and $`B^{}`$ are not acted on by Hamiltonian $`H`$. For the ancilla-assisted case the entanglement refers to the bipartite entanglement between systems $`A^{}A`$ and $`BB^{}`$.

Let us first consider the unassisted case, and the initial pure state of the system can always be Schmidt-decomposed as Peres

$$|\mathrm{\Psi }(0)=\underset{n}{}\sqrt{\lambda _n}|\psi _n|\varphi _n,$$
(6)
where $`\{|\psi _n\}`$ and $`\{|\varphi _n\}`$ are orthonormal sets of states, and $`\lambda _n>0n`$. As $`|\mathrm{\Psi }(0)`$ is normalized, $`_n\lambda _n=1`$. The state at time $`t`$ is described by the density operator $`\rho _{AB}(t)=U(t)\rho _{AB}(0)U^{}(t)`$, which satisfies $`\dot{\rho }_{AB}(t)=i[H,\rho _{AB}(t)]`$. Thus, the reduced density operator $`\rho _A(t)`$ satisfies

$$\dot{\rho }_A(t)=i\text{Tr}_B[H,\rho _{AB}(t)].$$
(7)
From Eq. (3), we know that the entanglement rate is Chi02

$$\mathrm{\Gamma }(t)=\text{Tr}_A[\dot{\rho }_A(t)\mathrm{log}\rho _A(t)].$$
(8)
Substituting Eq. (7) into Eq. (8), we obtain the entanglement rate at $`t=0`$ as

$$\mathrm{\Gamma }(t)|_{t0}=i\text{Tr}_A\{\text{Tr}_B[H,\rho _{AB}(0)]\mathrm{log}\rho _A(0)\},$$
(9)
which is general for arbitrary Hamiltonians and initial states. By varying initial states we can maximize the entanglement rate.

Equation (9) is solved by first obtaining the results

$$\text{Tr}_B[H,\rho _{AB}(0)]=\underset{mn}{}\sqrt{\lambda _m\lambda _n}(X_B)_{mn}[X_A,|\psi _n\psi _m|],$$
(10)
where $`(X_B)_{mn}=\varphi _m|X_B|\varphi _n`$, and

$$\mathrm{log}\rho _A(0)=\underset{n}{}\mathrm{log}\lambda _n|\psi _n\psi _n|.$$
(11)
The entanglement rate at $`t=0`$ for our Hamiltonian $`H`$ is then given by,

$$\mathrm{\Gamma }(t)|_{t0}=i\underset{mn}{}\sqrt{\lambda _m\lambda _n}\mathrm{log}\frac{\lambda _m}{\lambda _n}(X_A)_{mn}(X_B)_{mn},$$
(12)
where $`(X_A)_{mn}=\psi _m|X_A|\psi _n`$.

From Eq. (12), we obtain an upper bound for the entanglement rate

$`\mathrm{\Gamma }(t)|_{t0}`$ $`{\displaystyle \underset{mn}{}}\sqrt{\lambda _m\lambda _n}\left|\mathrm{log}{\displaystyle \frac{\lambda _m}{\lambda _n}}\right||(X_A)_{mn}||(X_B)_{mn}|`$
$`=`$ $`{\displaystyle \underset{mn}{}}(\lambda _m+\lambda _n)\sqrt{{\displaystyle \frac{\lambda _m}{\lambda _m+\lambda _n}}{\displaystyle \frac{\lambda _n}{\lambda _m+\lambda _n}}}\left|\mathrm{log}{\displaystyle \frac{\lambda _m}{\lambda _n}}\right|`$
$`\times |(X_A)_{mn}||(X_B)_{mn}|`$
$``$ $`{\displaystyle \frac{\beta }{2}}{\displaystyle \underset{mn}{}}(\lambda _m+\lambda _n)|(X_A)_{mn}||(X_B)_{mn}|`$
$``$ $`{\displaystyle \frac{\beta }{4}}{\displaystyle \underset{mn}{}}(\lambda _m+\lambda _n)(|(X_A)_{mn}|^2+|(X_B)_{mn}|^2)`$
$`=`$ $`\beta ,`$ (13)
where

$`\beta =`$ $`2\underset{x[0,1]}{\mathrm{max}}\sqrt{x(1x)}\mathrm{log}[x/(1x)]`$
$`=`$ $`2\sqrt{x_0(1x_0)}\mathrm{log}[x_0/(1x_0)]1.9123,`$ (14)
and $`x_0=0.9128`$. The first line simply follows from the triangle inequality, the second and third lines from Ref. Chi02 , and the fourth line from the inequality $`2|ab||a|^2+|b|^2`$ for any complex numbers $`a`$ and $`b`$. Finally, the last line results from the self-inverse properties of $`X_A`$ and $`X_B`$. To see this, let us examine the sum $`_{mn}\lambda _m|(X_A)_{mn}|^2=_{mn}\lambda _m(X_A)_{mn}(X_A)_{mn}^{}=_m\lambda _m_n(X_A)_{mn}(X_A)_{nm}=1`$, where the second equality results from the Hermitian property of $`X_A`$, and the last equality from the self-inverse property $`X_A=(X_A)^1`$. We have seen that the self-inverse property is essential to obtain the upper bound $`\beta `$. Another feature of Eq. (13) is that the result is applicable for any pure state with or without ancillas. Therefore, the upper bound $`\beta `$ is ancilla-independent. Next we show that the upper bound can be saturated by optimal input states.

Since the self-inverse operators $`X_i,i\{A,B\}`$ satisfy $`X_i^2=1`$, the eigenvalues of $`X_i`$ are $`\pm 1`$ with the corresponding eigenstates (possibly degenerate) denoted by $`|\pm _i`$. Then, we construct the optimal input states given by Referee

$`|\mathrm{\Psi }(0)={\displaystyle \frac{\sqrt{x_0}}{2}}(|+_A+|_A)(|+_B+|_B)`$
$`+{\displaystyle \frac{i\sqrt{1x_0}}{2}}(|+_A|_A)(|+_B|_B).`$ (15)
This state is of Schmidt form with Schmidt number 2. From Eqs. (12) and (14), we obtain $`\mathrm{\Gamma }(t)|_{t0}=\beta `$. Therefore, we can always find the optimal input states for which

$$E_H=E_H^{\text{anc}}=\beta .$$
(16)
Result (16) extends results for the Ising Hamiltonian Dur01 ; Chi02 , and we do not restrict the dimension of Hilbert spaces on which the self-inverse operators $`X_A`$ and $`X_B`$ act. The Hilbert space can be finite-dimensional or infinite-dimensional.

Self-inverse Hamiltonians not only exhibit elegant entanglement capability, but are also physically meaningful. As an example of a physical system with optimal input state and maximal entanglement rate $`\beta `$, we consider a spin-$`j`$ system with a $`(2j+1)`$-dimensional space spanned by basis states $`\{|n_j|j;nj,n=0,1,\mathrm{},2j\}`$. We define a number operator $`𝒩=J_z+j`$ of the spin-$`j`$ system satisfying $`𝒩|n_j=n|n_j`$. Here, commutators for $`J_z`$ and ladder operators $`J_\pm `$ satisfy the su(2) algebra $`[J_z,J_\pm ]=\pm J_\pm ,[J_+,J_{}]=2J_z`$. Then, using the number operators, we construct the self-inverse interaction Hamiltonian

$$H_1=(1)^𝒩(1)^𝒩,$$
(17)
where $`(1)^𝒩`$ is a parity operator of the system. This kind of Hamiltonian is realizable in physical systems Geometric . From Eq. (17), we see that Hamiltonian $`H_1`$ reduces to the Ising Hamiltonian $`H_{\text{ISING}}`$ for the case of $`j=1/2`$. When the dimension $`d=2j+1`$ is even, i.e., $`j`$ is a half integer, Hamiltonian $`H_1`$ is equivalent to the Ising Hamiltonian in the sense that the operator $`(1)^𝒩`$ can be written as $`𝟙_{𝕕/\mathrm{𝟚}}\sigma _𝕫`$, where $`𝟙_{}`$ denotes the $`N\times N`$ identity matrix. For integer $`j`$, Hamiltonian $`H_1`$ is not equivalent to the Ising Hamiltonian.

Let us consider the SU(2) entangled coherent state (ECS) ECS  $`|\eta _{\text{ECS}}`$ as an input state

$$|\eta _{\text{ECS}}=\sqrt{x_0}|\eta |\eta +i\sqrt{1x_0}|\eta |\eta ,$$
(18)
where $`|\eta =\mathrm{exp}(\eta J_+\eta ^{}J_{})|0_j`$ is the SU(2) or spin coherent state (SCS) SCS , $`\eta `$ is complex with unit modulus. From Eq. (12) and the identity $`(1)^𝒩|\eta =|\eta `$, it is straightforward to show that the SU(2) ECS is the optimal input state which yields the maximal entanglement rate $`\beta `$.

The above analysis is restricted to finite equal-dimensional composite systems, but this limitation is convenient, not essential. Consider the harmonic oscillator with infinite-dimensional Hilbert space. The operator $`a^{}a`$ is the number operator, where $`a^{}`$ ($`a`$) is the creation (annihilation) operator. We can have the Hamiltonians: $`H_2=(1)^𝒩(1)^{a^{}a}`$ and $`H_3=(1)^{a^{}a}(1)^{a^{}a}`$. By noticing that the SCS can be realized in the Fock space as the binomial state BS , the optimal input states for Hamiltonians $`H_2`$ and $`H_3`$ are directly obtained by appropriately replacing SCS with the binomial state in the SU(2) ECS (18).

Although in general one cannot expect an entangled state to be generated by Hamiltonian evolution from a product state, in this spin system the optimal input state $`|\eta _{\text{ECS}}`$ can be generated by the Hamiltonian $`H_1`$ from the product state $`|\eta |\eta `$. Mathematically, $`|\eta _{\text{ECS}}=\mathrm{exp}(iH_1t)|\eta |\eta `$ with $`x_0=\mathrm{cos}^2t`$. This fact and the ancilla-independence of the self-inverse Hamiltonian entanglement capability suggest that the entanglement capability depends more on the Hamiltonian and the associated unitary operator. Next, we demonstrate this character by introducing the entanglement rate for operator $`U(t)`$, and give the result that the entanglement capability of Hamiltonian $`H_1`$ is equal to the maximal operator entanglement rate for the associated unitary operator.

We review operator entanglement introduced in Ref. Zan01 . Notice that the linear operators over $`_d`$ (finite $`d`$) also form a $`d^2`$-dimensional Hilbert space denoted by $`_{d^2}^{\text{HS}}`$, and the corresponding scalar product between two operators $`X`$ and $`Y`$ is given by the Hilbert-Schmidt product $`X,Y:=\text{Tr}(X^{}Y)`$, and $`X_{\text{HS}}=X,X.`$ Then, the operator acting on $`_{d_1}_{d_2}`$can be considered as a state on the composite Hilbert space $`_{d_1^2}^{\text{HS}}_{d_2^2}^{\text{HS}}`$, and the operator entanglement can be defined as the entanglement of that state.

Any unitary operator $`V`$ acting on $`_{d_1}_{d_2}`$ may be Schmidt-decomposed as Nie02  $`V=_ns_nA_nB_n`$, where $`s_n>0n`$ and $`A_n`$ and $`B_n`$ are orthonormal operator bases for systems $`1`$ and $`2.`$ From the Schmidt form the entanglement of a unitary operator $`V`$ is determined to be $`(V)=_ls_l^2/(d_1d_2)\mathrm{log}\left[s_l^2/(d_1d_2)\right],`$ where the factor $`1/(d_1d_2)`$ arises from normalization of the unitary operator. We can think of the operator entanglement as a strength measure of the operator Nie02 .

If we consider the unitary operator $`U(t)`$, the operator entanglement $`[U(t)]`$ becomes a time-dependent function. Analogous to the definition of the state entanglement rate Dur01 , it is natural to define the operator entanglement rate at a certain time of interaction and the maximal entanglement rate

$$R(t):=\text{d}[U(t)]/\text{d}t,R_{\mathrm{max}}:=\underset{t}{\mathrm{max}}R(t),$$
(19)
respectively, where the maximization is over all time $`t`$. At certain times the entanglement rate becomes maximal.

The operator entanglement of an arbitrary unitary operator $`V`$ is equal to the entanglement of the state $`V(|\mathrm{\Phi }_{A^{}A}|\mathrm{\Phi }_{BB^{}})`$ Nie02 ,

$$(V)=E[V(|\mathrm{\Phi }_{A^{}A}|\mathrm{\Phi }_{BB^{}})],$$
(20)
where

$`|\mathrm{\Phi }_{A^{}A}`$ $`={\displaystyle \underset{n=0}{\overset{d_11}{}}}|n_A^{}|n_A_{d_1}_{d_1},`$ (21)
$`|\mathrm{\Phi }_{BB^{}}`$ $`={\displaystyle \underset{n=0}{\overset{d_21}{}}}|n_B|n_B^{}_{d_2}_{d_2}`$ (22)
are maximally entangled states. Equation (20) shows a direct relation between operator entanglement and state entanglement. Note that here dimension $`d_1`$ can differ from dimension $`d_2`$. From Eq. (20), we immediately have a relation between the operator entanglement rate and the state entanglement rate

$$R(t)=\mathrm{\Gamma }(t)=\dot{E}[U(t)(|\mathrm{\Phi }_{A^{}A}|\mathrm{\Phi }_{BB^{}})].$$
(23)
From Eq. (23) we find that the maximal operator entanglement rate gives a lower bound for the entanglement capability of a general Hamiltonian, which is the infinitesimal version of the one that the operator entanglement gives a lower bound of the entanglement capability of a unitary operator Nie02 .

Having defined the operator entanglement rate let us study the maximal operator entanglement rate $`R_{\mathrm{max}}`$ of the evolution operator associated with Hamiltonian $`H_1`$ with dimension $`d=2j+1`$ an even number Even . The associated unitary operator is in the Schmidt form

$$U_1(t)=\mathrm{cos}t𝟙_𝕕𝟙_𝕕𝕚\mathrm{sin}𝕥(\mathrm{𝟙})^𝒩(\mathrm{𝟙})^𝒩.$$
(24)
From the above equation, the operator entanglement and operator entanglement rate are given by

$`[U_1(t)]=`$ $`\mathrm{cos}^2t\mathrm{log}\mathrm{cos}^2t\mathrm{sin}^2t\mathrm{log}\mathrm{sin}^2t,`$ (25)
$`R(t)=`$ $`R[U_1(t)]=\mathrm{sin}(2t)\mathrm{log}(\mathrm{cot}^2t),`$ (26)
respectively. The above equation shows that the operator entanglement rate is a periodic function of time $`t`$; hence we can maximize over one period. It is straightforward to find that at time $`t=0.2932`$, the entanglement rate reaches its maximum value, i.e.,

$$R_{\mathrm{max}}=\beta .$$
(27)
Comparing Eqs. (16) and (27), we conclude that the entanglement capability of Hamiltonian $`H_1`$ is equal to the maximum operator entanglement rate of the associated unitary operator $`\mathrm{exp}(iH_1t)`$. Therefore, the entanglement capability of the self-inverse Hamiltonian is inherent in the evolution operator in the sense that it can be solely determined by the evolution operator, irrespective of states.

More generally, there is interest in the entanglement capability of quantum operators, not just Hamiltonians Mak00 ; Zan00 ; Kra01 ; Lei02 ; Ben02 ; Ber02 ; Wol02 ; WanSan02 . Here, we analyze the entanglement capability of the unitary operator generated by the self-inverse Hamiltonian $`H`$ (2). We quantify the entanglement capability of a unitary operator $`U`$ by the maximum entanglement which a unitary operator can create given an initial product state Kra01 :

$$E_U=\text{sup}_{|\gamma ,|\delta }E(U|\gamma |\delta ),$$
(28)
where $`|\gamma _{A^{}A}`$ and $`|\delta _{BB^{}}`$, namely, we include ancillas. For unitary operator $`U(t)`$ we immediately obtain

$`E_{U(t)}=`$ $`\text{sup}_{|\gamma ,|\delta }E[U(t)|\gamma |\delta ]=\text{sup}_{|\gamma ,|\delta }E[|\mathrm{\Psi }(t)],`$ (29)
$`|\mathrm{\Psi }(t)=`$ $`\mathrm{cos}t|\gamma |\delta i\mathrm{sin}t|\overline{\gamma }|\overline{\delta },`$ (30)
with normalized states $`|\overline{\gamma }=X_A|\gamma `$ and $`|\overline{\delta }=X_B|\delta `$.

The state (30) belongs to a class of bipartite entangled states discussed in Ref. Non . To quantify the entanglement of this state, we can use the standard entanglement measure, the entropy of entanglement. However, for convenience, we use the concept of the concurrence $`C`$ Conc  to quantify the entanglement since the state can be viewed as a state of two pseudo-qubits Wangjpa . Using the results of Ref. Wangjpa , the concurrence of state $`|\mathrm{\Psi }(t)`$ is obtained as

$$C[|\mathrm{\Psi }(t)]=|\mathrm{sin}(2t)|\sqrt{(1|\gamma |\overline{\gamma }|^2)(1|\delta |\overline{\delta }|^2)}$$
(31)
From the above equation, it is easy to find that $`E_{U(t)}=|\mathrm{sin}(2t)|,`$ which is the maximal entanglement which the operator $`U(t)`$ can generate from an arbitrary product state. As the states $`|\gamma _{A^{}A}`$ and $`|\delta _{BB^{}}`$, the entanglement capability $`E_{U(t)}`$ is ancilla-independent.

The unitary operaror $`U_1(t)`$ (24) is in fact a special case of $`U(t)`$ (2). Assuming that $`U_1(t)`$ act on $`_d_d`$ with even $`d`$, we find that the optimal input state for generating maximal entanglement is $`|\eta |\eta `$, and from Eq. (24) we find that

$$[U_1(t)]=C[U_1(t)]=E_{U_1(t)}=|\mathrm{sin}(2t)|.$$
(32)
Therefore, the operator entanglement of $`U_1`$ is equal to the entanglement capability of the unitary operator.

In conclusion, we have determined the entanglement capability of self-inverse Hamiltonian evolution. The self-inverse Hamiltonians studied here are physically relevant in their own right, and reduce to the important Ising Hamiltonian that is central to quantum information science. These Hamiltonians go beyond the two-qubit cases; namely, they can act on $`d_1\times d_2`$ composite system, where the dimension of the subsystem can be either finite or infinite. We introduced the concept of operator entanglement rate which is well-defined for composite finite systems, and for certain self-inverse Hamiltonians the maximal operator entanglement rate is equal to the entanglement capability. For a general Hamiltonian the maximal operator entanglement rate always gives a lower bound on the Hamiltonian entanglement capability.

Acknowledgments: The authors acknowledge helpful discussions with Paolo Zanardi, Dominic Berry, Stephen Bartlett, and Gerard Milburn. This project has been supported by an Australian Research Council Large Grant.