QED Interference in Charge Asymmetry Near the Z Resonance at Future Electron-Positron Colliders
Stanislaw Jadach Scott Yost
hep-ph








IFJPAN-IV-2017-11








    QED Interference in Charge Asymmetry

      Near the Z Resonance at Future
      Electron-Positron Colliders ^⋆ 




5mm


S. Jadach^a, S. Yost^b,a

2mm
^aInstitute of Nuclear Physics, Polish Academy of Sciences,

  ul. Radzikowskiego 152, 31-342 Kraków, Poland



^bThe Citadel, Charleston, SC, USA
 






The measurement of the charge asymmetry 
A_ FB(e^-e^+→μ^-μ^+) will play an important role
at the high-luminosity circular electron-positron collider
FCCee considered for construction at CERN.  In particular,
near the Z resonance, √(s)≃ M_Z ± 3.5 GeV,
A_ FB will provide a very precise value of
the pure electromagnetic coupling constant α_ QED(M_Z),
which is vitally important for overall tests of the Standard Model.
For this purpose, A_ FB will be measured at the FCCee with 
an experimental error better than
δ A_ FB≃ 3· 10^-5, at least
a factor of 100 more precisely than at past LEP experiments!
The important question is whether the effect of interference between photon
emission in the initial and final state (IFI) can be removed from 
the A_ FB data 
at the same precision level using perturbative QED calculations.  A first 
quantitative study of this problem is presented here, with the help of the
 program and a newly developed calculation based on soft photon 
resummation, matched with NLO and NNLO fixed-order calculations.
It is concluded that a factor of 10 improvement with respect to the LEP era 
is obtained.
We also present a clear indication that reducing the uncertainty
of charge asymmetry near the Z peak due to IFI down to
δ A_ FB≃ 3· 10^-5, 
i.e. the expected experimental precision at FCCee, is feasible.





^⋆This work is partly supported by

 the Polish National Science Center grant 2016/23/B/ST2/03927,
 the Citadel Foundation 
 and the CERN FCC Design Study Programme.














§ INTRODUCTION



At the future high-energy high-luminosity circular electron-positron collider
FCCee <cit.>
proposed for construction at CERN, 
the measurement of the muon charge asymmetry
A^μ_ FB=A_ FB(e^-e^+→μ^-μ^+) 
will play an additional important role.
Near the Z resonance, √(s_±)≃ M_Z ± 3.5 GeV, 
the measurement of A^μ_ FB may provide a very precise value of
the pure electromagnetic coupling constant α_ QED(M_Z),
which is vitally important for overall tests of the Standard Model (SM),
especially for the SM prediction of sin^2θ^l_ eff and 
M_W <cit.>, 
at a precision level at least a factor of 10 better than presently.
This kind of the measurement of α_ QED(M_Z) was proposed and 
analyzed in Ref. <cit.>.

In past measurements of the charge asymmetry A^μ_ FB (forward-backward
angular asymmetry) at electron-positron colliders LEP and SLC, 
the QED interference between photons emitted from the initial and final charged 
leptons had to be taken into account.
Let us abbreviate in the following A^μ_ FB≡.  In overall tests of the SM, 
the measurement of  contributed mainly to knowledge of the Z 
couplings and/or the effective electroweak mixing angle sin^2θ_ eff^l
<cit.>.

Thanks to the very high luminosity of the FCCee <cit.>, 
the charge asymmetry 
(M_Z ± 3.5  GeV) will possibly be measured with an error
δ/ ≃ 3· 10^-5 or even better 
<cit.>.
This immediately poses the question of whether the effect of
QED initial-final interference can be removed from the data
at the same precision level.  
How big is the effect of IFI in ?
Far from the resonance, it is about 2-3%
and it is even bigger for a tight cutoff on the total energy of the 
emitted photons.  At the top of the Z resonance, 
where  was measured most precisely in the past LEP experiments, 
the IFI effect is suppressed by the ratio Γ_Z/M_Z 
to the level of δ∼ 0.1%,
due to the long time separation between the creation
and the decay stages of the Z resonance, 
as elaborated in many LEP era works, see
Refs. <cit.>.
As we shall see in our analysis, at √(s)≃ M_Z ± 3.5 GeV,
the  same Γ_Z/M_Z suppression of IFI in A_ FB
still works to some extent,
but the IFI effect is nevertheless at the δ A_ FB∼ 1% level,
and growing for tight cutoffs, in spite of partial cancellations
in the difference between values at √(s)≃ M_Z ± 3.5,
as already noticed in Ref <cit.>.

This effect is huge with respect to the planned experimental precision at 
FCCee, and it would render measurement of the A_ FB (M_Z ± 3.5  GeV)
completely useless unless the theoretical evaluation of IFI is equally precise!
Note that in the LEP data analysis, a cutoff on the total
photon energy was imposed by requiring a minimum value of the effective mass of
the muon pair, M_μ^-μ^+,  or a maximum acollinearity angle, 
which was translated into an upper limit on the total photon energy
E_γ^ tot/s^1/2 varying between 0.5 to 0.998
(see Table 2.1 in Ref. <cit.>).
Due to the higher precision of the FCCee, a stronger cutoff probably will
be preferred
in order to minimize the background from hadronic and tau pair channels,
and for better control of the angular dependence of the muon detector 
efficiency[P. Janot, private communication.].
Also, theoretical control over IFI may be better for a stronger cutoff
(in spite of its larger size) thanks to the power of soft photon resummation
and a better elimination of the four fermion contributions.
For this reason, we will use a photon energy cutoff 
E_γ^ tot/s^1/2≤ 0.2 or stronger.

How precise are the theoretical evaluations 
of the effect of IFI in  presently available in perturbative QED?
In the pre-LEP era, O(α^1) fixed-order calculations
were quoted to provide ∼ 0.3-0.5% precision; 
see the review of Ref. <cit.>.
In the LEP1 phase near the Z resonance, thanks to  Γ_Z/M_Z suppression, 
the IFI effect in  at the Z peak was not a burning issue.
For instance, in the work of Ref. <cit.> used in
the final data analysis of LEP1 of Ref. <cit.>,
the calculations of the IFI effect were done using
ZFITTER <cit.> 
and TOPAZ0 <cit.> programs, 
crosschecking them with the KORALZ Monte Carlo 
<cit.>.

In all these calculations and programs, the treatment of IFI
was at the O(α^1) fixed-order level,
without soft photon resummation.
Pioneering work on the resummation of soft photon effects 
near a narrow resonance, including resummation of ln(Γ_Z/M_Z), 
was already done earlier by the Frascati group;
see Refs. <cit.>,
but it was not exploited in the above studies, 
mainly because they did not include hard photon effects in a realistic way.

Significant progress on the IFI issue was made just before the end of the LEP 
era, with the advent of new method of the soft photon resummation
matched with fixed-order QED corrections up to α^2
and electroweak (EW) corrections up to α^1,
the so-called coherent exclusive exponentiation <cit.> (CEEX)
and its implementation in the  program <cit.>.
The CEEX implementation in  has
included all the advances of soft photon resummation of the IFI contributions
of Refs. <cit.>
relevant for narrow resonances.
[In particular, resummation of ln(Γ_Z/M_Z) was included.]
The SM predictions of  for  and other experimental 
observables were possible for arbitrary event selections (cuts),
because  is a regular Monte Carlo (MC) event generator.
Correct matching of the O(α^1) IFI contributions
with other non-IFI corrections up to complete O(α^2) QED
and  O(α^1) electroweak
was implemented throughout the entire multiphoton phase space, 
including any number of soft and  hard photons.  

The CEEX/ calculation was instrumental in the analysis
of LEP2 data above the Z peak and near the WW threshold,
and helped to consolidate data analysis of e^-e^+→ ff processes
near the Z peak.
The precision of the IFI calculations quoted at the end of the LEP era
was δ≃ 0.1% at the Z peak
and δ≃ 0.3% far away from the Z resonance;
see Refs. <cit.>.
These papers represent the state of the art in the perturbative QED 
calculation of the IFI contributions to  until the present day.

The  precision tag on the IFI calculations, both near the Z peak and 
away from it,
was more than sufficient for analyzing all LEP experimental data at the end of 
the LEP era.  However, this precision was quite clearly underestimated,
i.e. most likely it was far better.
However, it was difficult to better quantify the theoretical 
uncertainty of the IFI prediction of , because there was 
no other calculation at a similar level of sophistication to compare with.
One of the main aims of this work will be to develop
a new alternative numerical calculation of the IFI contribution,
in order to compare with  and quantify
theoretical uncertainties of the IFI component in 
at the precision level higher than presently available.

Generally, one may be quite skeptical whether an improvement 
of the QED calculation of IFI
in  from the LEP-era δ∼ 10^-3
down to δ∼ 10^-5, i.e. by a factor of 100,
is feasible at all!
However, there is an interesting precedence –
the prediction of perturbative QED for the Z line shape
also progressed by a similar factor from the time before LEP started
until the end of the LEP era.
This was possible mainly due to soft photon resummation techniques.
The use of these techniques is again critical for
the present task of improving the QED calculation of IFI in .
The aim of this paper is to check how far we can advance on the road
to the precision required for FCCee.

Let us stress, that present work is not a progress report on the development
of , simply because  remains the same as in 1999.
It is, however, definitely a progress report on the understanding of the IFI
contribution, thanks to several newly developed analytic calculations 
implemented in the new computer code  and a wealth of numerical results
for various kinds of matrix elements, phase space integrations, cutoffs, etc.
This work will have to be continued in the future,
including possible upgrades of the matrix element in ,
and/or development of new MC programs even more sophisticated than .

Finally, in view of the growing interest in the higher order SM calculations
which would match the very high precision
of the FCCee experiment <cit.>,
it is important to note that the CEEX methodology of photon resummation 
and matching with fixed-order nonsoft QED and EW corrections also addresses
some important issues in the QED+EW perturbative calculations,
beyond what was typically done for the LEP data analysis, as quoted above.
The basic issue is that of the separation of QED and pure EW parts of the SM
in the perturbative expansion.
This is necessary, because QED corrections are larger, and 
their soft part has to be resummed to infinite order while the
nonsoft part must be included up to α^4,
while the perturbative series for more complicated EW corrections can be 
truncated earlier,
at α^2 or α^3 <cit.>.

In the calculations for LEP data analysis (see <cit.> 
and other Refs. quoted above),
the issue of separating QED and EW parts was not a critical one,
because resummed higher-order QED was typically confined to the initial-state radiation (ISR) effective 
radiator function, and in the remaining α^1 calculations, the
 QED and EW parts enter additively, and thus are well separated 
(except of negligible IFI which was controlled up to α^1.
Beyond α^1, the QED and EW parts often enter multiplicatively,
for instance in 2-loop graphs with one loop involving photon exchange 
and another loop with massive bosons or fermions,
or one loop of pure EW origin with a hard photon emission insertion.
The CEEX technique provides for clear methodology of separating/factorizing
QED and EW parts at any order.

Omitting algebra which can be found in Refs. <cit.>,
the main points of CEEX methodology can be summarized as follows:


  (i) In the first step of the factorization of the infrared (IR) factors
at the amplitude level,
for any group of multiloop graphs with one photon insertion,
the IR part is subtracted at the amplitude level using a well-known (1-loop) 
function defined in the classic Yennie-Frautschi-Suura 
work <cit.> times a finite contribution one order lower without 
a photon insertion.  (A similar procedure applies for multiloop corrections 
with two and more photon insertions.)
The remaining finite non-IR remnant will be used in the next step.
Similarly, for any group of real photon insertions into 
a given multiloop diagram with pure EW content, one subtracts, at the 
amplitude level, a well-known eikonal factor times a basic diagram with 
pure EW content.

A similar well-defined procedure applies for amplitudes with more 
real and virtual photon insertions.
The first step is finalized by constructing spin amplitudes for an arbitrary 
number of real photons distributed over the entire phase space in which 
non-IR remnants after IR subtractions are reinserted in a well-defined way, 
while IR virtual factors are exponentiated and explicit IR-divergent eikonal 
factors are ready for MC integration in the next step
[Collinear contributions giving rise to nonsoft mass logarithms
are included order-by-order in the present version of CEEX.].


  (ii) The second step of resummation, that of squaring spin 
amplitudes, spin summation and phase space integration, is done numerically
in the Monte Carlo event generator. (There is no possibility of doing it 
analytically.) In the above CEEX scheme, the bulk of the IFI contribution is 
present in the resummed/exponentiated real+virtual form-factor 
and in the interferences emerging from squaring multiphoton spin amplitudes. 
Smaller contributions will remain hidden in the nonsoft finite remnants
[For instance, spin amplitudes of the γ-Z box are split into an
IR-divergent part, which is moved to an exponentiated form-factor, and 
the remaining IR-finite remnants are incorporated in the multiphoton spin 
amplitudes.]
The treatment of the nonfactorizable γγ and 
γ Z α^1
boxes in the above resummation scheme can be seen explicitly 
in Eqs. <ref>, <ref> for the matrix element in and in Eqs. (21-24) of Ref. <cit.> for the CEEX matrix element 
in .
 

A somewhat more detailed overview of the CEEX factorization and 
resummation in QED is given in Secs. C.2.7 and C.3 of the recent 
review <cit.>, while complete details can be found in 
Refs. <cit.>.
The above CEEX factorization and resummation of the universal QED corrections
will work equally well for extensions of the SM (BSM), 
provided the BSM predictions are formulated at the amplitude level.
So far,  for the e^-e^+→ ff̅ +nγ process
is the only implementation of the CEEX scheme.
In Ref. <cit.>, it is argued that 
in the context of the FCCee project, it is urgent to implement it
also for the Bhabha process.

The so-called deconvolution of QED effects from LEP experimental data which was
instrumental in the construction of 
pseudo-observables <cit.>
can also be reorganized using CEEX technique.
In the LEP version, it was done using ZFITTER, TOPAZ0 and 
KORALZ, and it was proven to be acceptable within the LEP precision 
goals <cit.>, 
but the validity of this procedure is not automatically granted for 
FCCee precision.  In Sec. C.3 in Ref. <cit.>, 
a proposal is made for extending it to higher precision 
by exploiting the CEEX factorization and resummation scheme in the MC 
implementation.
An updated discussion on the same theme is also included in Sec. 5.7
of a more recent paper, Ref. <cit.>.
Validation of a more powerful scheme
of removing QED effects from experimental data at the precision 
level of the FCCee experiments will require a lot of numerical studies
of the type done in Ref. <cit.>,
and most likely the development of the MC programs even more powerful
and versatile than . 

The plan of the paper is the following:
Section <ref> explains the origin and character of the IFI effect
in the angular distribution of the e^-e^+→μ^- μ^+ process.

Section <ref> describes a new partly analytic, partly numerical 
calculation in the semisoft approximation and its software implementation, 
.  “Semisoft” means that the upper limit on the total photon energy 
is smaller than the total center-of-mass energy s^1/2, 
but in the presence of narrow resonance Z
it can be smaller or bigger than Γ_Z/M_Z.

Multiphoton spin amplitudes (spin amplitudes) are defined 
(Sec. <ref>),
in such a way that they reproduce the CEEX matrix element in the semisoft 
regime.  Squaring and summing spins is also done analytically, 
and finally, a phase space integration over photon angles and energies
is also performed analytically keeping total real photon energy fixed
(Secs. <ref> and <ref>).
Multiple sums over photons and phase space integrations are done 
(exactly) in a straightforward way, with a minimal use of Mellin-Fourier 
transforms
[Mellin transforms are used merely
  as generating functionals for reorganizing combinatorics of the multiple
  sums over photons.].

The IFI effect appears in the resulting muon angular distribution.
The final analytic result involves a fourfold convolution over radiator 
functions of the initial state radiation (ISR), final state radiation (FSR) 
and two functions due to initial-final state interference (IFI) 
(Sec. <ref>).
The remaining integration over the phase space is delegated to a numerical 
Monte Carlo method.  

Validity of the formula is formally extended to the 
full phase space, such that the ISR and FSR radiator functions can be upgraded 
with the known nonsoft QED contributions up to α^2 
(Sec. <ref>).
Electroweak α^1 corrections are also included at the same
level as in , that is using the  library <cit.>.
The remaining 5-dimension integration over the fourfold convolution and muon
angle is slightly reorganized for numerical integration 
(Sec. <ref>)
using the universal  MC tool <cit.> 
for integration and simulation.
The resulting MC generator  is ready for use in
 Sec. <ref>.

With all the above distributions and technicalities in place, 
a wealth of numerical results produced by  and  is presented 
in Sec. <ref>.
In particular, in Sec. <ref> 
it will be checked in the calibration exercise
that for the matrix elements with resummation and without IFI, the
three programs , and  agree within ∼ 10^-5 precision
for σ_ tot and .
Note that in the case where IFI is switched off, this kind of comparison
of  with the numerical tool KKsem based on analytic
exponentiation was already done in Ref. <cit.>.
The new thing here is the inclusion of the IFI.

In Sec. <ref>, the IFI effect in  will be examined 
for three energies √(s)=10,87.9,94.3GeV 
as a function of the cutoff on total photon energy,
comparing results of  and .
In Sec. <ref>, we shall focus on the difference of 
between √(s_+)=94.3GeV and √(s_-)=87.9GeV, which is directly
related to the measurement of α_QED(M_Z).
Secs. <ref> and <ref> will be devoted
to estimating higher-order QED uncertainties by means of comparing
results for  for several variants of the QED matrix element in .

Section <ref> summarizes the results,
focusing on the uncertainties in the QED calculation of the IFI effect 
in .

Three appendices include details of the analytic phase space integration
in the semisoft approximation (Appendix <ref>),
the kinematical mapping used in  for the four-fold convolution 
(Appendix <ref>) and a recollection of some old 
α^1 analytic formulas without exponentiation
(Appendix <ref>) to be used in the numerical studies in 
Sec. <ref>.







§ PHYSICS OF IFI
 







Any efficient evaluation of IFI in perturbative QED must
be based on a good understanding of the basic physics governing this phenomenon.
Let us consider the process e^-e^+→μ^- μ^+ accompanied with any number
of real and virtual photons, illustrated schematically in
Fig. <ref>.
In the case of final fermions emitted at wide angles, IFI can be neglected, 
and in the case of total photon energy (ISR+FSR) limited to K< E=√(s)/2,
the angular distribution is uniformly lowered by a θ-independent
Sudakov form factor
[The subscript “virt” appears because virtual corrections
 feature -∫_0^E, while real emissions add +∫_0^K,
 so that the uncompensated remnant -∫_K^E is of virtual origin.
]

    dσ/dcosθ(K) ≃dσ_ Born/dcosθexp[
      -∫_K^E dk^0/k^0( 2α/πlns/m_e^2
                    +2α/πlns/m_μ^2)_ virt]
    =dσ_ Born/dcosθ 
       e^-Δ(K/E).

The above relation is illustrated schematically in Fig. <ref>.






Photon emission is, however, suppressed in the small-angle limit θ→ 0,
as illustrated schematically in Fig. <ref>, simply because 
the outgoing muon inherits most of the electromagnetic field accompanying the
incoming electron of the same charge as the muon; hence there is no need for
the compensating action of bremsstrahlung.
In fact, bremsstrahlung dies out completely at θ=0,
and it is the IFI contribution which kills both ISR and FSR.
The virtual form factor in the angular distribution 
at t → 0, s-|t|-|u|=0, |u|→ s becomes 

    Δ=
     ∫_K^E dk^0/k^0( 2α/πlns/m_e^2
             +2α/πlns/m_μ^2
             -4α/πln|t|/|u|)
    →∫_K^E dk^0/k^0( 2α/πln|t|/m_e^2
             +2α/πln|t|/m_μ^2)
    ≃ 0.






On the other hand, in backward scattering,
illustrated schematically in Fig. <ref>, the situation 
is completely different.
The electromagnetic field accompanying e^- has to be replaced by
that of μ^+, hence the violent compensating action of the bremsstrahlung
is much stronger than for wide angles.
Here we have u → 0 (c→ -1 side),s-|t|-|u|=0, |t|→ s. 
Thus, IFI enhances the total QED correction by a factor of 2:

    Δ=
     ∫_K^E dk^0/k^0( 2α/πlns/m_e^2
             +2α/πlns/m_μ^2
             -4α/πln|t|/|u|)
    →∫_K^E dk^0/k^0( 4α/πlns/m_e^2
             +4α/πlns/m_μ^2),

creating a dip in the muon angular distribution for backward scattering
(in the presence of a cutoff on the total photon energy, as previously).




In reality, the distribution of cosθ far from the resonance
appears as shown in Fig. <ref>
for a relatively strong cutoff on total photon energy (2% of the beam energy).

The presence of a narrow resonance significantly changes the pattern of 
QED cancellations.
Let us analyze briefly how the real and virtual corrections combine
at a resonance position √(s)=M_Z. 
                         

  * For pure ISR, the virtual correction is
∼ -2α/πlns/m_e^2lnE/λ,
as without a resonance, while the real contribution
is cut by the resonance profile
∼ +2α/πlns/m_e^2lnΓ_Z/λ.
The resulting cross section σ(K) is suppressed by
the remnant virtual factor
[1-2α/πlnM_Z/Γ_Z]_ virt 
for any cut above the resonance width, K>Γ_Z.

  * The effect of FSR is the same as in the case without a resonance, i.e.
 σ(K) is suppressed by the remnant virtual factor
[1-2α/πlns/m_μ^2lnE/K]_ virt.

  * The case of IFI is most complicated. The virtual correction
∼ -4α/πlnt/ulnΓ_Z/λ 
is cut by the resonance (contrary to the ISR case).
The real correction 
∼ +4α/πlnt/ulnΓ_Z/λ 
is also cut by the resonance (similar to the ISR case).
The resulting dσ(K)/dΩ is strongly power-suppressed 
by the Γ_Z/M_Z factor for any cut above the resonance width, 
K>Γ_Z!  For an energy cut below the resonance width,  K<Γ_Z,
IFI starts to rise logarithmically,
i.e. the suppression factor is
∼  1-2α/πlnt/ulnΓ_Z/K.

Away from the resonance, IFI gradually changes
to the previous nonresonant case,
and in the entire neighborhood of the resonance,
a QED calculation including photon resummation
at the amplitude level (CEEX) is mandatory.




The above mechanism is clearly illustrated in Fig. <ref>,
where the IFI contribution to  is shown as a function of
v_max, the cutoff on the total photon energy in units of the beam energy.
[More precisely, v=1-M^2_μμ/s. 
Here, we have temporarily used the ALEPH definition of v.]
As we see, far from the Z resonance and for a loose
photon energy cutoff, ≃ 2% and grows for stronger cutoffs.
In the middle of the resonance, it is strongly suppressed, < 0.1%,
and starts to grow below v_max≃Γ_Z/M_Z ≃ 0.02.
Remarkably, at the other two energies √(s)≃ M_Z ± 3.5 GeV,
Γ_Z/M_Z suppression is still quite strong, more than factor 1/5.

On the methodology side, although we are interested mainly in
the IFI effect off the Z peak, at √(s)≃ M_Z± 3.5 GeV,
it is worth also keeping an eye on √(s)=M_Z and energies far away 
from the resonance.  This is simply because any technical problem or 
mistreatment of physics which may cause a small effect at 
√(s)≃ M_Z± 3.5 GeV could be magnified there, making it
easier to trace it back and eliminate.
This is why we shall often compare our principal results 
with the results at  √(s)=M_Z and  √(s)=10 GeV.






§ NEW CALCULATION IN THE SEMISOFT APPROXIMATION AND  MONTE CARLO



As outlined in the Introduction, the aim of this section is to
describe in a detail all ingredients in the newly developed 
 Monte Carlo which will be served in the next sections for 
validation of the  predictions for (s) at a precision at least a 
factor of 10 better than in the past.
 is not a true event generator because photon momenta 
are partly integrated out analytically. 
Nevertheless, it provides weighted MC events with explicit muon four-momenta;
hence angular distributions of muons with an arbitrary cutoff on the total
photon energy can be calculated.

In Sec. <ref>, the multiphoton 
matrix element (spin amplitudes) will be defined.
In Secs. <ref> and <ref>
the above matrix element will be squared, and spin summation and
and phase space integration will be done partly analytically.
The IFI effect appears in the resulting muon angular distribution.
The resulting formula involves a fourfold convolution over radiator functions
of the initial state radiation (ISR), final state radiation (FSR) and
two functions due to initial-final state interference (IFI).
(Further analytic integration is not possible.)

In addition, in Sec. <ref>, the phase space integration 
is extended to the full phase space, and matching with the 
known α^1 and α^1 results for the ISR and FSR 
radiator functions is performed.
The radiator functions are convoluted with the 
effective Born spin amplitudes in which EW corrections are included.

In Sec. <ref>, it is explained how the remaining 
5-dimensional integration over the four radiator functions (ISR, FSR and 
2×IFI) and the azimuthal angle θ of the muon is performed 
numerically using a Monte Carlo method.
It is not trivial due to the presence of new types of singularities 
in the IFI radiator functions, different from 
the standard ones of the ISR and FSR radiator functions.
The newly developed computer program  is new software tool,
methodologically completely independent of ,
although it exploits some building blocks of ,
for instance the  library of the EW corrections <cit.>.






 §.§ Matrix element of multi-soft-photon emission 
 in the semisoft approximation



Let us consider the matrix element of the process

    e^-(p_1) +e^+(p_2)→μ^-(q_1) +μ^+(q_2) +γ(k_1)+…+γ(k_n)

near the Z resonance in the soft photon limit.
The standard kinematic variables
s=(p_1+p_2)^2,   t=(q_1-p_1)^2,   u=(q_2-p_1)^2
will be used.
Around any narrow resonance, the notion of the soft photon limit 
has to be refined.  In the framework of the
standard Yennie-Frautschi-Suura (YFS) <cit.>
soft photon resummation, one starts with all photons being very soft, 
i.e. k_i^0 ≪Γ_Z ≪√(s)/2.
Near the resonance, however, it is worth considering a wider soft photon range,
with k_i^0 ≪√(s)/2, but allowing photon energies comparable to
or even greater than the resonance width Γ_Z.
In the following, we shall refer to this regime as the 
semisoft approximation.
Following the notation of Ref. <cit.>, in the semisoft regime,
the matrix element of our process reads as follows:

    ^μ_1,μ_2,...,μ_n(p_i,q_j,k_l)=
     ∑_V=γ,Z∑_ 
     e^α B_4^V(s_I,t,m_γ)∏_i∈ I j^μ_i_I(k_i) ∏_r∈ F j^μ_r_F(k_r) _V(s_I,t ),
    
       
     s_I= P_I^2,   P_I=p_1+p_2-∑_i∈ I k_j,  
       
     j_I^μ(k)=eQ_e ( p_1^μ/kp_1-p_2^μ/kp_2),   
     j_F^μ(k)=eQ_f ( q_1^μ/kq_1-q_2^μ/kq_2),
    
       α B_4^V(s,t,m_γ)= 
    α B_4(s,t,m_γ)+ αΔ B_4^V(s,t,M_V^2).

The above formula involves a sum over the set of 2^n partitions 
{} = { I,F}^n,

    {} ={
     (I,I,I,...,I),
     (F,I,I,...,I),
     (I,F,I,...,I),
     (F,F,I,...,I), ...,
     (F,F,I,...,F)},

of photons among the initial and final state.
The meaning of the shorthand notation
i∈ I is that ∏_i∈ I includes all photons
with _i=I and similarly _r=F for r∈ F.

The form factor  B_4(p_i,q_i,m_γ) is the standard one appearing
in YFS resummation <cit.> for four charged particles
in the scattering process.
As stressed in Refs. <cit.>,
in the semisoft regime, an additional term in the form factor

    αΔ B_4^Z(s,t,M^2)=
    -2Q_e Q_f α/πln( t/u)
      ln( M_Z^2-s/M_Z^2 ),  M^2= M_Z^2-iM_ZΓ_Z,  Δ B_4^γ≡ 0,

must be included, but only in the resonant component of the amplitude.
For γ exchange, only the standard α B_4 of the YFS scheme
is needed, and α B_4^V is not present.
Most important is that, in the semisoft approximation,
the energy argument of the resonance propagator in the
Born matrix element _V must be shifted
by the total energy lost to initial state photons j∈ I,
[In the strict YFS soft limit this energy shift may be neglected.
  In the semisoft regime it could also be neglected
  for the γ-exchange part.
  For the sake of a better treatment of the collinear (mass) singularities,
  it is best to keep it everywhere.]
because of its strong energy dependence.
The same additional dependence on s_I 
also enters into the form factor α B_4^V.
The summation over all partitions of n photons between 
the initial and final state {I,F} is mandatory
in order to obey Bose-Einstein symmetry and gauge invariance.
Fermion spinor indices are implicit in _V.
The standard YFS virtual form factor B_4
is usually regularized with a photon mass m_γ. The mass regulator can
be removed once the real and virtual calculations are combined.

In the framework of coherent exclusive
exponentiation (CEEX) <cit.>, the above 
matrix element 
represents a zeroth-order CEEX matrix element defined throughout 
the entire phase space, including hard photons. 
Higher orders are also defined in the CEEX scheme, and implemented
for a finite number of the hard photons; see Ref. <cit.>.

The same matrix element can be rewritten in a compact form using a generating
functional formulation (Mellin-Fourier transform):

    ^μ_1,μ_2,...,μ_n(p_i,q_j,k_1,...,k_n)=
    
          
    =∑_V=γ,Z∫ d^4 Q d^4 x/(2π)^4
    e^ix · (P-Q) 
    e^α B_4^V(Q^2,t,m_γ)[ ∏_i=1^n J^μ_i(x,k_i) ]
     _V(Q^2,t )
    
       
    J^μ(x,k) =  e^-ik · x  j^μ_I(k) + j^μ_F(k).

The corresponding total cross section reads:

    σ(s)=
       1/ flux(s)∑_n=0^∞1/n!∫d^3 q_1/ q_1^0d^3 q_2/ q_2^0 ∏_i=1^n ∫d^3 k_i/ k_i^0 δ(P-q_1-q_2-∑_i=1^n k_i ) 
       ×^μ_1,μ_2,...,μ_n(p,q,k_1,...,k_n)
    [ _μ_1,μ_2,...,μ_n(p,q,k_1,...,k_n)]^*,

where P=p_1+p_2.
Note that in the above expression, the standard Lorentz invariant phase space
integral extends over the entire phase space,
as it would in the Monte Carlo implementation,
i.e., energy conservation naturally limits photon energies from above.
A cutoff on the total photon energy will be imposed later 
in our analytic calculations, in order to perform phase
space integration analytically.

In the semisoft approximation, the matrix element 
of Eq. (<ref>) is simple enough
that in the absence of experimental cuts
one can perform an analytic integration over photon angles and energies
and sum explicitly over photon multiplicities.
This is what we refer to as an analytic exponentiation.

The first important step will be reorganization
of the multiphoton distributions
keeping phase space integration untouched.
Later on phase space
integrations will be done in the semisoft regime
and finally contributions from hard photon phase space
will be reintroduced in the standard matching procedure at the level
of partly integrated distributions.





 §.§ Reorganization of multiphoton distributions




One may perform analytic reorganization of multiphoton distributions,
a necessary first step in the analytic exponentiation,
by means of a combinatorial reorganization of
the sum over photons without using the generating functional
(Mellin-Fourier transform) formulation of Eq. (<ref>).
This method was developed in Ref. <cit.>,
albeit for the resonant component only.
Another alternative method of integration/summation over semisoft photons
would employ a coherent states technique.  
This method was used, for instance, in Refs. <cit.>.
Let us start from the generating functional form of Eq. (<ref>),
which was used in the original YFS paper <cit.>,
although it was applied there for the simpler nonresonant case.
Of course, all three methods lead to an identical final result.

In the first step, let us introduce the usual eikonal factors
for photon emission from the initial state, final state, and
their interference:

    S_I(k) = -j_I(k) · j_I(k),  
     S_F(k) = -j_F(k) · j_F(k),  
     S_X(k) = -j_I(k) · j_F(k),

and Fourier representations of the δ-functions
of the phase space:[
The overall 4-momentum conservation δ^(4)-function will be present implicitly
in the next steps.]

    σ(s)=
    
       
    =1/ flux(s)∑_n=0^∞1/n!∑_V,V'∫d^3 q_1/ q_1^0d^3 q_2/ q_2^0  d^4 Q d^4 x/(2π)^4 d^4 Q' d^4 x'/(2π)^4
    e^ix · (P-Q) -ix' · (P-Q') d^4 y/(2π)^4
    e^ iy · (P-q_1-q_2) 
        ×∏_i=1^n ∫d^3 k_i/ k_i^0 [ e^-ik_i · (y+x-x') S_I(k_i) 
         +e^-ik_i · (y+x)    S_X(k_i)
         +e^-ik_i · (y-x')   S_X(k_i)
         +e^-ik_i · y        S_F(k_i)
     ]
    
        ×_V(Q,t) ^*_V'(Q',t) 
    e^α B_4^V(Q^2,t,m_γ) +α B_4^*V'(Q'^2,t,m_γ)


In the above functional representation, the
summation over photon multiplicities (exponentiation) is trivial:
[Both the virtual functions B_4^V
  and real emission integrals over S-factors are regularized
  temporarily using a small photon mass m_γ, which will cancel
  in the final result.]

    σ(s)=
    1/ flux(s)∑_V,V'∫d^3 q_1/ q_1^0d^3 q_2/ q_2^0  d^4 Q d^4 x/(2π)^4 d^4 Q' d^4 x'/(2π)^4
    e^ix · (P-Q) -ix' · (P-Q') d^4 y/(2π)^4
    e^ iy · (P-q_1-q_2) 
        ×exp{∫d^3 k/ k^0 [ e^-ik · (y+x-x') S_I(k) 
         +e^-ik · (y+x)    S_X(k)
         +e^-ik · (y-x')   S_X(k)
         +e^-ik · y       S_F(k)
     ]
    }
        ×exp{α B_4^V(Q^2,s,t) +α( B_4^V'(Q'^2,s,t))^*
    } _V(Q,t) ^*_V'(Q',t)

The integrations over x,x' and y can be reorganized in order
to achieve a clear factorization into ISR, FSR, and IFI parts,
as shown in Appendix <ref>.

A slightly reorganized form of Eq. (<ref>)
with U^μ representing the
total photon momentum of pure FSR emission,
K^μ representing the total momentum of pure ISR emission, 
and with R^μ and R'^μ aggregating IFI photons
present in _V and _V'^* correspondingly, 
reads as follows

    σ(s)=
    1/ flux(s)∫d^3 q_1/ 2q_1^0d^3 q_2/ 2q_2^0 
          d^4 K  d^4 R  d^4 R'  d^4 U
    δ^4(P-q_1-q_2-K-R-R'-U)
    
       ×∫d^4z / (2π)^4 
    e^ iz· K +∫d^3 k/ k^0  e^-ik · z  S_I(k) ∫ d^4u/ (2π)^4 
    e^ iu· R +∫d^3 k/ k^0  e^-ik · u  S_X(k) 
       ×∫d^4u'/ (2π)^4 
    e^ iu'· R' +∫d^3 k/ k^0  e^-ik · u' S_X(k) ∫d^4y/ (2π)^4 
    e^  iy · U  +∫d^3 k/ k^0  e^-ik · y  S_F(k) 
       ×∑_V,V'=γ,Z_V(P-K-R) ^*_V'(P-K-R')
    
       ×exp{
      2α B_4(s,t,m_γ) 
      +αΔ  B_4^V((P-K-R)^2) +(αΔ  B_4^V'((P-K-R')^2) )^*
    };

see also the illustration in Fig. <ref>.






The role of the Mellin transform in the above
algebra was merely to provide compact bookkeeping of the complicated
sums in the multiphoton distributions,
without any modification of the underlying phase space integration.
At any step, we could go back to standard phase space without any cutoffs;
for instance Eq. (<ref>) can be rewritten as follows:

    σ(s)=
    1/ flux(s)∫d^3 q_1/ 2q_1^0d^3 q_2/ 2q_2^0 
          d^4 K  d^4 R  d^4 R'  d^4 U
    δ^4(P-q_1-q_2-K-R-R'-U)
    
       ×∑_n_1=0^∞1/n_1!∏_i_1=1^n_1d^3 k_i_1/ k^0_i_1 S_I(k_i_1)
    δ_K=∑_i_1=1^n_1 k_i_1 
    ∑_n_2=0^∞1/n_2!∏_i_2=1^n_2d^3 k_i_2/ k^0_i_2 S_X(k_i_2)
    δ_R=∑_i_2=1^n_1 k_i_2 
       ×∑_n_3=0^∞1/n_3!∏_i_3=1^n_3d^3 k_i_3/ k^0_i_3 S_X(k_i_3)
    δ_R'=∑_i_3=1^n_3 k_i_3 
    ∑_n_4=0^∞1/n_4!∏_i_4=1^n_4d^3 k_i_4/ k^0_i_4 S_F(k_i_4)
    δ_U=∑_i_4=1^n_4 k_i_4 
       ×∑_V,V'=γ,Z_V(P-K-R) ^*_V'(P-K-R')
    
       ×exp{
      2α B_4(s,t,m_γ) 
      +αΔ  B_4^V((P-K-R)^2) +(αΔ  B_4^V'((P-K-R')^2) )^*
    }.

This is a generalization of Eq. (88) in Ref. <cit.>, which
was obtained there using pure combinatorics,
without any use of the Mellin-Fourier transform.
Both virtual and real photon integrals are IR-regularized using finite
photon mass m_γ.

Another advantage of the compact Eq. (<ref>)
is that by means of adding and subtracting

    ∫_k^0≤ Ed^3 k/ k^0 [ S_I(k) +2S_X(k) +S_F(k) ],
       E=√(s)/2,

in the form-factor exponent, we obtain a manifestly IR-finite expression <cit.>:

    σ(s)=
    1/ flux(s)∫d^3 q_1/ 2q_1^0d^3 q_2/ 2q_2^0 
          d^4 K  d^4 R  d^4 R'  d^4 U
    δ^4(P-q_1-q_2-K-R-R'-U)
    
       ×∫d^4z / (2π)^4 
    e^ iz· K +∫d^3 k/ k^0  [e^-ik · z -θ(k^0<E)]  S_I(k) ∫ d^4u/ (2π)^4 
    e^ iu· R +∫d^3 k/ k^0  [e^-ik · u-θ(k^0<E) ]  S_X(k) 
       ×∫d^4u'/ (2π)^4 
    e^ iu'· R' +∫d^3 k/ k^0  [e^-ik · u' -θ(k^0<E)] S_X(k) ∫d^4y/ (2π)^4 
    e^  iy · U  +∫d^3 k/ k^0 [ e^-ik · y -θ(k^0<E) ] S_F(k) 
       ×∑_V,V'=γ,Z _V(P-K-R) ^*_V'(P-K-R'),
    
       ×exp{
      Y(p_i,q_i)
      +αΔ  B_4^V((P-K-R)^2) +(αΔ  B_4^V'((P-K-R')^2) )^*
    },

where (e^-ik · y -1)/k^0 is IR-finite for k^0→ 0
and the classic YFS form factor

    Y(p_i,q_i)=  2α B_4(s,t,m_γ)
          + ∫_k^0≤ Ed^3 k/ k^0 [ S_I(k) +2S_X(k) +S_F(k) ]

is also finite in the m_γ→ 0 limit.

Reintroducing an IR cutoff  on the real photon energies,
Eq. <ref> can be also rewritten in
the following equivalent form with the standard phase space integration
and without any Mellin transforms <cit.>:

    σ(s)=
    1/ flux(s)∫d^3 q_1/ 2q_1^0d^3 q_2/ 2q_2^0 
          d^4 K  d^4 R  d^4 R'  d^4 U
    δ^4(P-q_1-q_2-K-R-R'-U)
    
       ×∑_n_1=0^∞1/n_1!∫_ k^0_i_1>∏_i_1=1^n_1d^3 k_i_1/ k^0_i_1 S_I(k_i_1)
    δ_K=∑_i_1=1^n_1 k_i_1 
    ∑_n_2=0^∞1/n_2!∫_ k^0_i_2>∏_i_2=1^n_2d^3 k_i_2/ k^0_i_2 S_X(k_i_2)
    δ_R=∑_i_2=1^n_1 k_i_2 
       ×∑_n_3=0^∞1/n_3!∫_ k^0_i_3>∏_i_3=1^n_3d^3 k_i_3/ k^0_i_3 S_X(k_i_3)
    δ_R'=∑_i_3=1^n_3 k_i_3 
    ∑_n_4=0^∞1/n_4!∫_ k^0_i_3>∏_i_4=1^n_4d^3 k_i_4/ k^0_i_4 S_F(k_i_4)
    δ_U=∑_i_4=1^n_4 k_i_4 
       ×exp( -∫_<k^0<Ed^3 k/ k^0  S(k) ) ∑_V,V'=γ,Z_V(P-K-R) ^*_V'(P-K-R')
    
       ×exp{
      Y(p_i,q_i)
      +αΔ  B_4^V((P-K-R)^2) +(αΔ  B_4^V'((P-K-R')^2) )^*
    }.

The phase space integration in the above formula cannot be performed 
analytically.
(It is done numerically without any approximation in the  program.)
In the following, this phase space integration will be done analytically
in the semisoft approximation.




 §.§ Analytic integration over photon momenta



In the next step, we shall integrate over the photon angles
in Eq. (<ref>),
introducing the cutoff 
E_max = v_max 1/2 s^1/2,
v_max≪ 1,
on the total photon energy in order to simplify the phase space integral,
and staying within the semisoft approximation for the
multiphoton distributions,
as in Eq. (<ref>) and Eq. (<ref>).

Let us show how it is done for the 
initial state part of this multiphoton integral.
In the semisoft photon limit, the integrand of Eq. (<ref>)
has no dependence on the spatial components of K outside of the 
e^iz· K factor.  Typically, the Born matrix element and the resonant
form factor have a dependence on K^0 through

    (P-K-R)^2= P^2-2P·(K+R)+(K+R)^2 ≃ s-2√(s) (K^0+R^0),

but no dependence on the spatial components K⃗. Thus, the integral 
over K⃗ yields a factor δ^3(z⃗) leads to
[ See also Eq. (<ref>) in the next subsection
  for versions without a Mellin transform.
]

    ∫d^4 K  d^4z / (2π)^4 
    e^ iz· K +∫d^3 k/ k^0    S_I(k) [e^-ik · z -θ_k^0<E]   
       
    =∫  dK^0 ∫dz^0 /2π 
    e^ iz^0 K^0 +∫d k^0/ k^0 γ_I [e^-ik^0 z^0 -θ_k^0<K^0]   
    = ∫d K^0/K^0   F(γ_I) γ_I  
      ( K^0/E)^γ_I,

where the integration over photon angles resulted in

    γ_I=γ_I(s)=
     ∫d^3 k/ k^0 
          S_I(k) δ( 2k^0/√(s)-1).

The subtle point is that the elimination of
∫ d^3K⃗δ^3(K⃗-∑_i=1^n k⃗_i)
implies that we keep K⃗=∑_i=1^n k⃗_i
everywhere in the entire integrand.
Note that in , the above “recoil effect” in the Born matrix element 
and phase space integral is taken into account correctly 
for hard photons as well.  The function 

    F(γ)≡exp(-γ C_E )/Γ(1+γ)

is well known from YFS work (Eq. (2.44) in Ref. <cit.>) and 
is due to the competition of soft real photons for the available fixed total energy.

Similarly, we are able to integrate over FSR and IFI photons:

    ∫ d^4U d^4y/ (2π)^4 
    e^  iy · U  +∫d^3 k/ k^0  S_F(k) [e^-ik · y -θ_k^0<E ] 
    =  ∫dU^0/U^0 γ_F 
     ( K^0/E)^γ_F F(γ_F),
    
       ∫ d^4 R  d^4u/ (2π)^4 
    e^ iu· R +∫d^3 k/ k^0    S_X(k) [e^-ik · u-θ_k^0<E ]  
    = ∫dR^0/R^0 γ_X  
       ( K^0/E)^γ_X F(γ_X),
    
       ∫ d^4 R'd^4u'/ (2π)^4 
    e^ iu'· R' +∫d^3 k/ k^0  S_X(k) [e^-ik · u'-θ_k^0<E ] 
    = ∫dR'^0/R'^0 γ_X 
      ( K^0/E)^γ_X F(γ_X),

where

    γ_F=γ_F(s)=
     ∫d^3 k/ k^0 
          S_F(k) δ( 2k^0/√(s)-1),
    
       γ_X=γ_X(cosθ)=
     ∫d^3 k/ k^0 
          S_X(k) δ( 2k^0/√(s)-1),

and θ is the angle the between the momenta p_1 of e^- and 
q_1 of μ^-.

Inserting all the above into Eq. (<ref>),
we finally obtain a
compact elegant formula:

    σ(s,v_max)=
    3σ_0(s)/8∑_V,V'∫_0^1 dv_I   dv_F  dr  dr' ∫d cosθ dϕ/2 θ(v_max-v_I-r-r'-v_F)
    
       ×ρ(γ_I,v_I) ρ(γ_X,r) ρ(γ_X,r') ρ(γ_F,v_F) 
     e^Y(p_i,q_i) 
       ×1/4∑_τ{e^αΔ B_4^V(s(1-v_I-r))^V_τ(v_I+r,c) 
              [e^αΔ B_4^V'(s(1-v_I-r'))^V'_τ(v_I+r',c)]^*
        },

where the Born spin amplitudes of Appendix <ref> are used and we define

    ρ(γ,v) = F(γ)   γ    v^γ-1,  
      v_I=2K^0/√(s),  
      r=ln2R^0/√(s),  
      r'=ln2R'^0/√(s),  
      v_F=2U^0/√(s).

The appearance of the real part [_V _V'^*] 
has resulted from symmetrization over r and r'.
The overall structure of the above integral is illustrated 
in Fig. <ref>.

Note that the YFS function  ρ(γ,v) obeys the following nice 
convolution rule
(related to the fact that it represents a Markovian process):

    ∫ dv_1 dv_2 δ(v-v_1-v_2) ρ(γ_1,v_1) ρ(γ_2,v_2) 
     = ρ(γ_1+γ_2,v),

but this feature cannot be exploited to simplify 
the integral of Eq. (<ref>),
because of the peculiar dependence of the matrix element on r and r'.
Let us stress that  the double convolution over ISR photons, 
separately for the Born amplitude and its conjugate  
seen in Eq. (<ref>),
is the landmark feature of the semisoft exponentiation
pioneered in Refs. <cit.> 
and implemented in .




 §.§ Matching of analytic exponentiation with fixed orders




Any matching of analytic exponentiation with fixed-order calculations 
must address the inclusion of the hard photon phase space beyond the 
semisoft regime represented in Eq. (<ref>)
with the cutoff v_max≪ 1 on the total photon energy.
The above matching will follow past examples 
in Refs. <cit.>.
It will result in the formula valid for  0<v_max≤ 1.

In order to match analytic exponentiation
with known analytic α^1,2,3 QED results for ISR and FSR
and compare with  over the entire phase space,
let us extrapolate the formula of Eq. (<ref>) 
beyond the semisoft regime to the entire range of the variable

    v= 1-M^2_μμ / s,    v∈ (0,1),

replacing soft photon approximation

    v = v_I+v_F+r+r',

with a multiplicative ansatz guided by the collinear kinematics,

    1-v = (1-v_I)(1-v_F)(1-r)(1-r').

With this ansatz, Eq. (<ref>) takes the form

    σ^(0)(s,v_max)=
    3σ_0(s)/8∑_V,V'∫ dv  dv_I   dv_F  dr  dr' δ_1-v=(1-v_I)(1-v_F)(1-r)(1-r') θ_v_max>v 
       ×∫d cosθ dϕ/2 ρ(γ_I(s),v_I) ρ(γ_F(s(1-v_I)(1-v_F)),v_F) ρ(γ_X(c),r) ρ(γ_X(c),r') 
     e^Y(p_i,q_i) 
       ×1/4∑_τ{e^αΔ B_4^V(  s(1-v_I)(1-r))^V_τ(1-(1-v_I)(1-r), c) 
                    
             [e^αΔ B_4^V'(s(1-v_I)(1-r'))^V'_τ(1-(1-v_I)(1-r'),c)]^*
        },

where c=cosθ and he have
inserted also the Born spin amplitudes of Appendix <ref>,
From now on we may use 0<v_max≤ 1.

In the numerical comparison of the above O(α^0)_ exp formula
with , it is worth including numerically significant 
O(α^2) contributions from the trivial phase integration.
It was shown in Ref. <cit.> (see Eq. (206) there)
that the following substitution does the job:

    ρ(γ_I,v_I) →ρ_I^(0)(γ_I,v_I)=
      ρ(γ_I,v_I) 
      exp[ 1/4γ_I
               +α/π( -1/2 +π^2/3)  ]
       [ 1 -1/4γ_Iln(1-v_I) ],
    
       ρ(γ_F,v_F) →ρ_F^(0)(γ_F,v_F)=
      ρ(γ_F,v_F) 
    
                          exp[ 1/4γ_I
               +α/π( -1/2 +π^2/3)
               -γ_F/2ln(1-v_F)
          ]
      [  1 -1/4γ_Fln(1-v_F) ],

where γ_F = γ_F(s(1-v_I)(1-v_F)).
[We could also use  γ_F = γ_F(s(1-v)), but
we have checked that it leads to the same numerical results.]

In order to compare with  O(α^2)  calculations
(including non-IR contributions of IFI up to O(α^1))
it is also quite easy to upgrade the ISR and FSR radiator functions
in Eqs. (<ref>) to O(α^2):

    ρ(γ_I,v_I) →ρ_I^(2)(γ_I,v_I)= 
      ρ(γ_I,v_I) 
      exp[ 1/4γ_I
               +α/π( -1/2 +π^2/3)  ]
       [ 1+γ_I/4 +γ_I^2/8
                           
       +v_I( -1 +v_I/2) 
       +γ_I ( -v_I/2 -1+3(1-v_I)^2/4ln(1-v_I)  )
      ],
    
    
       ρ(γ_F,v_F) →ρ_F^(2)(γ_F,v_F)=
      ρ(γ_F,v_F) 
      exp[ 1/4γ_I
               +α/π( -1/2 +π^2/3)
               -γ_F/2ln(1-v_F)
          ]
    
                [  1+γ_F/4 +γ_F^2/8
      +v_F( -1 +v_F/2) 
           + γ_F ( -v_F/2 
                            +v_F(2-v_F)/8ln(1-v_I)  )
      ],

see Tables I and II  in Ref. <cit.>.

The resulting ISR+FSR+IFI formula of Ref. (<ref>)
with all the above upgrades of ISR and FSR radiator functions
(with the resummation of ln(Γ_Z/M_Z)) 
is now ready for the MC implementation.

We are also going to implement the following
formula in which IFI is completely neglected:

    σ^(0)_ no IFI(s,v_max)=
    3σ_0(s)/8∫ dv  dv_I   dv_F δ_1-v=(1-v_I)(1-v_F) θ_v_max>v 
       ×∫d cosθ dϕ/2 ρ(γ_I(s),v_I) ρ(γ_F(s(1-v_I)(1-v_F)),v_F) 
     e^Y(p_i,q_i) 1/4∑_τ|_τ(v_I, c) |^2,


In the final push towards inclusion of as many known fixed
order results as possible into the analytic exponentiation formula,
we include the complete O(α^1) virtual IFI contributions.
This amounts to adding the non-IR parts of the γγ and γ Z
box diagrams explicitly provided in Eqs. <ref> in
Appendix <ref>
to the Born spin amplitudes:

    _τ(s,t) →_τ (s,t)
       +_τ^{γγ}(s,t,m_γ)
       +_τ^{γ Z} (s,t,m_γ)
    
                         
       -2α B_4(s,t,m_γ) _τ(s,t)
       -αΔ B_4^Z(s,t) ^Z_τ(s,t).

This is done in the framework of the standard YFS-inspired reorganization
of the IR singularities, the same way as in the CEEX matrix element of ,
without any danger of double counting.
The additional subtraction of αΔ B_4^V prevents double counting
with the resummation/exponentiation of this term in the semisoft regime.

In addition electroweak and QCD corrections are also included
in coupling constants of Born amplitudes, the same way as in .
Both  and  use the  library of 
O(α^1) EW corrections <cit.>
(including some of O(α^2))
and the method in which EW corrections are inserted into Born-like
parts of the spin amplitudes in  is essentially the same as 
in ZFITTER <cit.>.
It is described in Eqs. (21-25) of Ref. <cit.>.

Finally, with all the above changes due to matching with  O(α^1)
and  O(α^2) known fixed-order corrections,
we are now ready to implement the results of analytic exponentiation
of Eqs. (<ref>) and (<ref>)
with radiator functions of Eq. (<ref>,<ref>)
and box insertions, using the Monte Carlo method.

Coming back to
the extension of the phase space in Eqs. (<ref>) and (<ref>)
we see that it has now a well defined meaning: for IFI switched off these formulas
coincide with the well known QED convolution formulas for the total cross section
<cit.>
including hard photons.
However,
for the angular distributions Eq. (<ref>) is not able to reproduce exactly
the α^1,2 angular distribution beyond the soft limit.
On the other hand, from the analysis of Ref. <cit.> it is known 
that it reproduces numerically very well
α^1,2 MC results without IFI near the Z resonance
for  calculated practically for any choices on the muon angle,
within realistic cutoffs including hard photon emissions.
Unfortunately,
the known nonsoft α^1 IFI contributions to angular distributions,
see Ref. <cit.> and Appendix <ref>
cannot be reproduced exactly
[It would cost adding two extra phase space integration variables
  in  in order to complete α^1. ]
by Eqs. (<ref>,<ref>).
The main aim of Eqs. (<ref>,<ref>) 
with radiator functions of Eqs. (<ref>,<ref>)
is to test soft limit of IFI implementation
in  in the presence of α^2 ISR, FSR and α^1 EW corrections.

How different are
the above analytic resummations of the IFI effect in the semisoft approximation
from the known similar calculations in
the literature <cit.>?
Although the starting point in terms of multiphoton amplitudes is the same,
the semisoft resummation in
Ref.<cit.> exploits techniques of coherent states and Mellin 
transform
for dealing with multiple sums over photons and phase space integration,
while our approach is based on the straightforward combinatorics 
and direct phase space integration.
References <cit.> attempt to do final phase space
integration analytically, while in our approach we perform them numerically,
gaining more flexibility in the matching with finite order results
and in numerical comparisons with .

As a parenthetical remark, let us remark that
one may try to do some extra ad hoc simplifications,
strictly speaking not justified in the semisoft regime,
which may have some practical advantages in the parametrization of the MC results or data.
One example is the following variant of Eq. (<ref>)
which implements IFI in the approximate form:
[ There are more variants of this formula,
  for instance setting u=0 in Born matrix element and form factor, etc.]

    σ^(0)(s,v_max)=
    3σ_0(s)/8∫ dv  dv_I   dv_F  du δ_1-v=(1-v_I)(1-v_F)(1-u) θ_v_max>v 
       ×∫d cosθ dϕ/2 ρ(γ_I(s),v_I) ρ(γ_F(s(1-v_I)(1-v_F)),v_F) ρ(2γ_X(c),u) 
     e^Y(s,c) 
       ×1/4∑_τ|∑_V
         e^2αΔ B_4^V(  s(1-v_I)(1-u))^V_τ(c,1-(1-v_I)(1-u), c) |^2.

The simplification is due to neglecting r and r' dependence
in the Born matrix element and keeping the integration over u=r+r'.
The quality of the above approximation can only be judged
using numerical tests.





 §.§ Numerical integration methodology




Our aim is to perform numerically the 5-dimensional and 3-dimensional 
integrals in Eqs. (<ref>) and (<ref>)
using the Monte Carlo integrator FOAM <cit.>.
This is not quite trivial because the integrands of 
Eqs. (<ref>) and(<ref>)
are singular and nonpositive.
Singularities due to ρ_I of ISR and ρ_F of FSR
can be easily eliminated with 
the following simple mapping of variables
[FOAM can cope with these singularities even
  without such a mapping.]

    v=x_max  y_1^1/γ_I,   
      u=x_max  y_2^1/γ_I,   y_i∈(0,1).

The variable x_max=0.999... is a technical cutoff
introduced to avoid numerical instabilities near v=1.
The main problem is the integration over the two
more strongly singular and nonpositive ρ_X factors.
This occurs when

    γ_X(θ)= 
        2Q_eQ_fα/πln( 1-cosθ/1+cosθ)

becomes negative: γ_X(θ) = -β <0 in the forward
hemisphere, where cosθ>0.

In fact, one may think that in such a case the integral 
of Eq. (<ref>) does not make sense at all, because the
singularity r^γ_X-1=r^-β-1 from ρ_X is even not integrable!
However, a closer examination of the multiphoton integral which has led to 
ρ_X(r,γ_X(θ)) reveals that
the original distribution is in fact regularized with the familiar
plus-prescription (...)_+.

In order to understand the problem better, it is worth examining
the generic YFS multiphoton integral

    ∫_K^0<Ed^4 K  d^4z / (2π)^4 
    e^ iz· K +∫d^3 k/ k^0    S(k) [e^-ik · z -θ_k^0<E]   
       
    =∫_K^0<E  d^4 K 
    e^ -∫d^3 k/ k^0   S(k)  θ_<k^0<E∑_n=0^∞1/n!∏_i=1^n ∫_<k^0<Ed^3 k_i/ k_i^0   S(k_i) δ^4(K-∑_i=1^n k_i)
    
       
    =∫_0^E  d K^0 
    e^ -∫γd k^0/ k^0 θ_<k^0<K^0∑_n=0^∞1/n!∏_i=1^n ∫_<k^0<K^0γd k_i^0/ k_i^0  δ(K^0-∑_i=1^n k_i^0)
    
       
    =∫_0^E  dK^0 ∫dz /2π 
     e^ iz K^0 +∫d k^0/ k^0 γ [e^-ik^0  z -θ_k^0<K^0]   
       
    = ∫_0^Ed K^0/K^0 γ   F(γ)
      ( K^0/E)^γ
    = ∫_0^1 dv γ v^γ-1   F(γ)
    = ∫_0^1 dv ρ(γ,v)
    =  F(γ).
 
It is easy to check that the above integral
is always finite and well-defined for any choice of S=S_I,S_F,S_X, 
even for negative S and for negative γ!
Obviously, for γ>0, the singularity v^γ-1 is integrable and does
not require any regulation.  Closer inspection of Eq. (<ref>)
with an explicit IR regulator ϵ≪ 1 reveals that for any γ,
including γ=-β<0, the following holds:

    ρ(γ,v)=
    e^ -∫γd k^0/ k^0 θ_<k^0<K^0∑_n=0^∞1/n!∏_i=1^n ∫_<k^0<K^0γd k_i^0/ k_i^0  δ(K^0-∑_i=1^n k_i^0)
    
       
    = δ(v) F(γ) [1- ∫_^1 dv' γv'^γ-1]
     +θ(v-) F(γ) γ v^γ-1
    = F(γ) [ δ(v) + (γ v^γ-1)_+ ].
 
The standard plus prescription can be formulated
either in a regulator-independent way

    ∫_0^1 dv ϕ(v) (γ v^γ-1)_+ =
      ∫_0^1 dv  [ϕ(v)-ϕ(0)] γ v^γ-1,
 
or with an explicit regulator ≪ 1

    (γ v^γ-1)_+ 
      = γ v^γ-1θ(v-)
           - δ(v) ∫_^1 dv' γv'^γ-1
      = γ v^γ-1θ(v-) -δ(v)[1-^γ].
 
Of course, for γ>0, it becomes simpler,
because for → 0, we get
(γ v^γ-1)_+ →γ v^γ-1 -δ(v)
and ρ(γ,v) →  F(γ) γ v^γ-1.
However, the explicit IR regulator remains mandatory for γ<0.

As a closing crosscheck, let us verify that for γ<0 and regularized

    ρ(-β,v) = F(-β)  
             [δ(v) ^-β -θ(v>) β v^-1-β] ,

the basic convolution rule of Eq. (<ref>) still holds:
[In the → 0 limit, of course.]

    ∫ dv_1 dv_2 δ(v-v_1-v_2) ρ(-β,v_1) ρ(γ,v_2) 
     = ρ(γ-β,v).

In terms of the Markovian process, the function ρ(γ,v) for γ>0
represents adding more (soft) photons.
The other  ρ(-β,v) function is undoing that (backward evolution).

In the context of the explanation of the physics of IFI in 
Sec. <ref>,
the presence of ρ(γ,v) with negative γ in the
forward hemisphere in Eq. (<ref>)
is now perfectly understandable:
ρ(γ_X,r) ρ(γ_X,r') is undoing
part of the ISR and FSR photon emission coming from
ρ(γ_I,v)ρ(γ_F,u)!
See the upward (blue) arrow in Fig. <ref> 
for the corresponding graphical illustration.





In the numerical MC integration, it is not difficult to introduce a small 
IR regulator  into ρ(γ_X,r) when γ_X<0.
In the integrand for FOAM, this is done as a part of the mapping of
the integration variables r and r' into internal variables of FOAM.

Another issue is that the integrand becomes negative for γ_X<0, 
for r>, or for r'>.
This is handled in a standard way using
using weighted MC events with a nonpositive weight.
In the actual integration by means of FOAM, the modulus of the integrand
is used during the exploration stage,
while in the following MC calculation of the integral,
the MC events are weighted with the true signed distribution.
The distribution of the MC weights in the second stage has two peaks:
[This entails a certain loss of integration precision, but it
 turns out to be affordable.]
the bigger one close to +1 and smaller one near -1;
see Fig. <ref>.
More details on the mappings used in the construction
of the integrand for FOAM in  program are given 
in Appendix <ref>.











§ NUMERICAL RESULTS FROM  AND 




In this section, we present results from
the updated v4.22 of , also referred to as KKMCee,
the non-MC integrator .
[ uses Gauss quadrature programs to integrate
 analytic formulas up to 3 dimensions.
 It was developed at the time of preparing Ref. <cit.>
]
and the newly developed  MC integrator (simulator)
program , 
based on the C++ version of FOAM  <cit.>.
implements the 5-dimensional integral 
of Eq. (<ref>) including IFI,
together with its 3-dimensional variant without IFI
of Eq. (<ref>).
They will be often nicknamed in the following 
as 5 and 3 correspondingly.

Another subgenerator in  taking care of 2-dimensional
integration over v and cosθ will be used for
reproducing and/or implementing old pure O(α^1) 
results without resummation.

In 5 and 3, one may choose ISR and FSR structure functions
with soft photon exponentiation and QED corrections up to
O(α^0), O(α^1) and O(α^2),
as defined in Tables I and II in Ref. <cit.>.
Pure QED nonlogarithmic O(α^2) corrections are <10^-5,
hence are neglected for ISR, FSR and IFI.
They should be included and evaluated more precisely in the future.

The Born cross section in both 5 and 3
is implemented using two types
of subprograms of ,
both of them using spin amplitudes: either calculated in the scheme
[This is a variant of Kleiss-Stirling method of Ref. <cit.>.]
of Ref. <cit.> and labeled with GPS or CEEX,
or using spin amplitudes of KORALZ <cit.>
and labeled as EEX.
Note that it is not possible to use EEX Born for IFI implementation;
hence in 5 only GPS/CEEX Born amplitudes are implemented.

Electroweak and QCD corrections are included in  in the rescaled
coupling constants of Born amplitudes,
both for CEEX/GPS and EEX type, the same way as in .
Contributions from nonfactorizable γγ and γ Z boxes
are also included; see also Eq. (<ref>) for the details.
Both  and  use DIZET library of 
the O(α^1) EW corrections <cit.>
(including some of O(α^2))
and the method in which EW corrections are inserted into Born-like
parts of the spin amplitudes in  is essentially the same as 
in ZFITTER <cit.>.
It is described in Eqs. (21-25) of Ref. <cit.>.
This method protects completeness of the 
O(α^1) content of the EW corrections.
If there is any bias introduced in this method, 
then it has to be of O(α^2).

It should be kept in mind that in 5, hard photon corrections
are included in integrated form in the structure functions up to 
O(α^2) for ISR and FSR, while for IFI they are not 
included–only the finite parts of the virtual O(α^1) 
IFI corrections (γ-Z boxes) are included there.
(In 3, IFI is completely absent.)

The immediate short-term aim in this section is to prove 
that these programs correctly calculate
σ(v<v_max) and (v_max)
with physical and technical precision δ∼ 10^-4
and δσ/σ∼ 3· 10^-4.
This is a factor of 10 better than at LEP, but still a factor of 10
short of what needed for FCCee near the Z resonance.
An additional cutoff |cosθ|<c_max will sometimes be imposed.
An analysis for more realistic cuts will be presented
in a separate publication.  
The IFI effect in  depends strongly on the
cutoff on the total photon energy v_max,
which will typically be varied between v_max=0.002 and v_max=0.200.
As already pointed out in the Introduction, 
such a cutoff, stronger that in typical LEP data analysis, 
may be necessary at FCCee for the sake of better control of backgrounds 
and higher order QED effects.
Moreover, the expectation is that semisoft photon resummation employed in 5
(taking into account the energy shift due to ISR in the Z propagator)
will work fairly well in this cutoff range near the Z pole.
[On the other hand, strict YFS soft photon approximation
  neglecting the ISR energy shift in the Z propagator is expected to be
  adequate for our precision requirements 
  only for v_max≤ 10^-4.]

In the following analysis, 
event selection will be examined in terms of two variables only, 
cosθ for the angle between e^- and μ^-
and v=1-M^2_μμ/s.
The variable v represents approximately the total energy
of all ISR and FSR photons, in units of the beam energy.
(More results for realistic selection cuts will be shown in the next paper.)
Of course, once harder photons are allowed,
the definition of  cosθ is no longer unique.
For the  results, we will use the cosθ 
definition of Ref. <cit.>  
unless otherwise stated–see Sec. <ref>
for more discussion of other choices of cosθ and their precise definitions.







 §.§ Outline of the numerical investigations




We have conducted numerical studies with three different 
programs, ,  and  
featuring several variants of QED matrix elements
and different types of phase space integration.
For the convenience of the reader, we summarize in Table <ref>
all types of programs and QED matrix elements (M.E.) used in them.
The CEEX matrix element of  for IFI component
is rated in the table as α^1 because
of missing nonsoft α^2 parts of the pentabox diagrams
specified in Fig. 5 in Ref. <cit.>,
but the soft/infrared parts of these diagrams are included thanks
to the semisoft resummation technique.
It would be desirable to include these pentaboxes 
in a future version of 
[Another urgent desirable update of the M.E. in  would
be inclusion of the α^3ln^3(s/m_e^2) corrections.].


Let us outline the plan of the following tests
which will lead to new estimates of the theoretical uncertainty
of the IFI calculation:


  * In Sec. <ref>, we shall find 
that the influence of the choice of the muon scattering angle on
the measurement of  is negligible.

  * Section <ref> is devoted to a calibration exercise
in which the correctness of the MC integration is checked by
comparing the cutoff dependence of σ(v_max) and (v_max)
from three programs, ,  and  with IFI switched off.
It is done first for a maximally simple variant of the QED matrix element 
with resummation and then for the best one.

  * In the next step, in Sec. <ref>, the
IFI effect in  is examined in the results for  (v_max)
from  and  for a maximally simple and the best QED matrix element
separately for three energies s^1/2=10,87.9,94.3GeV.

  * Section <ref> is devoted to the difference
Δ A_ FB^ IFI(v_max)=A_ FB^ IFI(v_max,s_+)-A_ FB^ IFI(v_max,s_-)
in which IFI effect cancels. 
Results from  and  for this difference will be compared.
Δ A_ FB is directly related to the measurement of α_QED(M_Z) at FCCee.

  * Finally, in Sec. <ref> results for 
the energy difference Δ A_ FB^ IFI(v_max) from will be analyzed for QED matrix elements with an increasing level of 
sophistication in order to estimate its theoretical uncertainty due to 
missing higher orders of QED.





 §.§ On the choice of the scattering angle θ



In the limit when all photons are very soft, the
momenta of the final muons are back to back and the scattering angle θ
between e^- and μ^- is unique.
Once at least one photon becomes energetic, the final muons are 
not back to back and there are many possible definitions
of the effective θ.
Using θ^(1)=∠(e^-,μ^-)
or    θ^(2)=∠(e^+,μ^+)
is not a favorable choice experimentally,
because it does not exploit fully the power of the tracker detector,
which detects both μ^± equally well. 

An example of a choice favorable for experiments,
taking full advantage of the very good
angular resolution of the muon detectors (trackers), 
which is much higher than the energy resolution,
is that of ref <cit.>

    cosθ^∙ = y_1cosθ_1 - y_2 cosθ_2,
    
       
      y_1 = sinθ_2/sinθ_1+sinθ_2,  
      y_2 = sinθ_1/sinθ_1+sinθ_2,  
      y_1+y_2 = 1.


However, it was shown in Ref. <cit.> 
that analytic evaluation
of the IFI effect according to the O(α^1) QED matrix element
can be easily done using

    cosθ^⋆ = x_1 cosθ^(1) +x_2 cosθ^(2),  
    x_i= q_i^0/(q_1^0+q_2^0),   x_1+x_2=1.

We use this choice for
most of the numerical results presented in this work,
unless otherwise stated.
Moreover, in Ref. <cit.> compact analytic results 
were obtained for a charge asymmetry defined using the first moment 

    Ã^*_ FB= 3/2∫_-1^1 cosθ^* dσ/σ

instead of the conventional forward-backward asymmetry
=(σ_F-σ_B)/σ.

For , the choice of cosθ is irrelevant as long as
all photons are sufficiently
[Sufficiently from the point of view of the FCCee precision]
soft.
Once at least one photon becomes
energetic, the O(α^1) contribution calculated for
a well-defined choice of cosθ should be included in .
So far, this is not yet done – 
it should be done in the next version.
Most likely, the preferred choice for  will be
cosθ^⋆.

On the other hand,  is a regular MC event generator providing four-momenta of
both muons (and all photons), 
hence it provides a prediction for  with any definition of cosθ.
Let us examine, using ,
how different the QED predictions for 
are for the above two choices of θ when v<0.2. 
Figure <ref> shows that the difference
between A^∙_ FB and A^*_ FB is below expected FCCee experimental 
precision of δ∼ 3· 10^-5,
i.e. all of our analysis for cosθ^∙
is valid for cosθ^⋆ and vice versa.

Using , it is easy to examine the difference between
Ã^*_ FB and A^*_ FB.
Figure <ref> shows that such a difference
might be sizable, up to ∼ 1%.
However, the difference in the IFI component could cancel between two 
calculations – 
for instance, we have checked that it does cancel in the difference
between  and , for IFI switched on.












 §.§ Baseline calibration, ISR+FSR without IFI



Let us start with the baseline calibration of the MC tools at 
the precision level ∼ 10^-5
at √(s_-)= 87.9 GeV and √(s_+)= 94.3 GeV.
Although KKsem does not include IFI, it is still useful
for checking the normalization of both  and .
Of course, normalization is irrelevant for our main observable, 
(v_max), but it is still profitable to test it,
simply because some technical problem 
that would be evident in σ(v_max) could
produce a small annoying effect in  as well.
Thus, it is better to keep an eye on both of these.

As already underlined, our main aim in the present study 
is a precise prediction for  at two energies
√(s_±)= 87.9 near Z resonance.
However,  in order to get better confidence in the implementation
of the QED matrix element,
we will also check  at √(s)= 10 GeV,
where the Z resonance is negligible, and at √(s)=M_Z,
where the suppression of IFI due to the long life time of the Z is maximal.

Let us start with a purely technical test with IFI off 
at √(s)=94.3, 87.9 and 10 GeV,
presented in Fig. <ref>.
In the LHS plot, all cross sections σ(v_max)
are divided by the reference cross section from .
All calculations are at simplistic exponentiated O(α^0) 
QED including ISR and FSR, but without IFI.
Different types of the Born matrix element, EEX or GPS, are used.
Very good agreement is seen, up to statistical error
δσ/σ∼ 3· 10^-5.
The agreement for  is also very good,
essentially up to statistical error
δ∼ 1· 10^-5 at √(s_±).
The above equality of the  and  results
is very important, because it illustrates/proves the quality of the MC integrators
– it should be kept in mind that
for IFI off they integrate exactly the same 3-dimensional integrand.
Even more significant is the agreement,
to within δσ/σ∼ 1· 10^-5 near the Z resonance,
 of  with the other two
programs for the simplified EEX0 matrix element.
This is because
for MC statistics of 2· 10^10 events one may expect problems with
rounding errors in the accumulation of the weights in the 
histograms.[Running in parallel on 100 nodes and combining the
  histograms afterwards helps to reduce this problem.]
The slightly bigger discrepancy 
beyond statistical error of δ∼ 3· 10^-5
for 10 GeV is not yet statistically significant
and not so important for our aims.

We conclude that the technical precision of the MC numerical integration
in all three programs, ,  and , is satisfactory for our needs.
Moreover, the above test is also important due to the fact that the IFI effect
is added in  by reweighting MC events generated without IFI.
Hence, the technical precision established for the non-IFI mode persists 
when IFI is switched on.



















In Fig. <ref>, we continue baseline testing without IFI,
now with O(α^2) exponentiated ISR and FSR.
The relative differences δσ/σ between  and versus  are examined.
It is done for the CEEX/GPS and EEX Born matrix element.
The relative difference δσ/σ∼ 3· 10^-4 for confirms all older tests in Ref. <cit.>,
rated at the ∼ 1 · 10^-3 level.
[
  This is not a high priority, but we shall try later 
  to find the source of these differences in the normalization.]
On the other hand, the differences in  between and  or  are again of the order of the statistical error,
which is ∼ 3· 10^-5, except √(s)=10 GeV, where it is slightly bigger.

The main result of the tests presented in
Figs. <ref> and  <ref> is that the basic technical
precision (in the MC integration) of  and 3 
near the Z resonance is generally better than δ∼ 3· 10^-5.
The implementation of QED photonic corrections for ISR and FSR (no IFI)
up to O(α^2) 
was also tested at this precision level.
[ This is not true for , where subprograms
 with and without IFI are independent modules generating 
 their own different MC events.]




 §.§ IFI contribution to  from  and 



Let us now increase the level of sophistication by one important step –
including IFI.  
This will be done first in the simpler case
(A) for ISR, FSR and IFI at the level O_ exp(α^0) with exponentiation,
and next in the case 
(B) for exponentiated IFI at the level O_ exp(α^1),
accompanied by ISR and FSR up to exponentiated O_ exp(α^2).

In case (A), results for (v_max) from  and  are
shown in Fig. <ref>,
while in case (B), the results are shown in Fig. <ref>,
for energies √(s)= 87.9,  94.3,  10 GeV in both cases.
The absolute predictions for  from  and  are
seen in the LHS plots of the these figures.
The differences in  due to switching on the IFI contribution are
quite sizable and rising quickly for ≤ 0.06,
up to 5% for ≤ 0.002.

The IFI contribution to  is shown more clearly in the RHS plots of 
Figs. <ref> and <ref>, 
where the absolute predictions for the IFI effect in 
from  and  are presented.
The most important result is the one represented by the
red curve (c) in the RHS figures in Fig. <ref>.
It represents the difference between  and  for the IFI contribution.
This crucial difference is up to
δ∼ 5· 10^-4.
(It will be analyzed carefully one more time in the next section.)
It is definitely above the technical precision level
δ < 3· 10^-5,
established previously in case of IFI switched off.

















How can we understand the above result?
In the case of Fig. <ref> where both   and  
are at the same O_ exp(α^0) level for ISR, FSR and IFI,
with semisoft resummation of IFI, the source of the difference is a
different treatment of the matrix element 
far away from the infrared point v_max=0.
Remembering that
the energy shift in the Z-resonance propagator is properly taken into account
in the semisoft approximation,
the difference between  and  [curve (c)]
should be proportional to v_max and should vanish for v_max→0.
This is what we see in Fig. <ref>.
[The slight difference at v_max→ 0 for √(s)=10 GeV
 can be traced to small spikes in the 0.99<|cosθ|<1 range,
 to be examined separately. It goes away for realistic experimental cutoffs.]

In the case of Fig. <ref>, the difference between  and 
 should reflect the fact that in  the entire
O_ exp(α^1) real and virtual contributions are included,
while in 5, the O_ exp(α^1) real contribution is 
incomplete.  
This could increase the difference between  and .
In fact it changes sign and increases the difference by at most a factor of 2.  
This can be seen as unexpected.
In order to have an idea how big the O_ exp(α^1) 
real photon IFI contribution can be, 
we have also included this contribution
(curve (e)) in  Fig. <ref>,
subtracting the soft component in an ad-hoc manner.
As we see, curve (e) typically has the same sign as 
the difference between  and  shown in curve (c),
but is a factor of 3-4 bigger.
Apparently,  includes most of the O(α^1) 
hard photon IFI contribution.
[It would be interesting to include this missing 
  O_ exp(α^1) real photon IFI contribution in .]

The inclusion of QED  O_ exp(α^1) virtual corrections and
box diagrams was done in following the prescription of Eq. (<ref>).
The pure O(α^1) numerical results
in  Fig. <ref> were reproduced
using analytic formulas of Refs. <cit.>,
which are collected and tested numerically
one more time in Appendix <ref>.

In spite of the incompleteness 
of the O(α^1) IFI in , the above result
makes us confident that we are quite close to reaching our first
intermediate goal of controlling the IFI effect in 
at the level of δ∼ 10^-4
in the semisoft resummation regime (v_max≤ 0.06).

Let us also finally show just one example of the entire angular
distribution dσ/dcosθ from  and ,
simply because agreement in  does not necessarily imply 
agreement in the angular distributions.
In Fig. <ref>, such a comparison is done for a
relatively mild cutoff v_max=0.02 on the total photon energy.
The angular distributions agree to within 0.005% as expected.
 











 §.§ (s_±) from  and  in presence of IFI.



As explained in Ref. <cit.>, the QED coupling constant
α_ QED(M_Z)  is closely related to
(s_±), but the exact relation is not straightforward
and we are not trying to reproduce it.
We limit our interest to the propagation of errors from (s_±)
to α_ QED(M_Z), which is simpler and can be read from Eq. (4.9)
in Ref. <cit.>.
For our purpose, it will be enough to use a simplified version of this equation,

    δα_ QED/α_ QED|_M_Z≃δ(s_+) - δ(s_-)/(s_+) - (s_-) ,

which is valid for small δ(s_±) and/or when there are no
strong cancellations between them.
[We thank P. Janot for pointing this out to us.]
This will be true in the following numerical examples,
 and we shall show typically 
the numerator δ(s_+) - δ(s_-)
along with the uncertainties  δ(s_±).

Having the above in mind we reexamine the comparisons 
between  and  of the previous section for this Δ.






From now on, we impose a realistic cutoff |cosθ|<0.9 in the tests;
however the cutoff has little influence on the resulting .
To start with, in Fig. <ref>
we show (v_max) from  at √(s_±)
with IFI switched on/off and with the best QED matrix element in .
The  changes sign between these two energies.
On the other hand, IFI keeps the same sign, hence we expect partial
cancellation of the IFI effect in the α_ QED(M_Z).
We do not pursue the reconstruction of  α_ QED(M_Z)
and only plot the difference of IFI effect between two energies
in the LHS of Fig. <ref> as a guide.

We have produced the same figure for , but we do not show it
here, because it looks essentially the same as Fig. <ref>.
What is more interesting is to reexamine the difference between
 and 
    δ A_ FB(s_±) = A_ FB(s_±)|_-A_ FB(s_±)|_,

already shown in curve (c) of Fig. <ref>,
and its difference between two energies √(s_±)

    Δδ A_ FB = δ A_ FB(s_+) -  δ A_ FB(s_-)

relevant for the uncertainty in the measurement of α_ QED(M_Z).
We are interested in the above quantity primarily for IFI switched on.
This quantity is shown in Fig. <ref>, see curve (c) there.
It turns out that Δδ A_ FB≤ 2· 10^-4 
within the interesting range of photon energy cutoff v_max≤ 0.1.
In Fig. <ref>, we have also included two
dashed lines
[The dashed lines of the band are at  
 ± 1.1 · 10^-4  |A_ FB(s_+)-A_ FB(s_+)|_v_max→ 0 = ±0.57 · 10^-4.
]
marking the band of the present uncertainty 
δα_ QED /  α_ QED(M_Z) = 1.1 · 10^-4
according to Ref. <cit.>.
The aim of FCCee is of course to get substantially smaller error than that.

The main contribution to Δδ A_ FB in curve (c) comes
from the uncertainty
in the IFI implementation (most likely in ),
as can be seen from curve (d) in Fig. <ref>,
which represents Δδ A_ FB for IFI switched off.
The aim of future work will be to get
Δδ A_ FB≤ 3· 10^-5 for IFI on,
that is to the same level as for IFI off,
in the semisoft regime v_max≤ 0.06.








The above 
Δδ A_ FB≤ 2· 10^-4
can be treated as an (over)conservative estimate of the
uncertainty of the IFI prediction for in the semisoft regime, which is much better
than the LEP-era estimate but still not up to the needs of FCCee.
A less conservative estimate will be provided in the next section.




 §.§ On  for O_ exp(α^i), i=0,1,2 in 



The differences between  and  provide much valuable information,
because the two programs differ quite a lot technically 
(MC soft photon phase space integration versus analytic integration),
while implementing the same physics of QED corrections.
However,  alone offers interesting insight into missing higher order
QED corrections related to IFI.
[Provided we trust the smallness of the technical precision 
error of .]

In , one may choose three types of the QED multiphoton matrix element
with resummation at increasing sophistication levels, 
O_ exp(α^i), i=0,1,2.
In Fig. <ref> we examine differences
in the IFI contribution A_ FB^IFI(v_max) 
between O(α^2) and O(α^1)
and also between O(α^1) and O(α^0).
In all of them, IFI may be switched on or off.
Complete non-IR  O(α^1) corrections are included in the
O(α^i),  i=1,2 case while in the O(α^0)
case, only the IR part of exponentiated IFI is implemented.
In the most sophisticated case of the O_ exp(α^2) QED 
matrix element in ,
only pure nonlogarithmic photonic corrections are missing.
[In particular, the non-IR parts of QED pentaboxes are missing; 
  see Fig. 5 in Ref. <cit.>.]






In the LHS of Fig. <ref>, we show plots of the IFI component
for all three cases O_ exp(α^i),  i=0,1,2,
while in the RHS we see the differences, for the two energy 
points √(s_±).
[Differences in Fig. <ref> are obtained using MC weights
 event per event, so statistical errors are grossly overestimated.
 This explains the lack of fluctuations among bins.]

The most important difference in Fig. <ref>, 
between  for
O(α^2) and O(α^1), is below the
statistical error of 10^-4.
This can be treated as a measure of the missing QED photonic higher order
corrections in the  predictions for  for this particular type of
experimental cutoff, v_max<0.2 and |cosθ|<0.9,
near Z resonance, |M_Z-√(s)| ≤ 3.5 GeV.

Finally, let us remark that the MC results for  presented here
with a statistical precision of 10^-4 were obtained using ∼ 10^10
MC events generated in parallel runs on PC farms.
Reducing the statistical error to 10^-5 will be feasible, but not trivial.
However, higher precision may be also feasible with less MC events
using the technique of recording differences of the MC weights, 
as it was done in some plots shown in the following.





 §.§ More on the uncertainty of the ISR effect in .



In this section, we will present a few results from  
which, in particular, will
give us more insight on the ISR effects in  when 
IFI is switched on and off.
 
In Fig. <ref>, we show differences
between O_ exp(α^i), i=1,2 results
from  with the CEEX matrix element in the case of IFI switched off –
that is pure ISR and FSR effects.
In fact, the ISR effect is dominant here.
The variation is smaller than 3· 10^-5 and cancels
between the two energies √(s_±).
[Such a cancellation of the ISR effect 
  was already noticed in Ref. <cit.>.]






The same phenomenon is seen in Fig. <ref>,
albeit the differences are smaller, as expected.
Note also that in both of the above cases, the effect of ISR
is completely negligible for v_max≤ 0.05,
that is for cutoffs on photon energy interesting experimentally!





Finally, we switch on IFI and examine again the differences
between O_ exp(α^i), i=1,2 results
from  with the CEEX matrix element in the case of IFI switched on.
The results are shown in Fig. <ref>.
This is the most interesting result, because it shows
the indirect influence of ISR on the IFI contribution to .
Curve (c) shows that for the difference in  between the two 
energies √(s_±), the first and second order results agree 
to within ≤ 2· 10^-5.
The disagreement is larger than was seen in the previous graph  with IFI off
in the semisoft region v_max≤ 0.06.





One may conclude that the above result provides us strong indication that
the QED uncertainty in  from  is 
of the order of the expected FCCee experimental error δ≃ 3· 10^-5.
In the above plots, statistical MC errors are negligible,
because all differences between the various QED matrix elements
are calculated using differences between the weights for the same sets of
weighted MC events.











§ SUMMARY AND OUTLOOK



The extensive numerical results
presented in this work and summarized in Table <ref> allow us to
conclude that the technical and physical uncertainty of the prediction of for the IFI component of  near the Z resonance is of order ∼ 10^-4
(row 3 in Table <ref>).
This is definitely better than the state of art in the LEP era 
∼ 2· 10^-3 - 3· 10^-3.
For IFI switched off the technical precision ∼ 10^-5 was obtained
(rows 1-2 in Table <ref>).
Some of the results presented here indicate that the physical precision (higher orders)
of  the IFI component in  near the Z resonance from  is in fact
at the level ∼ 1.5 · 10^-5 
(row 7 in Table <ref>).
I.e., what is needed in the FCCee experiment proposed to
measure the QED coupling constant at the scale M_Z with this precision.
This would allow α_ QED(M_Z) to be determined to a precision
significantly better than the present estimate of Ref. <cit.>,
which is δα_ QED/α_ QED≃ 1.1· 10^-4.
However, more work is needed for achieving better confidence
in the technical precision and higher order
photonic QED corrections in the  results.
More work is also needed to estimate other missing nonphotonic QED corrections
(e.g. pair emission) and electroweak corrections.
Extension of the presented analysis to more realistic
experimental selections (cuts) is also desirable.

The newly developed auxiliary MC program  was instrumental in
the above achievement.  For more precise tests of ,
it would be profitable to include in the phase space of  
the exact contribution from nonsoft real O(α^1) emission
matched with semisoft analytic resummation.




Acknowledgment 

The authors would like to thank Patrick Janot, Maciej Skrzypek
and Zbigniew Wa̧s
for critical reading of the manuscript and useful comments.



APPENDIX



§ FACTORIZING THE EXPONENTIATED FORMULA



Starting from Eq. (<ref>), let us introduce 
∫ d^4z δ^4(z-y-x+x')=1,
∫ d^4u δ^4(u-y+x)  =1 and
∫ d^4u' δ^4(u-y+x') =1
in the Fourier expression, obtaining

    σ(s)=
    1/ flux(s)∑_V,V'∫d^3 q_1/ q_1^0d^3 q_2/ q_2^0  d^4 Q  d^4 x /(2π)^4 d^4 Q' d^4 x'/(2π)^4
    e^ix · (P-Q) -ix' · (P-Q') d^4 y/(2π)^4 e^ iy · (P-q_1-q_2) 
         ×∫ d^4 K d^4 z/ (2π)^4  d^4 R d^4 u/ (2π)^4  d^4 R' d^4 u'/ (2π)^4 
    e^ i(z-y-x+x')· K + i(u-y-x) · R +i(u'-y +x')· R' 
        ×exp{∫d^3 k/ k^0 [ e^-ik · (y+x-x') S_I(k) 
         +e^-ik · (y+x)    S_X(k)
         +e^-ik · (y-x')   S_X(k)
         +e^-ik · y       S_F(k)
     ]
    }
        ×exp{α B_4^V(Q^2,t,m_γ) +α( B_4^V'(Q'^2,t,m_γ))^*
    } _V(Q,t) ^*_V'(Q',t)

The lowest-order spin amplitudes _V,  V=γ,Z 
are, up to a normalization constant,
equal to the amplitudes _τ^V defined in 
Appendix <ref>,
but fermion helicities are temporarily suppressed.

Thanks to the above reorganization, we may clearly factorize the result into
contributions due to the ISR, FSR, and IFI components of 
multiphoton emission:

    σ(s)=
    1/ flux(s)∫d^3 q_1/ q_1^0d^3 q_2/ q_2^0  d^4 Q d^4 x/(2π)^4 d^4 Q' d^4 x'/(2π)^4
    e^ix · (P-Q) -ix' · (P-Q')∑_V,V'
       ×∫ d^4 K 
         d^4 R 
         d^4 R' 
    e^ i(-x+x')· K + i(-x) · R +i( +x')· R' 
       ×∫d^4z / (2π)^4 
    e^ iz· K +∫d^3 k/ k^0  e^-ik · z  S_I(k) ∫ d^4u/ (2π)^4 
    e^ iu· R +∫d^3 k/ k^0  e^-ik · u  S_X(k) 
       ×∫d^4u'/ (2π)^4 
    e^ iu'· R' +∫d^3 k/ k^0  e^-ik · u' S_X(k) ∫d^4y/ (2π)^4 
    e^  iy · (P-q_1-q_2 -K -R -R')  +∫d^3 k/ k^0  e^-ik · y  S_F(k) 
       ×exp{α B_4^V(Q^2,t,m_γ) +α( B_4^V'(Q'^2,t,m_γ))^*
    } _V(Q) ^*_V'(Q')
    
    
    
    
    
    
    
       =
    1/ flux(s)∑_V,V'∫d^3 q_1/ q_1^0d^3 q_2/ q_2^0 
    d^4 Q   d^4 Q' 
       ×∫ d^4 K 
         d^4 R 
         d^4 R' δ^4(P-K-R-Q) δ^4(P-K-R'-Q')
    
       ×∫d^4z / (2π)^4 
    e^ iz· K +∫d^3 k/ k^0  e^-ik · z  S_I(k) ∫ d^4u/ (2π)^4 
    e^ iu· R +∫d^3 k/ k^0  e^-ik · u  S_X(k) 
       ×∫d^4u'/ (2π)^4 
    e^ iu'· R' +∫d^3 k/ k^0  e^-ik · u' S_X(k) ∫d^4y/ (2π)^4 
    e^  iy · (P-q_1-q_2 -K -R -R')  +∫d^3 k/ k^0  e^-ik · y  S_F(k) 
       ×exp{α B_4^V(Q^2,t,m_γ) +α( B_4^V'(Q'^2,t,m_γ))^*
    } _V(Q) ^*_V'(Q').






§ MAPPINGS IN THE  INTEGRAND



The regularized radiator distribution for the IFI component
in the semisoft photon analytical exponentiation

    ρ(γ,v) = F(γ)( δ(v) ^γ 
        + θ(v-) γ v^γ-1),    ∫_0^1 dv ρ(γ,v) =F(γ) ≡ F_γ,

is valid for both positive and negative γ.
The regulator  should be smaller that any scale dependence
in the Born cross section times the target precision of the calculation.
In our case it should be below Γ_Z/M_Z by factor of at least
10^-4, i.e. < 10^-5 is recommended.
[In the actual MC runs, we use = 10^-6.

]

The distribution for  should be positive in the exploration phase; hence,

    ρ̃(γ,v)=|ρ(γ,v)|
    
    	= F_γ [ δ(v) ^γ 
                       + θ(v-)  |γ| v^γ-1]

is used.
The mapping from v to the internal variable r∈ (0,1) of is chosen such that its Jacobian compensates exactly ρ̃(v).
More precisely, v(r) is the solution of the equation

    r ∫_0^1 dv' ρ̃(γ,v') =
      ∫_0^v dv' ρ̃(γ,v') = F_γ  R(v).

Note that for γ>0 we have R(1)=1,
while for  γ<0 we get R(1)=2 e^γ-1 > 1.
Differentiating Eq. (<ref>) we get
F_γ R(1) dr = ρ̃(γ,v) dv ;
hence the Jacobian is

    J(v) = |dv/dr| = F_γ R(1) (ρ̃(γ,v))^-1.


For γ>0, the mapping 
[with R(1)=1 and R()=^γ]
is simply

    v(r)=0,   for   r<R()=^γ,
    
       
    v(r)= r^1/γ,   for   r>R().

The corresponding Jacobian is

    J(v) = 1/R()= ^-γ  for   v=0   and  
    J(v) = F_γ  (ρ̃(γ,v))^-1=v/rγ  for   v>.


For γ<0 the mapping 
[with R(1)=2^γ-1 and R()=^γ]
is more complicated:

    v(r)=0,   for   
     r<R()/R(1)=^γ/2 ^γ-1,
    
       
    v(r)= [ 2R() -r R(1) ]^1/γ
        = [ 2^γ -r (2^γ-1) ]^1/γ   for    r> R()/R(1).

The corresponding Jacobian reads

    J(v) = R(1)/R()= 2^γ-1/^γ  for   v=0   and  
    J(v) = F_γ  R(1)/ρ̃(γ,v)   for   v>.


In the second simulation stage,  generates weighted MC events
with the IFI component being w=J(v) ρ(γ,v).
In the case of γ>0, the weight (component) in  will be w=1 
for any v, while for γ<0 it will be w=R(1)
for v=0 and w=-R(1) for v>.

In addition, special care has to be taken in the case of γ→ 0,
that is for cosθ≃ 0, because in this region
the above mappings can be numerically unstable due to the limited
range of the exponent in floating-point arithmetic.
Because of that, when |γln|<Δ≪ 1,
many of the above formulas have to be expanded accordingly.
[
 The value Δ=10^-4 used now looks OK,
 as the error of ∼Δ^2=10^-8 is more than acceptable.]

For |γln|<Δ≪ 1 and γ >0,
the expanded distribution, mapping, and Jacobian read: 

    ρ̃(v)=ρ(v)
      = F_γ[
      δ(v) (1+γln) +θ(v>) γ/v],  
      R()= (1+γln),  R(1)=1,
    
       
    v(r)=0   for   r<R(),    
    v(r)= exp[-1/γ(1-r) ]   for   r>R(),
    
       
    J(v) = 1/R()  for   v=0   and  
    J(v) = F_γ  (ρ̃(γ,v))^-1  for   v>.

For γ<0, the expanded expressions 
with R()= 1+γln,  R(1)    = 1+2γln>1,
read:

    ρ̃(v)=|ρ(v)|
     =F_γ[
        δ(v) (1+γln) -θ(v>) γ/v],  
       
    v(r)=0,   for   
     r < R()/R(1)=1+γln/1+2γln,
    
       
    v(r)= exp[ 1/γ(1-r) R(1) ] ,   for    r> R()/R(1),
    
       
    J(v) = R(1)/R()=     
          for   v=0   and  
    J(v) =  F_γ R(1)/ρ̃(γ,v)  for   v>.








§ ZERO AND FIRST ORDER AMPLITUDES WITHOUT RESUMMATION



For constructing the semisoft photon 
analytical resummation and matching
with the fixed-order O(α^1) result,
we need the zeroth and first order amplitudes and distributions
in analytical form.
In particular, we will need the differential cross section
of the final muons, integrated over photon angles,
but keeping control over the photon energy.
The relevant results are scattered over several papers
<cit.>.
See also Refs. <cit.>,
where they are sometimes incomplete,
or given in a form not suitable for our purposes;
hence it is worth collecting them once more in this appendix.

Following the notation of Ref. <cit.>, the Born cross section
and charge asymmetry read as follows:

    d σ^(0)(s(1-v))/dc = 
      3σ_0(s) /81/4∑_,τ=±| _τ(v,c)  |^2 
     = 3σ_0(s)/8[ (1+c^2) (v)  +2c (v) ],
    
       _τ(v,c)
     = ^γ_τ(v,c) +^Z_τ(v,c)
     = (τ + c)  D_,τ(v),
    
       
      D_,τ(v) 
     = D^γ_,τ(v)+ D^Z_,τ(v)
     = qq̃/1-v +  g_g̃_τ/ζ -v,
    
       ζ = s-M_Z^2+i Γ_Z M_Z/s,  
     g_τ = g_V +τ g_A,  g̃_τ = g̃_V +τg̃_A,  σ_0= 4 απ^2/3s,

where c = cosθ, q=Q_e, q̃=Q_μ are electric charges,
,τ=± are twice the helicity of e^- and μ^-, and

    (v)=
    1/4∑_τ| D_,τ(v) |^2 =
      c_0/ (1-v)^2
    +  2c_1 /(1-v)(ζ-v) 
    +  c_2 / |ζ-v |^2 ,
    
       (v)=
    1/4∑_ττ| D_,τ(v) |^2 =
       2d_1 /(1-v)(ζ-v) 
    +  d_2 / |ζ-v |^2 ,
    
       
    c_0 = (qq̃)^2,   
    c_1 = qq̃ g_v g̃_v,   
    c_2 = (g_v^2+g_a^2)( g̃_v^2+g̃_a^2) ),  
       
    d_1 = qq̃ g_a g̃_a,   
    d_2 = 4 g_v g_a  g̃_v g̃_a.


The integration over cosθ results in

    σ^(0)=
      σ_0 1/4∑_τ | D_,τ(0)|^2,  σ^ (0)= ∫ 2cosθ^ dσ^(1)
    = σ_0 1/4∑_ττ | D_,τ(0)|^2
    
       
     A_ FB^(0) = 3/4⟨ 2cosθ^⟩^(0)
       = 3/4∫ 2cosθ^  dσ^(0)/σ^(0)
       = 3/4 ∑_ττ | D_,τ(0)|^2 /∑_τ  | D_,τ(0)|^2 
       = 3/4 (0)/(0).

Following the notation of Ref. <cit.>,
the noninterference O(α^1)
results with implicit integration over photon angles
and explicit integration over photon energy up to x=v_max read:

    A_ FB^(1)(x)
     = 3/4σ^(1)(x) /σ^(1)(x),    σ^ (1)(x)
      =  ∫_v<x 2cosθ^ d σ^(1),
    
       σ^(1)(x)/σ_0
    = [1 +F(x) ] (0) + W (x), 
     σ^(1)(x)/σ_0
     =[1 +F^(x) ] (0) + W^(x),
    
       
     W( x)= ∫_0^x dv [γ_I(s) P(v) 
       + q^2 α/πΔ_sδ(v) ]  (1-v) (v),
    
       
    W^(x)=
    ∫_0^x dv [γ_I(s) P(v) 
       +q^2 α/πΔ_sδ(v) 
       -q^2 α/π v ]  (1-v) (v),


    F(x) = ∫_0^x dv  [γ_F(s(1-v)) P(v) 
      +q̃^2α/πΔ_sδ(v) ],    Δ_s =   -1/2 +π^2/3 ,
    
       
     F^(x) =
     ∫_0^x dv [γ_F(s(1-v)) P(v) 
      +q̃^2α/πΔ_sδ(v) 
      -q̃^2 α/π v  ],
    
       
     P(v)= (1+(1-v)^2/v)_+=
      -δ(v) 3/4ln1/
      +θ(v-) 1+(1-v)^2/v,

In Ref. <cit.>, analytical integrations over v were done,
but for the purpose of the present resummation, we are more 
interested in the above unintegrated version.

In Ref. <cit.>,
the contribution of IFI was added to the above charge asymmetry, 
 but in a version that was integrated over v.
The unintegrated version
[The unintegrated version of U^_,τ (x)
   was obviously used in Ref. <cit.>, but was not explicitly 
   shown there. Also, U_,τ (x) was not provided there.]
including ISR+FSR+IFI with complete O(α^1)
for v∈ (0,1),
needed for resummation is as follows:

    A_ FB^(1) (x) 
    = 3/4σ^ (1)(x)/σ^(1)(x) 
    = 3/4∫_v< x cosθ^  dσ^(1)/σ^(1)(x) ,
    
       σ^(1)(x)   = σ^(1)(x)     +σ_0 U(x),  
    σ^(1)(x)= σ^ (1)(x) +σ_0 U^(x),
    
       
     U(x)= ∫_0^x dv ρ_X^(1)(v)  (1-v)(v,0)
      + 3qq̃α/π{
               A_γ ^γ(0)
              +A_Z ^Z(0) },
    
        
     U^ (x)= ∫_0^x dv ρ_X^ (1)(v)  (1-v)(v,0)
      +  2 qq̃α/π{
               A^_γ ^γ(0)
             + A^_Z ^Z(0) },
    
       ρ_X^(1)(v)= 2qq̃α/π{δ(v) [ 3ln1/δ]
    +θ(v-δ) (-3) 2-v/2v}
       ρ_X^ (1)(v)= 2 qq̃α/π{δ(v) [ 5 ln1/δ]
    +θ(v-δ) (-1)/(2-v) v[ 10(1-v) + 3 v^2 ]
    },



The combined contributions to the total cross section 
from real soft emission (interference part)
and virtual γγ and γ Z boxes 
can be deduced from the k_max→ 0 limit of
formulas in Ref. <cit.>:

    A_γ= -1/2 ,  
      A_Z= -ln|1-ζ| -ζ +(1-ζ)(2-ζ) ln-ζ/1-ζ.
    
       ^γ(0)= c_0 + c_1 /ζ^* ,  ^Z(0)    =  c_1 /ζ^* + c_2 /ζζ^* ,


The analogous contributions to σ^ (1) can be obtained
from formulas in Ref. <cit.>:

    A^_γ= 65/36 -i2/3π,
    
       
    A^_Z=
       31/9ζ -9ζ^2+4ζ^3
       -ln(1-ζ)( 15/2 -13ζ+12ζ^2-4ζ^3 )
    
            
       +ln(-ζ) ( 5 -17/3ζ +2ζ^2 )
       +4ζ (1-ζ)^3  (  Li_2 ( -ζ/1-ζ)
               -π^2/6).
    
       ^γ(0)=  d_1 /ζ^* ,  ^Z(0)=  d_1 /ζ^* + d_2 /ζζ^*.


The following combinations of the Born amplitudes are involved:

    (v,u)=
     1/4∑_τ( D_,τ(v)^* D_,τ(u) )=
    
          
    ={ c_0/ (1-v)(1-u)
    + c_1 /(1-v)(ζ^*-u) 
    + c_1 /(ζ-v)(1-u) 
    + c_2 /(ζ-v)(ζ^*-u) }
       (v,u)=
     1/4∑_ττ ( D_,τ(v)^* D_,τ(u) ) =
    
          
    ={ d_1 / (1-v)(ζ^*-u) 
    + d_1 / (ζ-v)(1-u) 
    + d_2  / (ζ-v)(ζ^*-u) },
    
       ^V(0)=1/4∑_τ   D^V_,τ(0)^* D_,τ(0),  ^V(0)=
     1/4∑_ττ  D^V_,τ(0)^* D_,τ(0), 
    V=γ,Z,

Let us remark that the following relations hold:

    (v)= (v,v),  (v)= (v,v),


    ^γ(0)+ ^Z(0) = (0),  ^γ(0)+ ^Z(0) = (0),












We also need the virtual box and real soft contributions
before integration over c=cosθ.
Spin amplitudes for two γγ box diagram and two γ Z 
box diagram contributions, 
normalized the same way as the Born spin amplitudes, read as follows:

    ^{γγ}_τ 
    = (qq̃)^2 (τ X_1(c) +X_2(c)),
    
       ^{γ Z}_τ 
    = qq̃  g_ q_τ  (τ Z_1(c) +Z_2(c)).

Their interference with Born amplitudes leads to the following contributions
[
We use
 (1/4)∑_ττ  qq̃  g_g̃_τ = d_1
  = qq̃  g_a g̃_a and
(1/4)∑_τ g_g̃_τ = c_1 = qq̃  g_v g̃_v
]


    dσ^γγ/dc
    =3σ_0/81/4∑_τ
     2[ ^{γγ}_τ ^*_τ(0,c) ]
    
       
    =3σ_0/81/4∑_τ
    2{
    (qq̃)^2 
        [ X_1(c)+cX_2(c) +τ (cX_1 +X_2(c)) ] 
        D^*_τ(0)
       }
       
    =3σ_0/8 qq̃ 
    2{ (c_0 +c_1/ζ^* )[X_1(c)+cX_2(c)] 
               +d_1/ζ^* [cX_1(c) +X_2(c)]  }
       
    =3σ_0/8
    qq̃ 
    2{
       (c_0 +c_1/ζ^* ) [F^γγ(c) -F^γγ(-c)]
            +d_1/ζ^*   [F^γγ(c) +F^γγ(-c)]
       },

where

    X_1(c)+cX_2(c)= F^γγ(c) -F^γγ(-c),
    
       
      cX_1(c) +X_2(c)=F^γγ(c) +F^γγ(-c),  
      c_± = 1 ± c/2,
    
       
    F^γγ(c)=
    2α/π{
     2( lnm_γ^2/s +iπ) c_+^2 lnc_-/c_+
      -1/2 c   ( ln^2 c_- +2iπ)
                 + c_+ ( ln c_- +iπ)
    }.


Similarly, for the γ Z box, we have:

    dσ^γ Z/dc
    =3σ_0/81/4∑_τ
     2[ ^{γ Z}_τ ^*_τ(0,c) ]
    
       
    =3σ_0/81/4∑_τ
    2{
      qq̃[ Z_1(c)+cZ_2(c) +τ (cZ_1 +Z_2(c)) ] 
      D^*_τ(0)
    }
       
    =3σ_0/8 qq̃ 
    2{( c_1 +c_2/ζ^*)[Z_1(c) +cZ_2(c)] 
              +(  d_1 +d_2/ζ^*)[cZ_1(c) +Z_2(c)] }
       
    =3σ_0/8 qq̃ 
    2{( c_1 +c_2/ζ^*)
                           [F^γ Z(c) -F^γ Z(-c)]
              +( d_1 +d_2/ζ^*)
                           [F^γ Z(c) +F^γ Z(-c)] }

where F^γ Z(c) is related  in a simple way to f(s,t,u)
of Ref. <cit.>:

    F^γ Z(c)  = 2α/π  c_+^2 s f(s,t,u),  
       F^γ Z(-c) = 2α/π  c_-^2 s f(s,u,t).






















In the  code, the γ Z box of Ref. <cit.>
is programmed as follows:

    F^γ Z(c)=
     lnt/ulnm_γ^2/(tu)^1/2
     -2lnt/uln^2-s/^2 
     + Li_2( ^2+u/^2)
     - Li_2( ^2+t/^2)
    
       
    +(M^2-s)(u-t-M^2)/u^2(
       ln-t/slnM^2-s /M^2
      + Li_2( ^2+t/^2)
      - Li_2( ^2-s/^2)
     )
    
       
    +(^2-s)^2/usln^2-s/^2 
    +(^2-s)/uln-t/^2 ,

where ^2= M_Z^2-M_ZΓ_Z,  t=-(1-c)s and u=(1+c)s.

Finally, the above box contributions have to be combined with (interference)
the corresponding soft real emission contribution

    dσ^ soft_X/dc=
    dσ^(0)/dc  q q̃α/π δ^ soft_X(c),
    
       δ^ soft_X(c) =
      4lnc_-/c_+ln s^1/2ϵ/m_γ
     +ln^2 c_- -ln^2 c_+
     +2  Li_2(c_+)-2  Li_2(c_-)

such that the usual cancellation of the IR regulator m_γ occurs,
leaving out the IR cutoff on photon energy v≤ϵ≪ 1.

Let us finally define explicit relations between integrated
and unintegrated virtual+soft contributions:

    3A_γ=
    ∫ dc ( F^γγ(c) - F^γγ(-c)
                    + δ^ soft_X(c) ),
    
       
    2A^_γ=
    ∫ 2c dc ( F^γγ(c) + F^γγ(-c)
                    + δ^ soft_X(c) ),
    
       
    3 A_Z=
    ∫ dc ( F^γ Z(c) - F^γγ(-c)
                    + δ^ soft_X(c) ),
    
       
    2 A^_Z=
    ∫ 2c dc ( F^γ Z(c) + F^γγ(-c)
                    + δ^ soft_X(c) ),


Finally, in Figs. <ref> and <ref>
we crosscheck the old analytical results with ,
in which the integration over cosθ (virtual)
and over photon energy v (real photon) is done numerically.
As we see, there is perfect agreement between old analytical
formulas and new results using 2.















h-physrev3





