# Heisenberg Picture for Open Quantum Systems

Nachiket Karve  Department of Physics, Indian Institute of Technology Kanpur,
Kanpur 208016, Uttar Pradesh, India.  R. Loganayagam  International Centre for
Theoretical Sciences (ICTS-TIFR), Bangalore 560012, India.

###### Abstract

In this note, we develop a framework to describe open quantum systems in the
Heisenberg picture, i.e., via time evolving operator algebras. We point out
the incompleteness of the previous proposals in this regard. We argue that a
complete Heisenberg picture for an open quantum system involves multiple
_image Heisenberg operators_ for each system observable. For a given system
observable, the number of such image operators is equal to the dimension of
the environment Hilbert space. We derive a perturbative expression, accurate
upto arbitrary orders in the system environment coupling, for these image
operators in terms of a single _one point operator_. This expression depends
non-linearly on the state of the environment. This perturbative expression can
equivalently be thought of as deforming the operator product on the Hilbert
space of the open quantum system. In the Markovian limit, the one point
operator evolves by an adjoint Lindblad equation. We illustrate these ideas
using a simple spin system.

##  I Introduction

Any quantum system which interacts with the environment is called an open
quantum system. The aim of this note is to understand how such systems can be
described in Heisenberg picture. The conventional treatments of open quantum
systems proceed via Schrodinger picture [3], i.e., one starts with the full
Hilbert space formed by the direct product of the Hilbert spaces of the system
and that of the environment and then defines the state of the system via a
reduced density matrix obtained by tracing out the environment degrees of
freedom. The dynamics is then encoded in the evolution of such a reduced
density matrix. The question we are asking is as follows : can one instead
describe the same dynamics via an evolution of observables/operator algebra ?

The answer to this question is non-trivial for the following reason : given a
system and an environment interacting with each other, consider a Heisenberg
operator which starts off initially as a purely system operator (i.e., it is a
tensor product of an operator acting on the system Hilbert space and the
identity operator on the environment Hilbert space). Under Heisenberg
evolution, such an operator generically evolves into an operator which is no
more a tensor-product– this is just the statement of entanglement stated in
Heisenberg picture. It is hence unclear a priori how to project this evolution
into an evolution of a single system operator, the ‘reduced Heisenberg
operator’ so to speak.

A minimal requirement of this projection is this : any $N$ point function of
the Heisenberg evolved system operators computed in the total system should
equivalently be computable through $N$ point function of such a ‘reduced
Heisenberg operator’. As we will argue below, it is impossible to satisfy this
requirement : _there is no single system operator whose $N$ point functions
will reproduce the $N$ point functions of Heisenberg evolved system
operators_. We will however see that we can reproduce the required $N$ point
functions if we allow for multiple time-dependent system operators (the number
of such operators is equal to the dimension of the environment Hilbert space).
We will call these the image operators acting on the system Hilbert space.

While the above statement follows from a straightforward application of
quantum mechanics, the implications are somewhat perplexing. Consider the
example of a quantum particle coupled to an environment with an infinite
dimensional Hilbert space. A direct application of the statement above implies
that while there is a single position operator for this particle in the
Schrodinger picture, one requires an infinite number of Heisenberg position
operators to fully describe all correlation functions. A natural question
would be to ask how these multiple image operators are related to each other ?
How do they differ in their dynamics ? One objective of this note is to answer
these queries.

We have not found a discussion of the aforementioned conceptual difficulties
and their resolution in the existing literature. There have been previous
attempts to understand open quantum systems in the Heisenberg picture [6, 11,
5, 15]. Some attempts include calculating the radiative corrections to atomic
transition frequencies [1], and studying the non-Markovian nature of bosonic
and fermionic baths [17]. Here the adjective ‘non-Markovian’ denotes the fact
that the bath/environment retains memory, viz., the environment correlation
functions do not decay more quickly than the system correlations.

When the correlations in the environment do decay quickly, one can write down
a “master equation” which determines the Schrodinger evolution of the reduced
density matrix [19, 14]. These master equations have a Lindblad form [13, 10,
12]. It is natural then to shift this evolution on to the system operators and
obtain a Heisenberg picture with operators evolving via the adjoint master
equation[4, 3]. This approach identifies a single system operator that
reproduces the correct one point function of a Heisenberg evolved system
operator and focuses on the evolution of this operator. It should be noted
however that this operator, which we will call the _one point operator_ , does
not reproduce the correct higher point functions.

As mentioned above, a full Heisenberg picture for an open quantum system
requires the full set of image Heisenberg operators to reproduce all
correlation functions correctly. The central aim of this note is to describe
how these image Heisenberg operators are related to the one point operator
that appears in the discussion above. We will see that a precise relation can
be written down to arbitrary orders in system bath coupling. This relation
however turns out to depend non-linearly on the state of the environment. As
we will demonstrate later, an alternate way to interpret this relation is to
think of it as a deformation of operator products within the Hilbert space of
open quantum system.

The original motivation of this work is to study open quantum field
theories,i.e., dissipative effective field theories obtained by integrating
out fast degrees of freedom, fluid dynamics being the paradigmatic example in
this regard. The Heisenberg picture is natural and convenient in this context.
More generally, solving for the Schrodinger evolution of the full reduced
density matrix might often be a difficult endeavour whereas focusing on the
Heisenberg evolution of a subset of observables of experimental interest could
be tractable. This point of view already underlies the hierarchy of equations
approach[2], where the evolution of $N$-point correlation functions is
described in terms of higher point correlation functions – these are the
Schwinger-Dyson equations in the parlance of quantum field theory. Such
equations for correlators can be then derived from the operator formalism
described here.

To summarise, our goal here is to develop a general framework that can be used
to compute correlation functions and expectation values of system observables
in the simplest way possible. We would also like to describe and compute
Heisenberg picture operators corresponding to these observables without making
any assumptions about the nature of the system.

We will conclude this introduction with an outline: in section §II, we
systematically define a set of reduced Heisenberg operators corresponding to a
particular observable, which we call the image operators. We then introduce
and calculate the reduced $N$-point operators from the image operators, which
help us calculate correlation functions. In section §III, we determine the
time evolution of the image operators and the one point operator. In section
§IV, we show that this set of image operators corresponding to an observable
can be written in terms of the one point operator in a way that depends non-
linearly on the bath state. This gives an alternate prescription for the
computation of $N$-point functions in terms of a deformed product. In section
§V, we examine the Markovian limit, and get a Lindblad equation for 1-point
operators. In section §VI, we show that the number of computations required to
calculate a correlation function can be further reduced. Lastly, in section
§VII, we demonstrate the effectiveness of the techniques developed by solving
a simple spin system.

##  II Setting up the Problem

In any open system, we have a system of interest interacting with a large
bath. Let the Hamiltonian of the system under study be $H_{0}$, and that of
the bath be $H_{B}$. Further, we assume that they interact via the interaction
Hamiltonian $\lambda H_{I}$. Here $\lambda$ is a dimensionless parameter which
defines the strength of coupling. Our goal is to calculate the evolution of
various system observables in the limit that $\lambda$ is very small. The
Hamiltonian of the combined system is: $H=H_{0}+H_{B}+\lambda H_{I}$. This
full system is closed, and does not interact with any external sources. Thus,
any operator $A$ in the Hilbert space of the full system evolves according to
the Heisenberg equation of motion:

$\frac{d}{dt}O(t)=\frac{i}{\hbar}[H(t),O(t)],$ (1)

where we have assumed that none of the operators involved have explicit time
dependence.

Generally, a bath contains many more degrees of freedom than the system under
study. Solving the above equations involves keeping track of the bath as well.
Since we are not interested in the evolution of the bath, we would like to
reduce the problem to the system’s Hilbert space. While doing so, we must make
sure that the expectation values and correlation functions are preserved.

Let us assume that the interaction between the system and the bath was turned
on at $t=0$. Before turning on the interaction, the system and the bath were
decoupled. Therefore, the density matrix of the full system can be written as:
$\rho=\rho_{0}\otimes\rho_{B}$, where $\rho_{0}$ is the system’s density
matrix, and $\rho_{B}$ is the bath’s density matrix. Note that the density
matrix does not evolve in the Heisenberg picture, and therefore we do not have
to worry about the mixing of the system and bath density matrices.

The expectation value of any observable $O$ is calculated as:

$\displaystyle\langle O(t)\rangle$
$\displaystyle=\text{tr}\\{O(t)\rho_{0}\otimes\rho_{B}\\}=\text{tr}_{S}\\{\text{tr}_{B}\\{O(t)\rho_{B}\\}\rho_{0}\\}=\text{tr}_{S}\\{O_{S}(t)\rho_{0}\\},$
(2)

wherewe have defined $O_{S}(t)\equiv\text{tr}_{B}\\{O(t)\rho_{B}\\}$. Note
that the operator $O_{S}(t)$ is confined to the system’s Hilbert space, since
we have traced out the bath degrees of freedom. Further, we can calculate the
expectation value of the observable $O$ by taking a trace of
$O_{S}(t)\rho_{0}$ over the system. Thus, this operator can be thought of as
the reduced Heisenberg operator corresponding to the observable $O$, as far as
one point functions are concerned. We call this operator a 1-point operator.

Similarly, we can define an $N$-point operator as:

$(O_{1}(t_{1})\dots O_{N}(t_{N}))_{S}=\text{tr}_{B}\\{O_{1}(t_{1})\dots
O_{N}(t_{N})\rho_{B}\\}.$ (3)

However, we run into problems while multiplying two reduced operators. We note
that $(O_{1}(t_{1})O_{2}(t_{2}))_{S}\neq O_{1S}(t_{1})O_{2S}(t_{2})$ due to
the fact that the trace does not distribute over multiplication. Thus, the
information about the operator $(O_{1}(t_{1})O_{2}(t_{2}))_{S}$ is not fully
contained inside $O_{1S}(t_{1})$ and $O_{2S}(t_{2})$. This substantiates our
assertion in the introduction: the specification of $O_{S}(t)$ alone is not a
complete Heisenberg picture description of the dynamics since it misses the
information about higher point correlations.

To be able to construct any general $N$-point reduced operators, we will now
define another set of intermediate operators. Consider the Hilbert spaces
$\mathcal{H}_{0}$ and $\mathcal{H}_{B}$ of the system and the bath
respectively. Let $\\{\ket{i}|\ i\in I\\}$ be an orthonormal basis of
$\mathcal{H}_{0}$ and $\\{\ket{\alpha}|\ \alpha\in A\\}$ be an orthonormal
basis of $\mathcal{H}_{B}$. Thus, the full Hilbert space
$\mathcal{H}=\mathcal{H}_{0}\otimes\mathcal{H}_{B}$ is spanned by the basis
$\\{\ket{i\alpha}|\ i\in I,\alpha\in A\\}$.

Next, we define the projection operator
$T_{\alpha}:\mathcal{H}_{0}\rightarrow\mathcal{H}$ as
$T_{\alpha}=\ket{i\alpha}\bra{i}$; where we have used the Einstein summation
convention for summing over $i$. Note that this projection operator maps
states in the system’s Hilbert space $\mathcal{H}_{0}$ to states in the full
Hilbert space $\mathcal{H}_{0}\otimes\mathcal{H}_{B}$. These projection
operators obey the identity:

$T_{\alpha}T_{\alpha}^{\dagger}=\mathbbm{1}.$ (4)

Here, the symbol $\mathbbm{1}$ denotes the identity operator in the full
Hilbert space $\mathcal{H}_{0}\otimes\mathcal{H}_{B}$. Finally, we define the
operator $O_{\alpha\beta}(t)$ as:

$O_{\alpha\beta}(t)=T_{\alpha}^{\dagger}O(t)T_{\beta}.$ (5)

We call these operators as the _image operators_ $O_{\alpha\beta}(t)$ of the
observable $O$. These image operators act on the system Hilbert space and the
number of such image operators is equal to the dimension of the
bath/environment.

We can construct the 2-point image operators
$(O_{1}(t_{1})O_{2}(t_{2}))_{\alpha\beta}$ from the 1-point image operators
$O_{1\alpha\beta}(t_{1})$ and $O_{2\alpha\beta}(t_{2})$ using (4):

$(O_{1}(t_{1})O_{2}(t_{2}))_{\alpha\beta}=O_{1\alpha\gamma}(t_{1})O_{2\gamma\beta}(t_{2}).$
(6)

Similarly one can construct all $N$-point image operators and from them the
$N$-point reduced operators:

$(O_{1}(t_{1})\dots O_{N}(t_{N}))_{S}=(O_{1}(t_{1})\dots
O_{N}(t_{N}))_{\alpha\beta}\rho_{B\beta\alpha}.$ (7)

Inserting the identity (4) between operators in the Heisenberg equation (1),
we obtain the reduced Heisenberg equation for the image operators:

$\frac{d}{dt}O_{\alpha\beta}(t)=\frac{i}{\hbar}\Big{\\{}H_{\alpha\gamma}(t)O_{\gamma\beta}(t)-O_{\alpha\gamma}(t)H_{\gamma\beta}(t)\Big{\\}}.$
(8)

Note that this is an exact equation for a general open quantum system. We have
thus established that the image operators give a complete description of the
open quantum system in the Heisenberg picture, allowing us to compute
arbitrary higher point functions pf the system observables. The coupled
evolution of the image operators can also be then described exactly by the
above equation, provided the image Hamiltonians describing the open quantum
system are known.

##  III Perturbative Solution for System Operators

In this section, we will solve the reduced Heisenberg equation for a general
system operator in a perturbative expansion in the system bath coupling.As is
customary, we will do this by passing on to an interaction picture.

Before $t=0$, the system and bath do not interact and the system operator is
confined to the system’s Hilbert space,i.e.,

$O(0)=O\otimes\mathbbm{1}_{B},$ (9)

where $O$ is a Schrodinger picture operator and $\mathbbm{1}_{B}$ is the
identity operator in the bath Hilbert space. The image operators corresponding
to $O$ are then given by: $O_{\alpha\beta}(0)=O\ \delta_{\alpha\beta}$ at
$t=0$. Once the interaction is turned on, these image operators evolve and
become non-zero for $\alpha\neq\beta$.

From (1), it is easy to see that the Hamiltonian of the entire system, $H(t)$,
does not evolve in time:

$H(t)=H(0)=H_{0}+H_{B}+\lambda H_{I},$ (10)

where, $H_{0}$, $H_{B}$, and $H_{I}$ are Schrodinger picture operators. Thus
the reduced Heisenberg equation becomes:

$\frac{d}{dt}O_{\alpha\beta}(t)=\frac{i}{\hbar}\Big{\\{}H_{\alpha\gamma}O_{\gamma\beta}(t)-O_{\alpha\gamma}(t)H_{\gamma\beta}\Big{\\}},$
(11)

where $H_{\alpha\beta}$s are Schrodinger picture operators.

The Hamiltonian of the system $H_{0}$ is confined to the Hilbert space
$\mathcal{H}_{0}$ and hence the corresponding image operators are given by
$H_{0\alpha\beta}=H_{0}\ \delta_{\alpha\beta}$. On the other hand, the bath
Hamiltonian has the image operators $H_{B\alpha\beta}=H_{B\alpha\beta}\
\mathbbm{1}_{0}$ with $\mathbbm{1}_{0}$ being the identity operator in
$\mathcal{H}_{0}$ and $H_{B\alpha\beta}$ being the $(\alpha,\beta)^{th}$
element of $H_{B}$. Further, for simplicity, let us assume that the basis
$\\{\ket{\alpha}\\}$ of the Hilbert space $\mathcal{H}_{B}$ is the energy
eigen-basis of $H_{B}$ with $H_{B}\ket{\alpha}=E_{\alpha}\ket{\alpha}$. Then
the reduced Heisenberg equation takes the form

$\displaystyle\frac{d}{dt}O_{\alpha\beta}(t)=$
$\displaystyle\frac{i}{\hbar}[H_{0},O_{\alpha\beta}(t)]+\frac{i}{\hbar}(E_{\alpha}-E_{\beta})O_{\alpha\beta}(t)+\frac{i\lambda}{\hbar}\Big{\\{}H_{I\alpha\gamma}O_{\gamma\beta}(t)-O_{\alpha\gamma}(t)H_{I\gamma\beta}\Big{\\}}.$
(12)

We can pass on to the interaction picture by extending the standard steps to
image operators. We define the interaction picture operators
$\tilde{O}_{\alpha\beta}(t)$ and $\tilde{H}_{I\alpha\beta}(t)$ as:

$\tilde{O}_{\alpha\beta}(t)\equiv U_{0}(t)\ O_{\alpha\beta}(t)\
U_{0}^{\dagger}(t)\ e^{-\tfrac{i(E_{\alpha}-E_{\beta})t}{\hbar}},$ (13a)
$\tilde{H}_{I\alpha\beta}(t)\equiv U_{0}(t)\ H_{I\alpha\beta}\
U_{0}^{\dagger}(t)\ e^{-\tfrac{i(E_{\alpha}-E_{\beta})t}{\hbar}},$ (13b)

where $U_{0}(t)=e^{-\tfrac{iH_{0}t}{\hbar}}$ is the system evolution operator
before environmental interactions are turned on. Inserting these into the
above equation, we get:

$\frac{d}{dt}\tilde{O}_{\alpha\beta}(t)=\frac{i\lambda}{\hbar}\Big{\\{}\tilde{H}_{I\alpha\gamma}(t)\tilde{O}_{\gamma\beta}(t)-\tilde{O}_{\alpha\gamma}(t)\tilde{H}_{I\gamma\beta}(t)\Big{\\}}.$
(14)

This is the reduced Heisenberg equation in the interaction picture. This
equation can now be solved perturbatively in $\lambda$.

Up to first order in $\lambda$, we obtain

$\tilde{O}_{\alpha\beta}(t)=O\
\delta_{\alpha\beta}+\frac{i\lambda}{\hbar}\int_{0}^{t}dt_{1}[\tilde{H}_{I\alpha\beta}(t_{1}),O]+\mathcal{O}(\lambda^{2}).$
(15)

A solution to arbitrary order in $\lambda$ can be obtained via the Dyson
series:

$\tilde{O}_{\alpha\beta}(t)=\tilde{U}^{\dagger}_{I\gamma\alpha}(t)\ O\
\tilde{U}_{I\gamma\beta}(t),$ (16)

where the evolution operator $\tilde{U}_{I\alpha\beta}(t)$ has the
perturbative expansion

$\tilde{U}_{I\alpha\beta}(t)=\sum_{n=0}^{\infty}\Big{(}\frac{-i\lambda}{\hbar}\Big{)}^{n}\
\tilde{K}^{(n)}_{\alpha\beta}(t).$ (17)

Here, $\tilde{K}^{(n)}$ is the time-ordered product of interaction
Hamiltonians, as is familar from time-dependent perturbation theory:

$\tilde{K}^{(n)}_{\alpha\beta}(t)=\int_{0}^{t}dt_{1}\dots\int_{0}^{t_{n-1}}dt_{n}\
\tilde{H}_{I\alpha\gamma_{1}}(t_{1})\dots\tilde{H}_{I\gamma_{n-1}\beta}(t_{n}),$
(18)

In particular, we have
$\tilde{K}^{(0)}_{\alpha\beta}(t)=\delta_{\alpha\beta}$. We can now revert
back to the Heisenberg picture to write down a perturbative expansion for the
image operators

$\displaystyle O_{\alpha\beta}(t)$
$\displaystyle=U^{\dagger}_{I\gamma\alpha}(t)\ U_{0}^{\dagger}(t)\ O\
U_{0}(t)\ U_{I\gamma\beta}(t)$
$\displaystyle=\sum_{n=0}^{\infty}\Big{(}\frac{\lambda}{\hbar}\Big{)}^{n}\sum_{r=0}^{n}i^{n-2r}\
K^{(n-r)\dagger}_{\gamma\alpha}(t)U_{0}^{\dagger}(t)OU_{0}(t)K^{(r)}_{\gamma\beta}(t),$
(19)

where we have defined

$K^{(n)}_{\alpha\beta}(t)\equiv
e^{\frac{i(E_{\alpha}-E_{\beta})t}{\hbar}}U_{0}^{\dagger}(t)\tilde{K}^{(n)}_{\alpha\beta}(t)U_{0}(t).$
(20)

In order to avoid clutter in our subsequent expressions, we define the super-
operators $P^{(n)}_{\alpha\beta}$, and $P_{S}^{(n)}$ acting on a system
operator $A(t)$ as:

$P^{(n)}_{\alpha\beta}A(t)\equiv\sum_{r=0}^{n}i^{n-2r}K^{(n-r)\dagger}_{\gamma\alpha}(t)A(t)K^{(r)}_{\gamma\beta}(t),$
(21)

and

$P^{(n)}_{S}A(t)\equiv P^{(n)}_{\alpha\beta}A(t)\rho_{\beta\alpha}.$ (22)

Contracting $O_{\alpha\beta}$ with the bath density matrix
$\rho_{B\beta\alpha}$ gives a compact expression for $O_{S}(t)$ in terms of
these super-operators:

$\displaystyle O_{S}(t)=$
$\displaystyle\sum_{n=0}^{\infty}\Big{(}\frac{\lambda}{\hbar}\Big{)}^{n}P^{(n)}_{S}U_{0}^{\dagger}(t)OU_{0}(t).$
(23)

##  IV The Deformed Operator Product : Image Operators In Terms of One Point
Operators

In this section, we show that all $N$-point operators can be expressed solely
in terms of one point operators. We will then argue that these expressions can
equivalently be interpreted as deforming the operator product in the Hilbert
space of the open quantum system.

Using the super operators defined in the last section, the required
expressions can be obtained as follows. First, we rewrite the equation (23) in
the form

$\displaystyle O_{S}(t)=$ $\displaystyle\
\Big{\\{}1+\sum_{n=1}^{\infty}\Big{(}\frac{\lambda}{\hbar}\Big{)}^{n}P^{(n)}_{S}\Big{\\}}\
U_{0}^{\dagger}(t)\ O\ U_{0}(t),$ (24)

and then invert it via multinomial expansion to give

$\begin{split}U_{0}^{\dagger}(t)\ O\ U_{0}(t)&=\
\Big{\\{}1+\sum_{n=1}^{\infty}\Big{(}\frac{\lambda}{\hbar}\Big{)}^{n}P^{(n)}_{S}\Big{\\}}^{-1}O_{S}(t)\\\
&=\sum_{k=0}^{\infty}\sum_{n_{1},\dots,n_{k}=1}^{\infty}(-1)^{k}\left(\frac{\lambda}{\hbar}\right)^{n_{1}+\dots+n_{k}}P^{(n_{1})}_{S}\dots
P^{(n_{k})}_{S}O_{S}(t).\end{split}$ (25)

In the second step, we can use the equation (19) to express
$O_{\alpha\beta}(t)$ in terms of $O_{S}(t)$:

$\displaystyle
O_{\alpha\beta}(t)=\sum_{n=0}^{\infty}\sum_{k=0}^{\infty}\sum_{n_{1},\dots,n_{k}=1}^{\infty}(-1)^{k}\left(\frac{\lambda}{\hbar}\right)^{n+n_{1}+\dots+n_{k}}P^{(n)}_{\alpha\beta}P^{(n_{1})}_{S}\dots
P^{(n_{k})}_{S}O_{S}(t).$ (26)

This is the central result of this note. As advertised, the above expression
shows that all the image operators can be expressed in terms of one point
operators.

As all $N$ point operators can be reconstructed using image operators, it
follows that reduced $N$ point operators can also be written in terms of the
one point operators. However, as we noted before, the reduced $N$ point
operators are not simple products of the one point operator. The one point
operators have to be dressed into the image operators using the above formula
and then multiplied at the level of image operators with appropriate index
contractions. This procedure can be summarised by saying that the rule for
operator product on the system Hilbert space has to be modified.

Let us define a new operator product on the Hilbert space of an open quantum
state such that

$\begin{split}O_{S}(t_{1})\star O_{S}(t_{2})\ldots\star O_{S}(t_{N})\equiv
O_{\alpha_{1}\alpha_{2}}(t_{1})\ O_{\alpha_{2}\alpha_{3}}(t_{2})\ldots
O_{\alpha_{N}\alpha_{N+1}}(t_{N})\ \rho_{B\alpha_{N+1}\alpha_{N}}\
.\end{split}$ (27)

The way to interpret this expression is to think of the RHS using the equation
(26) which then defines a deformed operator product on the system Hilbert
space. This shows that the new product depends non-trivially on the state of
the bath. Once such a deformed operator product is defined, we can, in
principle, remove all references to the image operators at the cost of working
with a complicated product.

One can further write down an equation that expresses the evolution of all
$N$-point operators in terms of 1-point operators, that is local in time. For
example, $\frac{d}{dt}O_{S}(t)$ can be expressed as:

$\displaystyle\frac{d}{dt}O_{S}(t)=$
$\displaystyle\frac{i}{\hbar}[H_{0},O_{S}(t)]$
$\displaystyle+\sum_{n=1}^{\infty}\sum_{k=0}^{\infty}\sum_{n_{1},\dots,n_{k}=1}^{\infty}(-1)^{k}\Big{(}\frac{\lambda}{\hbar}\Big{)}^{n+n_{1}+\dots+n_{k}}\mathcal{D}_{t}P^{(n)}_{S}P^{(n_{1})}_{S}\dots
P^{(n_{k})}_{S}O_{S}(t),$ (28)

where, the super-operator $\mathcal{D}_{t}P^{(n)}_{S}$ acting on an operator
$A(t)$ is defined as:

$\displaystyle\mathcal{D}_{t}P^{(n)}_{S}A(t)=$
$\displaystyle\sum_{r=0}^{n}i^{n-2r}U_{0}^{\dagger}(t)\frac{d}{dt}\Big{[}U_{0}(t)K^{(n-r)\dagger}_{\gamma\alpha}(t)U_{0}^{\dagger}(t)\Big{]}U_{0}(t)A(t)K^{(r)}_{\gamma\beta}(t)\rho_{\beta\alpha}$
$\displaystyle+\sum_{r=0}^{n}i^{n-2r}K^{(n-r)\dagger}_{\gamma\alpha}(t)A(t)U_{0}^{\dagger}(t)\frac{d}{dt}\Big{[}U_{0}(t)K^{(r)}_{\gamma\beta}(t)U_{0}^{\dagger}(t)\Big{]}U_{0}(t)\rho_{\beta\alpha}.$
(29)

##  V Markovian Limit

A Markovian system is one in which the system does not remember its previous
state. In the Schordinger picture, the density matrix obeys the Lindblad
equation in the Markovian limit [13, 10, 12]. Under certain approximations,
the reduced one point operator obeys a similar equation, called the adjoint
master equation [3]. We would like to understand which assumptions lead to a
Lindblad equation in our formalism.

A general interaction Hamiltonian can be written as:

$H_{I}=\sum_{i}R^{i}\otimes S^{i},$ (30)

where $R^{i}$ is an operator in the system’s Hilbert space, and $S^{i}$ is an
operator in the bath’s Hilbert space. Thus, the corresponding image operators
of $H_{I}$ are:

$H_{I\alpha\beta}=\sum_{i}R^{i}S^{i}_{\alpha\beta}.$ (31)

Note that since $S^{i}$ is confined to $\mathcal{H}_{B}$,
$S^{i}_{\alpha\beta}$ is the $\alpha$, $\beta$ th element of $S^{i}$. From
(13b), we see that:

$\tilde{H}_{I\alpha\beta}(t)=\sum_{i}\tilde{R}^{i}(t)\tilde{S}^{i}_{\alpha\beta}(t),$
(32)

where, $\tilde{R}^{i}(t)=U_{0}(t)R^{i}U_{0}^{\dagger}(t)$, and
$\tilde{S}^{i}_{\alpha\beta}(t)={S}^{i}_{\alpha\beta}\
e^{-\frac{i(E_{\alpha}-E_{\beta})t}{\hbar}}$. We claim that $O_{S}(t)$ obeys a
Lindblad equation under the following assumptions:

  * •

The first moment of the interaction Hamiltonian is zero. That is,
$\text{tr}_{B}\\{\tilde{H}_{I}(t)\rho_{B}\\}=0$.

  * •

The two point bath correlation functions
$\langle\tilde{S}^{i}(t)\tilde{S}^{j}(t-\tau)\rangle_{B}$ decay rapidly with
$\tau$.

  * •

The two point bath correlation functions
$\langle\tilde{S}^{i}(t)\tilde{S}^{j}(t-\tau)\rangle_{B}$ are independent of
$t$.

The first assumption implies that:

$\tilde{S}^{i}_{\alpha\beta}(t)\rho_{B\beta\alpha}=0.$ (33)

Differentiating this with respect to $t$, we get:

$\sum_{\alpha\beta}(E_{\alpha}-E_{\beta})\tilde{S}^{i}_{\alpha\beta}(t)\rho_{B\beta\alpha}=0.$
(34)

Similarly, the third assumption implies that:

$\sum_{\alpha\beta}(E_{\alpha}-E_{\beta})\tilde{S}^{i}_{\alpha\gamma}(t)\tilde{S}^{i}_{\gamma\beta}(t-\tau)\rho_{B\beta\alpha}=0.$
(35)

Using these assumptions, equation (28) simplifies to:

$\displaystyle\frac{d}{dt}O_{S}(t)=$ $\displaystyle\
\frac{i}{\hbar}H_{0}O_{S}(t)+\Big{(}\frac{i\lambda}{\hbar}\Big{)}^{2}\sum_{\omega,\omega^{\prime}}\sum_{ij}J^{ij}(\omega)\\{A^{i\dagger}_{\omega}A^{j}_{\omega^{\prime}}O_{S}(t)-A^{i\dagger}_{\omega}O_{S}(t)A^{j}_{\omega^{\prime}}\\}+\text{h.c.}+\mathcal{O}(\lambda^{3}),$
(36)

where,

$J^{ij}(\omega)=\int_{0}^{\infty}d\tau\
e^{-i\omega\tau}\tilde{S}^{i}_{\alpha\gamma}(0)\tilde{S}^{j}_{\gamma\beta}(-\tau)\rho_{B\beta\alpha},$
(37)

and, $A^{i}$s are the Fourier coefficients of $\tilde{R}^{i}(t)$:

$\tilde{R}^{i}(t)=\sum_{\omega}e^{i\omega t}A^{i}_{\omega}.$ (38)

The above equation is of the Lindblad form.

##  VI Structure of N-Point Operators

In section §IV, we showed that $N$-point operators can be constructed from one
point operators, which reduces the number of calculations one has to perform
to calculate $N$-point correlation functions. However, one can further reduce
the number of calculations required to compute an $N$-point operator. One
finds that $N$-point operators have a very peculiar structure which can be
exploited to reduce the complexity of any calculation.

Using the definitions (21) and (22), we can expand out equation (26). We
describe a systematic way of writing down the $n$th order term in this
expansion: We partition $n$ into an even number of partitions
$\\{n_{1},m_{1},\dots,n_{k},m_{k}\\}$, such that:

$\sum_{i=1}^{k}(n_{i}+m_{i})=n,$ (39a) $n_{1}+m_{1}\geq 0,$ (39b)
$n_{i}+m_{i}>0,\ i=2,\dots,k.$ (39c)

We construct the term corresponding to this partition as follows: Begin with
$n_{k}$ and $m_{k}$. We multiply $O_{S}(t)$ by the operator
$\Big{[}\Big{(}\frac{-i\lambda}{\hbar}\Big{)}^{n_{k}}K^{(n_{k})}_{\gamma_{k}\alpha_{k}}(t)\Big{]}^{\dagger}$
from the left, and by the operator
$\Big{[}\Big{(}\frac{-i\lambda}{\hbar}\Big{)}^{m_{k}}K^{(m_{k})}_{\gamma_{k}\beta_{k}}(t)\Big{]}$
from the right. Here, a sum over $\gamma_{k}$ is implied. We contract this
with the bath density matrix $\rho_{B\beta_{k}\alpha_{k}}$. We do a similar
operation for $n_{k-1}$ and $m_{k-1}$. We then continue in a similar fashion
until $n_{2}$ and $m_{2}$. For $n_{1}$ and $m_{1}$, we perform the same
operation, except that we do not contract with the bath density matrix this
time. At the end, we add a factor of $(-1)^{k-1}$. Finally, we get the term:

$\begin{split}&(-1)^{k-1}\Big{[}\Big{(}\frac{-i\lambda}{\hbar}\Big{)}^{n_{1}}K^{(n_{1})}_{\gamma_{1}\alpha}(t)\Big{]}^{\dagger}\Big{[}\Big{(}\frac{-i\lambda}{\hbar}\Big{)}^{n_{2}}K^{(n_{2})}_{\gamma_{2}\alpha_{2}}(t)\Big{]}^{\dagger}\dots\Big{[}\Big{(}\frac{-i\lambda}{\hbar}\Big{)}^{n_{k}}K^{(n_{k})}_{\gamma_{k}\alpha_{k}}(t)\Big{]}^{\dagger}O_{S}(t)\\\
&\Big{[}\Big{(}\frac{-i\lambda}{\hbar}\Big{)}^{m_{k}}K^{(m_{k})}_{\gamma_{k}\beta_{k}}(t)\Big{]}\dots\Big{[}\Big{(}\frac{-i\lambda}{\hbar}\Big{)}^{m_{2}}K^{(m_{2})}_{\gamma_{2}\beta_{2}}(t)\Big{]}\Big{[}\Big{(}\frac{-i\lambda}{\hbar}\Big{)}^{m_{1}}K^{(m_{1})}_{\gamma_{1}\beta}(t)\Big{]}\rho_{B\beta_{2}\alpha_{2}}\dots\rho_{B\beta_{k}\alpha_{k}},\end{split}$
(40)

corresponding to the partition $\\{n_{1},m_{1},\dots,n_{k},m_{k}\\}$. We sum
over all such terms corresponding to even partitions of $n$ to get the $n$th
order term in the expansion of the image operator $O_{\alpha\beta}(t)$.

We can make a quick check of the validity of this method. When
$O_{\alpha\beta}(t)$ is contracted with $\rho_{B\beta\alpha}$, we get
$O_{S}(t)$ by definition. We expect the above method to give us the same
result. Consider the two partitions of $n$:
$\\{n_{1},m_{1},\dots,n_{k},m_{k}\\}$, and
$\\{0,0,n_{1},m_{1},\dots,n_{k},m_{k}\\}$, such that $n_{1}+m_{1}>0$. The
terms in the expansion of $O_{\alpha\beta}(t)$ corresponding to these
partitions can be calculated by using (40). If we contract both these terms
with $\rho_{B\beta\alpha}$, then they differ only by a negative sign. Thus,
both of them cancel out each other. The only remaining term is the one
corresponding to the partition $\\{0,0\\}$. Contracting this with
$\rho_{B\beta\alpha}$ gives us $O_{S}(t)$, as expected.

For an $N$ point operator, we define its irreducible part,
$I[O_{1}(t_{1}),\dots,O_{N}(t_{N})]$ as the sum of the terms which cannot be
broken into simpler terms. Diagrammatically, we represent this irreducible
part by:

$I[O_{1}(t_{1}),\dots,O_{N}(t_{N})]=\leavevmode\hbox to51.61pt{\vbox
to51.61pt{\pgfpicture\makeatletter\hbox{\hskip 25.8073pt\lower-25.8073pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\definecolor{pgffillcolor}{rgb}{0.85,0.85,0.85}\pgfsys@color@gray@fill{0.85}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{25.6073pt}{0.0pt}\pgfsys@curveto{25.6073pt}{14.1427pt}{14.1427pt}{25.6073pt}{0.0pt}{25.6073pt}\pgfsys@curveto{-14.1427pt}{25.6073pt}{-25.6073pt}{14.1427pt}{-25.6073pt}{0.0pt}\pgfsys@curveto{-25.6073pt}{-14.1427pt}{-14.1427pt}{-25.6073pt}{0.0pt}{-25.6073pt}\pgfsys@curveto{14.1427pt}{-25.6073pt}{25.6073pt}{-14.1427pt}{25.6073pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{
{}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} {
}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-16.11113pt}{-2.5pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{1, \ldots, N}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{
{}{}{}{}{}}{{{}}{{}}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>(41)

We claim that all time ordered $N$ point operators can be broken down into a
sum of irreducible operators. We demonstrate this by analysing the two point
and three point operators. The two point image operators can be constructed
from the one point image operators as:

$(O_{1}(t_{1})O_{2}(t_{2}))_{\alpha\beta}=O_{1\alpha\gamma}(t_{1})O_{2\gamma\beta}(t_{2}).$
(42)

Thus, the $n$th order term in the two point image operator is formed by the
$2k_{1}+2k_{2}$ partition of $n$. We will denote this partition by
$\\{n^{(1)}_{1},m^{(1)}_{1},\dots,n^{(1)}_{k_{1}},m^{(1)}_{k_{1}};n^{(2)}_{1},m^{(2)}_{1},\dots,n^{(2)}_{k_{2}},m^{(2)}_{k_{2}}\\}$.
As before, this partition must satisfy the conditions (39b) and (39c). We
contract each partition with $\rho_{B\beta\alpha}$ to get the expansion of the
two point operator, $(O_{1}(t_{1})O_{2}(t_{2}))_{S}$. Table 1 lists all
possible non-zero partitions and their corresponding diagrams. All other
partitions cancel out each other. This gives us the relation:

$(O_{1}(t_{1})O_{2}(t_{2}))_{S}=O_{1S}(t_{1})O_{2S}(t_{2})+I[O_{1S}(t_{1}),O_{2S}(t_{2})].$
(43)

Note that this means that the irreducible part of the two point operator is
nothing but the second order cumulant:

$I[O_{1S}(t_{1}),O_{2S}(t_{2})]=(O_{1}(t_{1})O_{2}(t_{2}))_{S}-O_{1S}(t_{1})O_{2S}(t_{2}).$
(44)

Partition | Term | Diagram
---|---|---
$n^{(1)}_{i}=m^{(1)}_{i}=n^{(2)}_{i}=m^{(2)}_{i}=0$ | $O_{1S}(t_{1})O_{2S}(t_{2})$ | 12
$n^{(1)}_{1}+m^{(1)}_{1}>0$, $n^{(2)}_{1}+m^{(2)}_{1}>0$ | $I[O_{1S}(t_{1}),O_{2S}(t_{2})]$ | 1, 2
Table 1: Two Point Operator

The three point image operators can be constructed from the one point image
operators as:

$(O_{1}(t_{1})O_{2}(t_{2})O_{3}(t_{3}))_{\alpha\beta}=O_{1\alpha\gamma}(t_{1})O_{2\gamma\delta}(t_{2})O_{3\delta\beta}(t_{3}).$
(45)

Again, as before, the $n$th order term in the three point image operator is
formed by $2k_{1}+2k_{2}+2k_{3}$ partitions of $n$. They will be denoted by
$\\{n^{(1)}_{1},m^{(1)}_{1},\dots,n^{(1)}_{k_{1}},m^{(1)}_{k_{1}};n^{(2)}_{1},m^{(2)}_{1},\dots,\\\
n^{(2)}_{k_{2}},m^{(2)}_{k_{2}};n^{(3)}_{1},m^{(3)}_{1},\dots,n^{(3)}_{k_{3}},m^{(3)}_{k_{3}}\\}$.
To get the three point operator, $(O_{1}(t_{1})O_{2}(t_{2})O_{3}(t_{3}))_{S}$,
we contract the three point image operators with $\rho_{B\beta\alpha}$. We
list all partitions and their corresponding terms and diagrams in table 2.

Partition | Term | Diagram
---|---|---
$n^{(k)}_{i}=m^{(k)}_{i}=0$, $k=1,2,3$ | $O_{1S}(t_{1})O_{2S}(t_{2})O_{3S}(t_{3})$ | 123
$n^{(1)}_{1}+m^{(1)}_{1}>0$, $n^{(2)}_{1}+m^{(2)}_{1}>0$, $n^{(3)}_{i}=m^{(3)}_{i}=0$ | $\mathcal{W}_{1,2,3}\\{O_{3S}(t_{3})I[O_{1S}(t_{1}),O_{2S}(t_{2})]\\}$ | 31, 2
$n^{(1)}_{1}+m^{(1)}_{1}>0$, $n^{(2)}_{i}=m^{(2)}_{i}=0$, $n^{(3)}_{1}+m^{(3)}_{1}>0$ | $\mathcal{W}_{1,2,3}\\{O_{2S}(t_{2})I[O_{3S}(t_{3}),O_{1S}(t_{1})]\\}$ | 23, 1
$n^{(1)}_{i}=m^{(1)}_{i}=0$, $n^{(2)}_{1}+m^{(2)}_{1}>0$, $n^{(3)}_{1}+m^{(3)}_{1}>0$ | $\mathcal{W}_{1,2,3}\\{O_{1S}(t_{1})I[O_{2S}(t_{2}),O_{3S}(t_{3})]\\}$ | 12, 3
$n^{(1)}_{1}+m^{(1)}_{1}>0$, $n^{(2)}_{1}+m^{(2)}_{1}>0$, $n^{(3)}_{1}+m^{(3)}_{1}>0$ | $I[O_{1S}(t_{1}),O_{2S}(t_{2}),O_{3S}(t_{3})]$ | 1, 2, 3
Table 2: Three Point Operator

Note that we have introduced the operator $\mathcal{W}_{1,2,3}$. This operator
makes sure that the operator product is ordered such that $O_{1}$ comes before
$O_{2}$, and $O_{2}$ comes before $O_{3}$. For example, the third term in
table 2 is:

$\mathcal{W}_{1,2,3}\\{O_{2S}(t_{2})I[O_{3S}(t_{3}),O_{1S}(t_{1})]\\}=(O_{1}(t_{1})O_{2S}(t_{2})O_{3}(t_{3}))_{S}-O_{1S}(t_{1})O_{2S}(t_{2})O_{3S}(t_{3}).$
(46)

Once again, we find that the irreducible part of the three point operator is
the third order cumulant:

$\displaystyle I[O_{1S}(t_{1}),O_{2S}(t_{2}),O_{3S}(t_{3})]=$ $\displaystyle\
(O_{1}(t_{1})O_{2}(t_{2})O_{3}(t_{3}))_{S}$
$\displaystyle-[(O_{1}(t_{1})O_{2}(t_{2})O_{3S}(t_{3}))_{S}-O_{1S}(t_{1})O_{2S}(t_{2})O_{3S}(t_{3})]$
$\displaystyle-[(O_{1}(t_{1})O_{2S}(t_{2})O_{3}(t_{3}))_{S}-O_{1S}(t_{1})O_{2S}(t_{2})O_{3S}(t_{3})]$
$\displaystyle-[(O_{1S}(t_{1})O_{2}(t_{2})O_{3}(t_{3}))_{S}-O_{1S}(t_{1})O_{2S}(t_{2})O_{3S}(t_{3})]$
$\displaystyle-O_{1S}(t_{1})O_{2S}(t_{2})O_{3S}(t_{3}).$ (47)

##  VII Example: Two Qubit System

We demonstrate the usefulness of the methods developed above by solving a
simple system. Consider a system of two spin-half particles with spin
operators $\overline{S}_{1}$ and $\overline{S}_{2}$ interacting via the
Hamiltonian:

$H=\lambda\ \overline{S}_{1}\cdot\overline{S}_{2}.$ (48)

We consider the first particle to be our system and the second particle to be
the bath. We will also assume that the bath has a density matrix of the form:

$\displaystyle\rho_{B}$
$\displaystyle=(1-c)\ket{\uparrow}\bra{\uparrow}+c\ket{\downarrow}\bra{\downarrow}=\begin{pmatrix}1-c&0\\\
0&c\end{pmatrix}.$ (49)

This system can be solved exactly, and hence it is easy to verify the results
of our perturbative method.

The interaction Hamiltonian of the system is
$H_{I}=\overline{S}_{1}\cdot\overline{S}_{2}$. Thus, the image operators of
$H_{I}$ are:

$\displaystyle H_{I\uparrow\uparrow}=\frac{\hbar^{2}}{4}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix},H_{I\uparrow\downarrow}=\frac{\hbar^{2}}{2}\begin{pmatrix}0&0\\\
1&0\end{pmatrix},H_{I\downarrow\uparrow}=\frac{\hbar^{2}}{2}\begin{pmatrix}0&1\\\
0&0\end{pmatrix},H_{I\downarrow\downarrow}=\frac{\hbar^{2}}{4}\begin{pmatrix}-1&0\\\
0&1\end{pmatrix}.$ (50)

Using this, we can calculate $K^{(1)}_{\alpha\beta}(t)$ and
$K^{(2)}_{\alpha\beta}(t)$:

$\displaystyle
K_{\uparrow\uparrow}^{(1)}(t)=\frac{\hbar^{2}t}{4}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix},K_{\uparrow\downarrow}^{(1)}(t)=\frac{\hbar^{2}t}{2}\begin{pmatrix}0&0\\\
1&0\end{pmatrix},$ $\displaystyle
K_{\downarrow\uparrow}^{(1)}(t)=\frac{\hbar^{2}t}{2}\begin{pmatrix}0&1\\\
0&0\end{pmatrix},K_{\downarrow\downarrow}^{(1)}(t)=\frac{\hbar^{2}t}{4}\begin{pmatrix}-1&0\\\
0&1\end{pmatrix},$ (51) $\displaystyle
K_{\uparrow\uparrow}^{(2)}(t)=\frac{\hbar^{4}t^{2}}{32}\begin{pmatrix}1&0\\\
0&5\end{pmatrix},K_{\uparrow\downarrow}^{(2)}(t)=\frac{\hbar^{4}t^{2}}{8}\begin{pmatrix}0&0\\\
-1&0\end{pmatrix},$ $\displaystyle
K_{\downarrow\uparrow}^{(2)}(t)=\frac{\hbar^{4}t^{2}}{8}\begin{pmatrix}0&-1\\\
0&0\end{pmatrix},K_{\downarrow\downarrow}^{(2)}(t)=\frac{\hbar^{4}t^{2}}{32}\begin{pmatrix}5&0\\\
0&1\end{pmatrix}.$ (52)

We can contract $K^{(1)}_{\alpha\beta}(t)$ and $K^{(2)}_{\alpha\beta}(t)$ with
the density matrix $\rho_{B\beta\alpha}$ to obtain $K^{(1)}_{S}(t)$ and
$K^{(2)}_{S}(t)$:

$K^{(1)}_{S}(t)=(1-2c)\frac{\hbar^{2}t}{4}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix},$ (53a)
$K^{(2)}_{S}(t)=\frac{\hbar^{4}t^{2}}{32}\begin{pmatrix}1+4c&0\\\
0&5-4c\end{pmatrix}.$ (53b)

Let us now use our method to calculate the operator $S_{1xS}(t)$, which is the
$x$-component of the first particle’s spin. We will be using equation (23). Up
to second order, we have:

$\displaystyle S_{1xS}(t)=$ $\displaystyle\
S_{1x}+\frac{i\lambda}{\hbar}[K^{(1)}_{S}(t),S_{1x}]+\Big{(}\frac{i\lambda}{\hbar}\Big{)}^{2}\Big{\\{}K^{(2)\dagger}_{S}(t)S_{1x}+S_{1x}K^{(2)}_{S}(t)$
$\displaystyle-K^{(1)\dagger}_{\gamma\alpha}(t)S_{1x}K^{(1)}_{\gamma\beta}(t)\rho_{B\beta\gamma}\Big{\\}}+\mathcal{O}(\lambda^{3})$
$\displaystyle=$
$\displaystyle\frac{\hbar}{2}\begin{pmatrix}0&1+\frac{1}{2}(1-2c)i\lambda\hbar
t-\frac{1}{4}(\lambda\hbar t)^{2}\\\ 1-\frac{1}{2}(1-2c)i\lambda\hbar
t-\frac{1}{4}(\lambda\hbar t)^{2}&0\end{pmatrix}+\mathcal{O}(\lambda^{3}).$
(54)

If we calculate $S_{1xS}(t)$ exactly, we find that:

$S_{1xS}(t)=\frac{\hbar}{2}\begin{pmatrix}0&\frac{1}{2}[1+\cos(\lambda\hbar
t)+i(1-2c)\sin(\lambda\hbar t)]\\\ \frac{1}{2}[1+\cos(\lambda\hbar
t)-i(1-2c)\sin(\lambda\hbar t)]&0\end{pmatrix},$ (55)

which agrees with our perturbative calculations.

Now that we have calculated the one point operator $S_{1xS}$ up to second
order, we immediately have the two point operator up to first order. Note that
equation (43) can be written as:

$(S_{1x}(t_{1})S_{1x}(t_{2}))_{S}=S_{1xS}(t_{1})S_{1xS}(t_{2})+\mathcal{O}(\lambda^{2}).$
(56)

We can directly use the result (54) and obtain:

$\displaystyle(S_{1x}(t_{1})S_{1x}(t_{2}))_{S}=$
$\displaystyle\frac{\hbar}{2}\begin{pmatrix}1+\frac{1}{2}i(1-2c)(t_{1}-t_{2})\lambda&0\\\
0&1-\frac{1}{2}i(1-2c)(t_{1}-t_{2})\lambda\end{pmatrix}+\mathcal{O}(\lambda^{2}).$
(57)

We can immediately verify that this agrees with the exact two point operator.

##  VIII Conclusion

In this work, we have developed a technique to handle open quantum systems in
Heisenberg picture. We started by defining $N$-point operators and image
operators for system observables. We then perturbatively solved the reduced
Heisenberg equation, and figured out the $N$-point operators and the image
operators as a function of time. We then showed that all $N$-point operators
can be expressed solely in terms of 1-point operators in a way that is
dependent on the bath state. This expression can then be used to deform the
operator product on the Hilbert space of the open quantum system.

The Heisenberg picture described in this note is novel in the way it
associates multiple system operators to a single system observable. It is
natural to speculate that these multiplicity of operators play a crucial role
in say describing decoherence in terms of observables rather than states. It
is also clear that these image operators and their inter-relations should
encode the entanglement of the system with its environment. It would be nice
to make these statements precise. We also hope that the description in terms
of the Heisenberg picture throws new light on discussions about
interpretations of quantum mechanics. In particular, it would be interesting
to explore the relation, if any, between the description of the open system in
terms of multiple Heisenberg image operators and the Everett’s many-world
interpretation[8, 7].

As we mentioned in the introduction, a main motivation of this work is to
develop a systematic framework for open quantum field theories. The formalism
described here can readily be adopted to field theories, i.e., quantum systems
with infinite number of degrees of freedom. The Heisenberg evolution in QFTs
famously involve UV divergences which need to be regulated and counter-termed
in order to study renormalisation within these theories. A related question is
about operator products within open quantum field theories using the deformed
product introduced in this work. We hope to report on these issues in future.

## Acknowledgements

We would to thank Nilay Kundu and Mukund Rangamani for their valuable comments
regarding this work. RL would like to acknowledge his debt to the people of
India for their sustained and generous support to research in the basic
sciences.

## References

  * Ackerhalt et al. [1973] Jay R. Ackerhalt, Peter L. Knight, and Joseph H. Eberly.  Radiation reaction and radiative frequency shifts.  _Phys. Rev. Lett._ , 30:456–460, Mar 1973.  doi: 10.1103/PhysRevLett.30.456.  URL https://link.aps.org/doi/10.1103/PhysRevLett.30.456.
  * Alonso and de Vega [2007] Daniel Alonso and Inés de Vega.  Hierarchy of equations of multiple-time correlation functions.  _Phys. Rev. A_ , 75:052108, May 2007.  doi: 10.1103/PhysRevA.75.052108.  URL https://link.aps.org/doi/10.1103/PhysRevA.75.052108.
  * Breuer and Petruccione [2002] H.P. Breuer and F. Petruccione.  _The Theory of Open Quantum Systems_.  Oxford University Press, 2002.  ISBN 9780198520634.  URL https://books.google.co.in/books?id=0Yx5VzaMYm8C.
  * Breuer et al. [1998] H.P. Breuer, B. Kappler, and F. Petruccione.  Heisenberg picture operators in the stochastic wave function approach to open quantum systems.  _The European Physical Journal D - Atomic, Molecular, Optical and Plasma Physics_ , 1998.  doi: 10.1007/s100530050058.  URL https://doi.org/10.1007/s100530050058.
  * Clark et al. [2010] S R Clark, J Prior, M J Hartmann, D Jaksch, and M B Plenio.  Exact matrix product solutions in the heisenberg picture of an open quantum spin chain.  _New Journal of Physics_ , 12(2):025005, feb 2010\.  doi: 10.1088/1367-2630/12/2/025005.  URL https://doi.org/10.1088%2F1367-2630%2F12%2F2%2F025005.
  * de Vega and Alonso [2017] Inés de Vega and Daniel Alonso.  Dynamics of non-markovian open quantum systems.  _Rev. Mod. Phys._ , 89:015001, Jan 2017.  doi: 10.1103/RevModPhys.89.015001.  URL https://link.aps.org/doi/10.1103/RevModPhys.89.015001.
  * DeWitt and Graham [1973] Bryce DeWitt and Neill Graham, editors.  _The Many Worlds Interpretation of Quantum Mechanics_.  Priceton University Press, 1973.
  * Everett [1957] Hugh Everett.  Relative state formulation of quantum mechanics.  _Rev. Mod. Phys._ , 29:454–462, Jul 1957.  doi: 10.1103/RevModPhys.29.454.  URL https://link.aps.org/doi/10.1103/RevModPhys.29.454.
  * Feynman and Vernon [1963] R.P Feynman and F.L Vernon.  The theory of a general quantum system interacting with a linear dissipative system.  _Annals of Physics_ , 24:118 – 173, 1963.  ISSN 0003-4916.  doi: https://doi.org/10.1016/0003-4916(63)90068-X.  URL http://www.sciencedirect.com/science/article/pii/000349166390068X.
  * Gorini et al. [1976] Vittorio Gorini, Andrzej Kossakowski, and E. C. G. Sudarshan.  Completely positive dynamical semigroups of n‐level systems.  _Journal of Mathematical Physics_ , 17(5):821–825, 1976.  doi: 10.1063/1.522979.  URL https://aip.scitation.org/doi/abs/10.1063/1.522979.
  * Kheirandish [2019] Fardin Kheirandish.  Open quantum systems in heisenberg picture, 2019.
  * Kossakowski [1972] A. Kossakowski.  On quantum statistical mechanics of non-hamiltonian systems.  _Reports on Mathematical Physics_ , 3(4):247 – 274, 1972.  ISSN 0034-4877.  doi: https://doi.org/10.1016/0034-4877(72)90010-9.  URL http://www.sciencedirect.com/science/article/pii/0034487772900109.
  * Lindblad [1976] G. Lindblad.  On the generators of quantum dynamical semigroups.  _Communications in Mathematical Physics_ , 1976.  doi: 10.1007/BF01608499.  URL https://doi.org/10.1007/BF01608499.
  * Nakajima [1958] Sadao Nakajima.  On Quantum Theory of Transport Phenomena: Steady Diffusion.  _Progress of Theoretical Physics_ , 20(6):948–959, 12 1958.  ISSN 0033-068X.  doi: 10.1143/PTP.20.948.  URL https://doi.org/10.1143/PTP.20.948.
  * Pan et al. [2014] Yu Pan, Hadis Amini, Zibo Miao, John Gough, Valery Ugrinovskii, and Matthew R. James.  Heisenberg picture approach to the stability of quantum markov systems.  _Journal of Mathematical Physics_ , 55(6):062701, 2014.  doi: 10.1063/1.4884300.  URL https://doi.org/10.1063/1.4884300.
  * Schlosshauer [2007] Maximilian A. Schlosshauer.  _Decoherence and the Quantum-To-Classical Transition_.  Springer-Verlag Berlin Heidelberg, 2007.  doi: 10.1007/978-3-540-35775-9.
  * Yang et al. [2013] Li-Ping Yang, C. Y. Cai, D. Z. Xu, Wei-Min Zhang, and C. P. Sun.  Master equation and dispersive probing of a non-markovian process.  _Phys. Rev. A_ , 87:012110, Jan 2013.  doi: 10.1103/PhysRevA.87.012110.  URL https://link.aps.org/doi/10.1103/PhysRevA.87.012110.
  * Zoller and Gardiner [1997] Peter Zoller and C. W. Gardiner.  Quantum noise in quantum optics: the stochastic schrödinger equation, 1997.
  * Zwanzig [1960] Robert Zwanzig.  Ensemble method in the theory of irreversibility.  _The Journal of Chemical Physics_ , 33(5):1338–1341, 1960.  doi: 10.1063/1.1731409.  URL https://doi.org/10.1063/1.1731409.


# Analytic forms for the $e^{+}e^{-}$ annihilation cross sections around a
resonance including initial state radiation

Baoxin Liu  Zhenyu Zhang<EMAIL_ADDRESS>Xiang Zhou
<EMAIL_ADDRESS>Hubei Nuclear Solid Physics Key Laboratory, School of
Physics and Technology,
Wuhan University, Wuhan, Hubei 430072, P. R. China

###### Abstract

The exact analytic form of cross sections including initial state radiation
with Kuraev-Fadin radiative function are obtained for $e^{+}e^{-}$
annihilation around a resonance. Despite accounting for vacuum polarization
and center-of-mass energy spread effects, the precision remains below 0.1%,
meeting the accuracy requirements of quantum electrodynamics corrections up to
$\mathcal{O}(\alpha^{2})$. The analytic forms lead to an enhancement in the
precision of experimental measurements of physical parameters, such as
branching fractions, and demonstrate significantly improved computational
efficiency in the regression procedures.

Keywords

Initial state radiation, cross section, branching fraction, analytic form.

##  I Introduction

The $e^{+}e^{-}$ annihilation experiments provide cleaner experimental
environments for quarkonium decays than $p\bar{p}$ and hadron decays
experiments. The measured quarkonium decay widths, or the corresponding
branching fractions, can serve as inputs in phenomenological models or be used
for comparison with theoretical predictions to test our understanding of
quantum chromodynamics (QCD) [1]. There is an unavoidable background from the
continuum process, which directly produces the final state via $e^{+}e^{-}$
annihilation, i.e., $e^{+}e^{-}\to\gamma^{*}\to f$ [2]. It has been shown that
the cross section from the interference term will lead to imprecise branching
fractions for the broad resonance above open heavy flavor threshold, such as
$\psi(3770)$ [3].

A novel study shows that the ratio of the cross section from the interference
term with respect to the resonance are surprisingly large compared to the
precision of the current experiments even for the narrow resonances below the
open heavy flavor threshold, such as $J/\psi$ and $\psi(3686)$ [4]. Therefore,
to achieve precision measurements of the quarkonium decay widths or the
corresponding branching fractions, the cross sections around the resonance at
least three energies must be measured [4].

Initial state radiation (ISR) is an essential quantum electrodynamics (QED)
correction for the precise measurement of cross sections in $e^{+}e^{-}$
annihilations [5]. In the ISR process, the electron or positron emits one or
more photons before annihilation, thereby reducing the center-of-mass (c.m.)
energy. While the Born cross sections of the final state are of interest, it
is their corresponding ISR-corrected cross sections that are measured
experimentally [6]. The ISR-corrected cross section can be obtained by the
structure function (SF) method, which is an integral transformation of the
kernel cross section by the radiative function [7]. The SF method, with
corrections up to $\mathcal{O}(\alpha^{2})$, can achieve the required accuracy
of about $0.1\%$ for c.m. energy ranging from 0.2 to 10 GeV [7].

The parameters of the Born cross section can be obtained through the integral
method and the inverse transformation method. The integral method involves
fitting the experimental observed cross sections with the integral of the
radiative function multiplied by the model function of Born cross sections
[7]. On the other hand, the inverse transformation method utilizes the Born
cross section model function to fit the transformed cross sections from the
experimental observed ones through an iterative procedure [6, 8] or by solving
numerical integral equations [9]. The integral method is straightforward but
time-consuming, as it requires numerical integration in the regression
iterations [10]. Conversely, the inverse transformation method is time-
efficient but introduces an additional uncertainty from the transformation
process [6].

For the narrow resonance processes, such as $J/\psi$ or $\psi(2S)$, the widths
of the resonances are significantly less than the c.m. energy spread and then
the experimental observed cross sections become strongly correlated.
Therefore, the parameters of the Born cross section can only be obtained
through the integral method [11]. Moreover, the consuming time of the
regression process becomes much more burdensome since an additional
convolution with the c.m. energy spread must be included in the numerical
integration [10].

The integral method will be significantly accelerated if the analytic form of
the ISR-corrected cross sections is obtained [10]. The radiative function
proposed by Kuraev and Fadin (KF) consists of an exponentiated part and a
finite-order leading-logarithmic (LL) part, accounting for soft multi-photon
emission and hard collinear bremsstrahlung, respectively [7]. In 1987, R. N.
Cahn was the first to derive the analytic form with the exponentiated part in
the KF radiative function [12]. The approximate analytic form, including the
LL parts in the KF radiative function and the upper limit correction by the
exponential expansions, has been developed in subsequent works [13, 14, 11,
10]. However, it is shown below that for the ISR-corrected cross sections
around $J/\psi$, the accuracy of the approximated analytic form is within a
few percent and propagates in the same order to the estimated value by
regression for the hadronic branching fraction. In a sense, the approximated
analytic form introduces an imperceptible yet non-negligible systematic
uncertainty, which is comparable to the statistical or even the total
systematic uncertainty for experiments.

In this paper, we first provide the formalisms of the ISR corrected Born cross
sections and validate the precisions in Section II where the formulas of the
Born cross sections around a resonance and the corresponding exact analytic
forms of the ISR-corrected cross sections incorporating the KF radiative
function, referred as the KF analytic forms, are provided in Appendices A and
B, respectively. Section III and Section IV introduce the formalisms that
incorporate the vacuum polarization effect and the c.m. energy spread effect,
respectively, demonstrating the corresponding precisions. The fit tests of the
toy Monte Carlo (MC) samples in the vicinity of $J/\psi$ are presented in
Section V. Finally, Section VI gives discussions and conclusions.

##  II ISR corrected Born cross sections

In the vicinity of a resonance, the amplitude $A_{\mathrm{tot.}}^{f}$ of the
final state $f$ in $e^{+}e^{-}$ colliders is the coherent sum of both
resonance $A_{R}^{f}$ and continuum amplitudes $A_{C}^{f}$. The Born cross
section can be written as [4]

$\displaystyle\sigma_{\mathrm{Born}}^{f}(W)$ $\displaystyle\propto$
$\displaystyle\left|A_{\mathrm{tot.}}^{f}(W)\right|^{2}=\left|A_{C}^{f}(W)\right.$
(1) $\displaystyle\left.+A_{R}^{f}(W)e^{i\phi}\right|^{2},$

where $\phi$ is the relative phase between the continuum amplitude $A_{c}^{f}$
and the resonance amplitude $A_{R}^{f}$. Therefore, the Born cross section is
a sum of three parts as

$\sigma_{\mathrm{Born}}^{f}(W)=\sigma_{B_{C}}^{f}(W)+\sigma_{B_{R}}^{f}(W)+\sigma_{B_{I}}^{f}(W),$
(2)

where $\sigma_{B_{C}}^{f}(W)\propto\left|A_{C}^{f}(W)\right|^{2}$,
$\sigma_{B_{R}}^{f}(W)\propto\left|A_{R}^{f}(W)\right|^{2}$ and
$\sigma^{f}_{B_{I}}\propto 2\Re\\{A_{c}^{f}(W)A_{R}^{f}(W)\\}$ denote the Born
cross sections from continuum, resonance and interference contributions,
respectively.

When the kernel cross section is the Born cross section, the ISR corrected
cross section $\sigma_{\mathrm{ISR}}^{f}$ is an integral of the Born cross
section $\sigma_{\mathrm{Born}}^{f}$ times the radiation function [4], that is

$\displaystyle\sigma_{\mathrm{ISR}}^{f}(W)$ $\displaystyle=$
$\displaystyle\int_{0}^{1-\left(\frac{W_{\mathrm{min}}}{W}\right)^{2}}\mathrm{d}xF(x,W)$
(3) $\displaystyle\times\sigma_{\mathrm{Born}}^{f}(W\sqrt{1-x}),$

where $W$ is the c.m. energy of $e^{+}e^{-}$ annihilation and $W_{\text{min}}$
is the threshold energy equal to the invariant mass of the final states or the
experimental cut off energy. The KF radiative function $F(x,W)$ has the form
[7]

$\displaystyle F(x,W)$ $\displaystyle=$ $\displaystyle
x^{\beta-1}(1+\delta)-\beta\left(1-\frac{x}{2}\right)+\frac{1}{8}\beta^{2}\left[4(2-x)\right.$
(4)
$\displaystyle\left.\times\log\frac{1}{x}-\frac{1+3(1-x)^{2}}{x}\log(1-x)\right.$
$\displaystyle\left.-6+x\right],$

with
$\delta=\frac{3}{4}\beta+\frac{\alpha}{\pi}\left(\frac{\pi^{2}}{3}-\frac{\pi^{2}}{12}\right)$
and $\beta=\frac{2\alpha}{\pi}\left(2\log\frac{W}{m_{e}}-1\right)$.

The ISR corrected cross section integral $\sigma_{\mathrm{ISR}}^{f}$ is a sum
of three parts, which is

$\sigma_{\mathrm{ISR}}^{f}(W)=\sigma_{C}^{f}(W)+\sigma_{R}^{f}(W)+\sigma_{I}^{f}(W),$
(5)

where

$\sigma_{C}^{f}(W)\equiv\int_{0}^{1-\left(\frac{W_{\mathrm{min}}}{W}\right)^{2}}\mathrm{d}xF(x,W)\sigma_{B_{C}}^{f}(W\sqrt{1-x}),$
(6)
$\sigma_{R}^{f}(W)\equiv\int_{0}^{1-\left(\frac{W_{\mathrm{min}}}{W}\right)^{2}}\mathrm{d}xF(x,W)\sigma_{B_{R}}^{f}(W\sqrt{1-x}),$
(7)

and

$\sigma_{I}^{f}(W)\equiv\int_{0}^{1-\left(\frac{W_{\mathrm{min}}}{W}\right)^{2}}\mathrm{d}xF(x,W)\sigma_{B_{I}}^{f}(W\sqrt{1-x}),$
(8)

are the integrals of the continuum, resonance and interference contributions,
respectively.

For the process $e^{+}e^{-}\to\mu^{+}\mu^{-}$ in the vicinity of $J/\psi$, the
Born cross section is [11]

$\displaystyle\sigma_{\mathrm{Born}}^{\mu^{+}\mu^{-}}(W)$ $\displaystyle=$
$\displaystyle\frac{4\pi\alpha^{2}}{3W^{2}}\left|1+\frac{W^{2}}{M}\right.$ (9)
$\displaystyle\left.\times\frac{3\sqrt{\Gamma_{ee}\Gamma_{\mu\mu}}}{\alpha\left(W^{2}-M^{2}+iM\Gamma\right)}e^{i\phi}\right|^{2},$

where $\alpha$ is the fine structure constant, $M$ and $\Gamma$ are the mass
and total decay width of $J/\psi$, $\Gamma_{ee}$ and $\Gamma_{\mu\mu}$ are the
decay widths of $J/\psi\to e^{+}e^{-}$ and $\mu^{+}\mu^{-}$, respectively. For
the processes $e^{+}e^{-}\to 2(\pi^{+}\pi^{-})\pi^{0}$ and
$\eta\pi^{+}\pi^{-}$ with $\eta\to\pi^{+}\pi^{-}\pi^{0}$ in the vicinity of
$J/\psi$, the Born cross sections can be written in a general form as [11]

$\displaystyle\sigma_{\mathrm{Born}}^{5\pi}(W)$ $\displaystyle=$
$\displaystyle\left(\frac{\mathcal{A}}{W^{2}}\right)^{2}\frac{4\pi\alpha^{2}}{3W^{2}}\left|1+\right.$
(10)
$\displaystyle\left.\frac{3W^{2}\sqrt{\Gamma_{ee}\Gamma_{\mu\mu}}\mathcal{C}_{1}e^{i\phi}\left(1+\mathcal{C}_{2}e^{i\Phi}\right)}{\alpha
M\left(W^{2}-M^{2}+iM\Gamma\right)}\right|^{2},$

where $\frac{\mathcal{A}}{W^{2}}$ is the form factor, $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$ are the ratios of amplitudes, $\Phi$ is the phase between
the strong and electromagnetic decays from $J/\psi$. For the process of
$e^{+}e^{-}\to 2(\pi^{+}\pi^{-})\pi^{0}$, $\mathcal{C}_{1}$ and $\phi$ are
assumed to be 1 and 0, respectively [11].

Therefore, it can be derived for the cross sections from continuum, resonance
and interference contributions. For example, the Born cross section from the
continuum contribution of $e^{+}e^{-}\to\mu^{+}\mu^{-}$ is

$\sigma_{B_{C}}^{\mu^{+}\mu^{-}}(W)=\frac{4\pi\alpha^{2}}{3W^{2}}.\\\ $ (11)

The Born cross sections are composed by rational functions, while the KF
raditive function contains not only exponential functions but also logarithmic
functions, such as $\log x$, which make the ISR corrected cross section
integral has an integrable singularity at the lower limit $x=0$ in Eq. (3).
Fortunately, the integral still has an analytic form. For example, the
analytic form for the improper integral of $\log\frac{1}{x}$ times
$\sigma_{B_{C}}^{\mu^{+}\mu^{-}}(W)$ is

$\displaystyle\mathcal{I}(W)$ $\displaystyle=$
$\displaystyle\int_{0}^{1-\left(\frac{W_{\mathrm{min}}}{W}\right)^{2}}\mathrm{d}x\beta^{2}\log\frac{1}{x}\sigma_{B_{C}}^{\mu^{+}\mu^{-}}(W\sqrt{1-x})$
(12) $\displaystyle=$
$\displaystyle\int_{0}^{1-\left(\frac{W_{\mathrm{min}}}{W}\right)^{2}}\mathrm{d}x\beta^{2}\log\frac{1}{x}\frac{4\pi\alpha^{2}}{3W^{2}(1-x)}$
$\displaystyle=$
$\displaystyle-\beta^{2}\frac{4\pi\alpha^{2}}{3W^{2}}\int_{0}^{1-\left(\frac{W_{\mathrm{min}}}{W}\right)^{2}}\mathrm{d}x\frac{\log
x}{1-x}$ $\displaystyle=$
$\displaystyle-\beta^{2}\frac{4\pi\alpha^{2}}{3W^{2}}\left[\mathrm{Li}_{2}\left(\frac{W_{\mathrm{min}}^{2}}{W^{2}}\right)-\frac{\pi^{2}}{6}\right],$

where $\mathrm{Li}_{2}(x)$ is the Spence’s function which is defined as

$\mathrm{Li}_{2}(z)=\int_{1}^{1-z}\frac{\log(x)}{1-x}\text{d}x.$ (13)

Since the analytic forms are tedious, we put them into the appendix B.

In Fig. 1, we compare the ISR corrected cross section
$\sigma^{\mu^{+}\mu^{-}}_{\mathrm{NI}}$ obtained through the numerical
integration (NI) with the KF analytic form
$\sigma^{\mu^{+}\mu^{-}}_{\mathrm{KF}}$, and the approximated analytic form
$\sigma^{\mu^{+}\mu^{-}}_{\mathrm{Approx}}$ for the
$e^{+}e^{-}\to\mu^{+}\mu^{-}$ process around the $J/\psi$ resonance. The
relative deviation between $\sigma^{\mu^{+}\mu^{-}}_{\mathrm{KF}}$ and
$\sigma^{\mu^{+}\mu^{-}}_{\mathrm{NI}}$ is less than $10^{-10}$, coinciding
with the numerical integration’s error tolerance, thereby affirming the
precision of our analytic form. Moreover, the relative deviation between
$\sigma_{\text{NI}}^{\mu^{+}\mu^{-}}$ and
$\sigma_{\text{Approx}}^{\mu^{+}\mu^{-}}$ reaches approximately 0.3% at the
incoherent phase $\phi=90^{\circ}$ where $\sigma_{B_{I}}^{\mu^{+}\mu^{-}}$ and
then $\sigma_{I}^{\mu^{+}\mu^{-}}$ vanish. Therefore, the precision is about
0.3% for the sum of the approximated analytic forms
$\sigma_{C}^{\mu^{+}\mu^{-}}$ and $\sigma_{R}^{\mu^{+}\mu^{-}}$ where the
precision of the approximated analytic form $\sigma_{R}^{\mu^{+}\mu^{-}}$ is
only 0.1% [13]. Conversely, at $\phi=0^{\circ}$ or $180^{\circ}$, the relative
deviation between $\sigma_{\text{NI}}^{\mu^{+}\mu^{-}}$ and
$\sigma_{\text{Approx}}^{\mu^{+}\mu^{-}}$ is about 3%. It is evident that the
deviation comes from $\sigma_{I}^{\mu^{+}\mu^{-}}$ is 3% in the approximated
analytic form.

Figure 1:  The comparisons of ISR corrected cross sections for
$e^{+}e^{-}\to\mu^{+}\mu^{-}$ at $\phi=0^{\circ}$ at $\phi=0^{\circ}$, 90∘ and
180∘. The Born cross section parameters are set to be PDG world averages [15]
and $W_{\mathrm{min}}=2.0$ GeV. In the top row, black lines represent KF
method results, blue dashed lines represent Approx method results and red
dotted lines represent NI method results. In the middle row, red dashed lines
represent the results of
$\left|\frac{\sigma^{\mu^{+}\mu^{-}}_{\text{Approx}}}{\sigma^{\mu^{+}\mu^{-}}_{\text{NI}}}-1\right|$.
In the bottom row, black lines represent the results of
$\left|\frac{\sigma^{\mu^{+}\mu^{-}}_{\text{KF}}}{\sigma^{\mu^{+}\mu^{-}}_{\text{NI}}}-1\right|$.

##  III Vacuum polarization effect

When the vacuum polarization (VP) effect is considered, the Born cross section
in the vicinity of a resonance is corrected as [16, 14]

$\displaystyle\sigma_{\mathrm{Born-VP}}^{f}(W)$ $\displaystyle=$
$\displaystyle\sigma_{B_{C}^{\mathrm{VP}}}^{f}(W)+\sigma_{R_{B}^{\mathrm{VP}}}^{f}(W)+\sigma_{B_{I}^{\mathrm{VP}}}^{f}(W)$
(14) $\displaystyle=$
$\displaystyle\frac{\sigma_{B_{C}}^{f}(W)}{\left|1-\Pi_{0}(W)\right|^{2}}+\sigma_{B_{R}}^{f}(W)$
$\displaystyle+\frac{\sigma_{B_{I}}^{f}(W)}{\left|1-\Pi_{0}(W)\right|},$

where $\Pi_{0}(W)$ is the non-resonant VP factor [16].

The ISR corrected cross section with VP effect is also a sum of three parts
which is

$\displaystyle\sigma_{\mathrm{ISR-VP}}^{f}(W)$ $\displaystyle=$
$\displaystyle\sigma_{C-\mathrm{VP}}^{f}(W)+\sigma_{R-\mathrm{VP}}^{f}(W)$
(15) $\displaystyle+\sigma_{I-\mathrm{VP}}^{f}(W),$

where

$\displaystyle\sigma_{R-\mathrm{VP}}^{f}(W)\equiv\sigma_{R}^{f}(W)$
$\displaystyle=$
$\displaystyle\int_{0}^{1-\left(\frac{W_{\mathrm{min}}}{W}\right)^{2}}\mathrm{d}xF(x,W)$
(16) $\displaystyle\times\sigma_{B_{R}}^{f}(W\sqrt{1-x}),$
$\displaystyle\sigma_{C-\mathrm{VP}}^{f}(W)$ $\displaystyle\equiv$
$\displaystyle\int_{0}^{1-\left(\frac{W_{\mathrm{min}}}{W}\right)^{2}}\mathrm{d}xF(x,W)$
(17)
$\displaystyle\times\frac{\sigma_{B_{C}}^{f}(W\sqrt{1-x})}{\left|1-\Pi_{0}(W\sqrt{1-x})\right|^{2}},$
$\displaystyle\simeq$
$\displaystyle\frac{\sigma_{C}^{f}(W)}{\left|1-\Pi_{0}(W)\right|^{2}},$ (18)

and

$\displaystyle\sigma_{I-\mathrm{VP}}^{f}(W)$ $\displaystyle\equiv$
$\displaystyle\int_{0}^{1-\left(\frac{W_{\mathrm{min}}}{W}\right)^{2}}\mathrm{d}xF(x,W)$
(19)
$\displaystyle\times\frac{\sigma_{B_{I}}^{f}(W\sqrt{1-x})}{\left|1-\Pi_{0}(W\sqrt{1-x})\right|},$
$\displaystyle\simeq$
$\displaystyle\frac{\sigma_{I}^{f}(W)}{\left|1-\Pi_{0}(W)\right|},$ (20)

are integrals of the resonance, continuum, and interference contributions,
respectively. Approximations are used in Eqs. (18) and (20) because the non-
resonance VP factors $\Pi_{0}(W)$ are numerically obtained [17].

In Fig. 2, we compare the ISR corrected cross section including VP effect
$\sigma_{\mathrm{NI-VP}}^{\mu^{+}\mu^{-}}$ obtained through NI with
$\sigma_{\mathrm{KF-VP}}^{\mu^{+}\mu^{-}}$ obtained using the KF analytic form
and $\sigma_{\mathrm{Approx-VP}}^{\mu^{+}\mu^{-}}$ obtained using the
approximated analytic form. The relative deviation between
$\sigma_{\mathrm{NI-VP}}^{\mu^{+}\mu^{-}}$ and $\sigma_{\mathrm{Approx-
VP}}^{\mu^{+}\mu^{-}}$ is at the order of $10^{-4}$, primarily arising from
the simplifications made in Eq. (18) and (20). The deviation meets the 0.1%
precision requirement of the SF method with corrections up to
$\mathcal{O}(\alpha^{2})$ [7], thereby fulfilling the criteria for practical
applicability. Furthermore, the relative deviation between
$\sigma_{\mathrm{NI-VP}}^{\mu^{+}\mu^{-}}$ and $\sigma_{\mathrm{Approx-
VP}}^{\mu^{+}\mu^{-}}$ remains about 0.3% at $\phi=90^{\circ}$, and 3% at
$\phi=0^{\circ}$ and $180^{\circ}$. It is evident that the precision of
$\sigma_{\mathrm{Approx-VP}}^{\mu^{+}\mu^{-}}$ directly inherits from the one
of $\sigma_{\mathrm{Approx}}^{\mu^{+}\mu^{-}}$.

Figure 2:  The comparisons of ISR corrected cross sections with VP effect for
$e^{+}e^{-}\to\mu^{+}\mu^{-}$ at $\phi=0^{\circ}$, $90^{\circ}$ and
$180^{\circ}$. The Born cross section parameters are set to be PDG world
averages [15] and $W_{\mathrm{min}}=2.0$ GeV. In the top row, black lines
represent KF method results, blue dashed lines represent Approx method results
and red dotted lines represent NI method results. In the bottom row, red
dashed lines represent the results of
$\left|\frac{\sigma^{\mu^{+}\mu^{-}}_{\text{Approx-
VP}}}{\sigma^{\mu^{+}\mu^{-}}_{\text{NI-VP}}}-1\right|$ and black lines
represent the results of $\left|\frac{\sigma^{\mu^{+}\mu^{-}}_{\text{KF-
VP}}}{\sigma^{\mu^{+}\mu^{-}}_{\text{NI-VP}}}-1\right|$.

##  IV c.m. energy spread effect

For the narrow resonances, such as $J/\psi$ whose decay width is 92.6 keV, the
c.m. energy spreads of $e^{+}e^{-}$ are much larger than the resonance widths.
For example, the c.m. energy spread around $J/\psi$ in BEPCII is less than 1
MeV [11]. Therefore, the effect for the c.m. energy spread must be considered
in the experimentally observed cross section $\sigma_{\mathrm{exp}}^{f}(W)$ by
a Gaussian convolution with $\sigma_{\mathrm{ISR-VP}}^{f}(W)$ which is

$\displaystyle\sigma_{\mathrm{ISR-exp}}^{f}(W)$ $\displaystyle=$
$\displaystyle\int_{W-n\S_{E}}^{W+nS_{E}}\frac{1}{\sqrt{2\pi}S_{E}}\exp\left(\frac{-\left(W-W^{\prime}\right)^{2}}{2S_{E}^{2}}\right)$
(21) $\displaystyle\times\sigma_{\mathrm{ISR-
VP}}^{f}\left(W^{\prime}\right)dW^{\prime}.$

In Fig. 3, we compare ISR corrected cross sections including both the VP and
c.m. energy spread effects $\sigma^{\mu^{+}\mu^{-}}_{\mathrm{NI-exp}}$
obtained through the NI with $\sigma^{\mu^{+}\mu^{-}}_{\mathrm{KF-exp}}$
obtained using the KF analytic form and
$\sigma^{\mu^{+}\mu^{-}}_{\mathrm{Approx-exp}}$ obtained using the
approximated analytic form. The relative deviation between
$\sigma^{\mu^{+}\mu^{-}}_{\mathrm{NI-exp}}$ and
$\sigma^{\mu^{+}\mu^{-}}_{\mathrm{KF-exp}}$ remains at the order of $10^{-4}$.
Moreover, the relative deviation between $\sigma_{\mathrm{ISR-
VP}}^{\mathrm{NI}}$ and $\sigma_{\mathrm{ISR-VP}}^{\mathrm{PE}}$ is still 0.3%
at $\Phi=90^{\circ}$ and 3% at $\Phi=0^{\circ}$ and $\Phi=180^{\circ}$. It is
evident that the Gaussian convolution of the c.m. energy spread does not
impact on the precision.

Figure 3: The comparisons of ISR corrected cross section with VP effect and
c.m. energy spread effect for $e^{+}e^{-}\to\mu^{+}\mu^{-}$ at
$\phi=0^{\circ}$ at $\phi=0^{\circ}$, $90^{\circ}$ and $180^{\circ}$. The Born
cross section parameters are set to be PDG world averages [15] and
$W_{\mathrm{min}}=2.0$ GeV. In the top row, black lines represent KF method
results, blue dashed lines represent Approx method results and red dotted
lines represent NI method results. In the bottom row, red dashed lines
represent the results of $\left|\frac{\sigma^{\mu^{+}\mu^{-}}_{\text{Approx-
exp}}}{\sigma^{\mu^{+}\mu^{-}}_{\text{NI-exp}}}-1\right|$ and black lines
represent the results of $\left|\frac{\sigma^{\mu^{+}\mu^{-}}_{\text{KF-
exp}}}{\sigma^{\mu^{+}\mu^{-}}_{\text{NI-exp}}}-1\right|$.

In Table 1, we present a comparison of the consuming time for
$\sigma^{f}_{\mathrm{KF-exp}}$, $\sigma^{f}_{\mathrm{Approx-exp}}$, and
$\sigma^{f}_{\mathrm{NI-exp}}$ for the processes $e^{+}e^{-}\to\mu^{+}\mu^{-}$
and $2(\pi^{+}\pi^{-})\pi^{0}$. Both the KF and approximated analytic forms
significantly reduce the consuming time. Notably, $\sigma^{f}_{\mathrm{KF-
exp}}$ yields the shortest computing time for the hadronic decay process.

Table 1:  The comparison of the consuming times by parallel computation on 12 CPU cores within Python-based codes.  | $e^{+}e^{-}\to\mu^{+}\mu^{-}$ | $e^{+}e^{-}\to 5\pi$
---|---|---
$N_{\text{points}}$ | $T_{\text{KF-exp}}\mathrm{(s)}$ | $T_{\text{Approx-exp}}(\mathrm{(s)}$ | $T_{\text{NI-exp}}\mathrm{(s)}$ | $T_{\text{KF-exp}}\mathrm{(s)}$ | $T_{\text{Approx-exp}}\mathrm{(s)}$ | $T_{\text{NI-exp}}\mathrm{(s)}$
10 | 0.18 | 0.13 | 9.90 | 0.15 | 0.59 | 15.24
100 | 0.48 | 0.34 | 32.59 | 0.34 | 2.06 | 50.73
1000 | 2.03 | 1.35 | 307.99 | 2.44 | 17.51 | 469.61

##  V Fit results using toy MC samples

We use the toy MC samples for $e^{+}e^{-}\to\mu^{+}\mu^{-}$ and $e^{+}e^{-}\to
2(\pi^{+}\pi^{-})\pi^{0}$ around $J/\psi$ to test our analytic form and
compare it to the approximated one. Firstly, the toy MC samples for the ISR
corrected cross sections of $e^{+}e^{-}\to\mu^{+}\mu^{-}$ are generated with
1% uncertainty at 20 c.m. energies around $J/\psi$. The minimized
$\chi^{2}_{\mu^{+}\mu^{-}}$ is performed with the free parameters $M$, $S_{E}$
and $\phi$. The $\chi^{2}_{\mu^{+}\mu^{-}}$ is Eq. (22)

$\chi^{2}_{\mu^{+}\mu^{-}}(M,S_{E},\phi)=\sum_{i=1}^{20}\left[\frac{\sigma^{\mu^{+}\mu^{-}}_{i}-\sigma_{\text{ISR-
exp}}^{\mu^{+}\mu^{-}}(W_{i};M,S_{E},\phi)}{\Delta\sigma^{\text{MC}}_{i}}\right]^{2},$
(22)

where $\sigma^{\mu^{+}\mu^{-}}_{i}$ is the cross section at every energy point
$W_{i}$, $\Delta\sigma^{\mu^{+}\mu^{-}}_{i}$ is the corresponding uncertainty.
The estimators $\widehat{M}$ and $\widehat{S_{E}}$ and their corresponding
uncertainties $\sigma_{\widehat{M}}$ and $\sigma_{\widehat{S_{E}}}$ are taken
as the information of the $e^{+}e^{-}$ collider.

Secondly, the toy MC samples for the ISR corrected cross sections of
$e^{+}e^{-}\to 2(\pi^{+}\pi^{-})\pi^{0}$ are then generated with 1%
uncertainty at the same 20 c.m. energies around $J/\psi$. The minimized
$\chi^{2}_{5\pi}$ is performed with the free parameters $\mathcal{A}$,
$\mathcal{C}_{2}$ and $\Phi$, and the nuisance parameters $\widehat{M}$ and
$\widehat{S_{E}}$. Therefore, the estimators $\widehat{\mathcal{A}}$,
${\widehat{\mathcal{C}}}_{2}$ and $\widehat{\Phi}$ can be obtained. The
$\chi^{2}_{5\pi}$ is Eq. (23)

$\displaystyle\chi^{2}_{5\pi}(\mathcal{A},\mathcal{C}_{2},\Phi)$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{20}\left[\frac{\sigma^{5\pi}_{i}-\sigma_{\text{ISR-
exp}}^{5\pi}(W_{i};\mathcal{A},\mathcal{C}_{2},\Phi,M,S_{E})}{\Delta\sigma^{5\pi}_{i}}\right]^{2}$
(23)
$\displaystyle+\left(\frac{M-\widehat{M}}{\sigma_{\widehat{M}}}\right)^{2}+\left(\frac{S_{E}-\widehat{S_{E}}}{\sigma_{\widehat{S_{E}}}}\right)^{2}.$

The branching fraction for $J/\psi\to 2(\pi^{+}\pi^{-})\pi^{0}$ satisfies [11]

$\mathrm{Br}=\left(\frac{\mathcal{A}}{W^{2}}\right)^{2}\Gamma_{\mu\mu}\left|(1+\mathcal{C}_{2}e^{i\Phi})\right|^{2}.$
(24)

Finally, the above toy MC sample generations and $\chi^{2}$ minimizations are
repeated 10000 times. The fit results of $e^{+}e^{-}\to\mu^{+}\mu^{-}$ are
shown in Fig. 4. The estimators of $\widehat{M}$ and $\widehat{S_{E}}$ are
centered around the MC truth values for both KF and approximated analytic
forms. Only the estimator $\widehat{\phi}$ of the KF analytic form is centered
around the MC truth value. While the distribution of the estimator
$\widehat{\phi}$ of the approximated analytic form deviate more than 2$\sigma$
from the MC truth value. The $\chi^{2}_{\mu^{+}\mu^{-}}$ of the KF analytic
form is better than that of the approximated analytic form.

Figure 4:  The histograms of the estimators obtained from toy MC samples for
$e^{+}e^{-}\to\mu^{+}\mu^{-}$. (a), (b), (c) and (d) show the distributions of
$\Delta{\widehat{M}}=\widehat{M}-M^{\mathrm{MCTruth}}_{J/\psi}$,
$\widehat{S_{E}}$, $\widehat{\phi}$ and $\chi^{2}_{\mu^{+}\mu^{-}}$,
respectively. Histograms with black edges represent the KF analytic form, blue
edges represent Approx analytic form and red lines are the values of MC truth.

The fit results for the process $e^{+}e^{-}\to 2(\pi^{+}\pi^{-})\pi^{0}$ are
depicted in Fig. 5. Two solutions, labeled as Solution I and II, are obtained.
The distributions of $\chi^{2}_{5\pi}$ are identical for both Solution I and
II. The $\chi^{2}_{5\pi}$ of the KF analytic form is better than that of the
approximated analytic form. Only in Solution I the estimator $\widehat{\Phi}$
is centered around the MC truth value of 90∘. While in Solution II,
$\widehat{\Phi}$ is centered around $-90^{\circ}$. Furthermore, in Solution I,
only the estimator $\widehat{\mathrm{Br}}$ of the KF analytic form is centered
around the MC truth value, whereas the estimator of the approximated analytic
form deviates significantly. The toy MC test shows that the KF analytic form
avoids a subtle yet non-negligible uncertainty that arises from employing the
approximated analytic form.

Figure 5:  The histograms of the estimators obtained from toy MC samples for
$e^{+}e^{-}\to 2(\pi^{+}\pi^{-})\pi^{0}$. (a), (b), (c) and (g) show the
distributions of $\widehat{\mathcal{A}}$, $\widehat{\mathcal{C}}_{2}$, $\Phi$
and $\widehat{\mathrm{Br}}$ for Solution I, respectively; (d), (e), (f) and
(h) show the distributions of $\mathcal{A}$, $\mathcal{C}_{2}$, $\Phi$ and
$\widehat{\mathrm{Br}}$ for Solution II, respectively; (i) shows
$\chi^{2}_{5\pi}$ distributions for both Solution I and II. Histograms with
black edges represent the KF analytic form, blue edges represent Approx
analytic form and red lines are the values of MC truth.

##  VI Discussions and Conclusions

The ISR correction is essential for precise cross section measurements in
$e^{+}e^{-}$ annihilation. The structure functions method provides an integral
transformation from the Born cross sections to the experimentally observed
ones. For the first time, we present the exact analytic form for the ISR
corrected cross section aorund a resonance with the KF radiative function.
Even including the vacuum polarization (VP) and c.m. energy spread effects,
the precision remains matches the accuracy requirement of the structure
functions (SF) method up to $\mathcal{O}(\alpha^{2})$ [7]. The analytic form
can accelerate the regression procedure used to extract physical parameters,
such as the branching fraction, by over 150 times. Utilizing the toy MC
samples, we reveal a non-negligible few percent systematic uncertainty caused
by the approximated analytic form.

For the currently running $e^{+}e^{-}$ collision experiments, the BESIII
experiment has been accumulated 10 billion $J/\psi$ events [18] and 3 billion
$\psi(3686)$ events [19] which are at least 6 times larger than those used in
previous BESIII measurements. Furthermore, the Belle II experiment plans to
take about 500 fb-1 data for each vector bottomnium state [20] which are at
least ten or huandreds of times larger than the Belle experiment. The KF
analytic form will be helpful to handle the ISR correction properly in
currently running experiments at BESIII and Belle II, as well as in planned
ones, such as the super-tau-charm factories (STCF) [21] and the super-$J/\psi$
factory [22].

###### Acknowledgements.

X. Z. would like to acknowledge useful conversations with Y. D. Wang and K.
Zhu. This work has been supported by the National Natural Science Foundation
of China (NSFC) under Contract Nos. U2032114 and 12192265).

##  Appendix A Analytic forms for Born Cross Sections

The Born cross section of the process $e^{+}e^{-}\to\mu^{+}\mu^{-}$ is [11]

$\displaystyle\sigma_{\mathrm{Born}}^{\mu^{+}\mu^{-}}(W)$ $\displaystyle=$
$\displaystyle\frac{4\pi\alpha^{2}}{3W^{2}}\left|1+\frac{W^{2}}{M}\right.$
(25)
$\displaystyle\left.\times\frac{3\sqrt{\Gamma_{ee}\Gamma_{\mu\mu}}}{\alpha\left(W^{2}-M^{2}+iM\Gamma\right)}e^{i\phi}\right|^{2},$

which can be taken apart into three parts

$\sigma^{\mu^{+}\mu^{-}}_{B_{c}}(W)=\frac{4\pi\alpha^{2}}{3W^{2}},$ (26)
$\sigma^{\mu^{+}\mu^{-}}_{B_{R}}(W)=\frac{12\pi
W^{2}\Gamma_{ee}\Gamma_{\mu\mu}}{M^{2}((W^{2}-M^{2})^{2}+M^{2}\Gamma^{2})},$
(27) $\displaystyle\sigma^{\mu^{+}\mu^{-}}_{B_{I}}(W)$ $\displaystyle=$
$\displaystyle\frac{8\pi\alpha\sqrt{\Gamma_{ee}\Gamma_{\mu\mu}}}{M((W^{2}-M^{2})^{2}+M^{2}\Gamma^{2})}$
(28)
$\displaystyle\times\left((W^{2}-M^{2})\text{cos}\phi+M\Gamma\text{sin}\phi\right).$

The Born cross section of the process $e^{+}e^{-}\to 2(\pi^{+}\pi^{-})\pi^{0}$
is [11]

$\displaystyle\sigma_{\mathrm{Born}}^{5\pi}(W)$ $\displaystyle=$
$\displaystyle\left(\frac{\mathcal{A}}{W^{2}}\right)^{2}\frac{4\pi\alpha^{2}}{3W^{2}}\left|1+\frac{W^{2}}{M}\right.$
(29)
$\displaystyle\left.\times\frac{3\sqrt{\Gamma_{ee}\Gamma_{\mu\mu}}\mathcal{C}_{1}e^{i\phi}\left(1+\mathcal{C}_{2}e^{i\Phi}\right)}{\alpha\left(W^{2}-M^{2}+iM\Gamma\right)}\right|^{2},$

which can also be taken apart into three parts:

$\sigma_{B_{c}}^{5\pi}(W)=\left(\frac{\mathcal{A}}{W^{2}}\right)^{2}\frac{4\pi\alpha^{2}}{3W^{2}},$
(30) $\displaystyle\sigma_{B_{R}}^{5\pi}(W)$ $\displaystyle=$
$\displaystyle\mathcal{C}_{1}^{2}(1+\mathcal{C}_{2}^{2}+2\mathcal{C}_{2}\text{cos}\Phi)$
(31)
$\displaystyle\times\frac{12\pi\mathcal{A}^{2}\Gamma_{ee}\Gamma_{\mu\mu}}{M^{2}W^{2}((W^{2}-M^{2})^{2}+M^{2}\Gamma^{2})},$
$\displaystyle\sigma_{B_{I}}^{5\pi}(W)$ $\displaystyle=$
$\displaystyle\frac{8\pi\alpha\mathcal{A}^{2}\sqrt{\Gamma_{ee}\Gamma_{\mu\mu}}}{W^{4}M((W^{2}-M^{2})^{2}+M^{2}\Gamma^{2})}$
(32)
$\displaystyle\times\left[(W^{2}-M^{2})\mathcal{C}_{1}(\text{cos}\phi+\mathcal{C}_{2}\text{cos}\phi\text{cos}\Phi\right.$
$\displaystyle\left.-\mathcal{C}_{2}\text{sin}\phi\text{sin}\Phi)+M\Gamma\mathcal{C}_{1}(\text{sin}\phi+\right.$
$\displaystyle\left.\mathcal{C}_{2}\text{cos}\phi\text{sin}\Phi+\mathcal{C}_{2}\text{sin}\phi\text{cos}\Phi)\right].$

##  Appendix B Analytic forms of ISR corrected cross sections

The ISR corrected cross section from the continuum part of
$e^{+}e^{-}\to\mu^{+}\mu^{-}$ is

$\displaystyle\sigma_{C}^{\mu^{+}\mu^{-}}$ $\displaystyle=$
$\displaystyle\frac{4\pi\alpha^{2}}{3W^{2}}I_{0}.$ (33)

The ISR corrected cross section from the resonance part of
$e^{+}e^{-}\to\mu^{+}\mu^{-}$ is

$\displaystyle\sigma^{\mu^{+}\mu^{-}}_{R}$ $\displaystyle=$
$\displaystyle\frac{6\pi\Gamma_{ee}\Gamma_{\mu\mu}W^{2}}{M^{5}\Gamma}\left[\left(1+\frac{A}{B}\right)I_{1}-\frac{1}{B}I_{2}\right]$
(34) $\displaystyle+\mathrm{c.c.}.$

The ISR corrected cross section from the interference part of
$e^{+}e^{-}\to\mu^{+}\mu^{-}$ is

$\displaystyle\sigma^{\mu^{+}\mu^{-}}_{I}$ $\displaystyle=$
$\displaystyle\left[\frac{4\pi\alpha\sqrt{\Gamma_{ee}\Gamma_{\mu\mu}}}{M^{3}\Gamma}(\Gamma\text{sin}\phi-M\text{cos}\phi)+\right.$
(35)
$\displaystyle\left.\frac{4\pi\alpha\sqrt{\Gamma_{ee}\Gamma_{\mu\mu}}}{M^{4}\Gamma}W^{2}\text{cos}\phi\left(1+\frac{A}{B}\right)\right]I_{1}$
$\displaystyle~{}~{}-\frac{I_{2}W^{2}}{B}\frac{4\pi\alpha\sqrt{\Gamma_{ee}\Gamma_{\mu\mu}}}{M^{4}\Gamma}\text{cos}\phi+\mathrm{c.c.}.$

The ISR corrected cross section from the continuum part of $e^{+}e^{-}\to
2(\pi^{+}\pi^{-})\pi^{0}$ is

$\displaystyle\sigma^{5\pi}_{C}$ $\displaystyle=$
$\displaystyle\frac{4\pi\alpha^{2}\mathcal{A}^{2}}{3W^{6}}I_{3}.$ (36)

The ISR corrected cross section from the resonance part of $e^{+}e^{-}\to
2(\pi^{+}\pi^{-})\pi^{0}$ is

$\displaystyle\sigma^{5\pi}_{R}$ $\displaystyle=$
$\displaystyle\mathcal{C}_{1}^{2}(1+\mathcal{C}_{2}^{2}+2\mathcal{C}_{2}\text{cos}\Phi)\frac{6\pi\mathcal{A}^{2}\Gamma_{ee}\Gamma_{\mu\mu}}{\Gamma
M^{5}W^{2}}$ (37)
$\displaystyle\left(\frac{1}{A+B}I_{0}+\frac{B}{A+B}I_{1}\right)+\mathrm{c.c.}.$

The ISR corrected cross section from the interference part of $e^{+}e^{-}\to
2(\pi^{+}\pi^{-})\pi^{0}$ is

$\displaystyle\sigma^{5\pi}_{I}$ $\displaystyle=$
$\displaystyle\frac{4\pi\alpha\mathcal{A}^{2}\sqrt{\Gamma_{ee}\Gamma_{\mu\mu}}}{M}\left\\{\left[\frac{\mathcal{C}_{1}}{W^{4}M^{2}}\right.\right.$
(38)
$\displaystyle\left.\left.\times(\text{sin}\phi+\mathcal{C}_{2}\text{cos}\phi\text{sin}\Phi+\mathcal{C}_{2}\text{sin}\phi\text{cos}\Phi)\right.\right.$
$\displaystyle\left.\left.-\frac{\mathcal{C}_{1}(\text{cos}\phi+\mathcal{C}_{2}\text{cos}\phi\text{cos}\Phi-\mathcal{C}_{2}\text{sin}\phi\text{sin}\Phi)}{\Gamma
W^{4}M}\right]\right.$
$\displaystyle\left.\times\left[\frac{B}{(A+B)^{2}}I_{0}+\left(\frac{B}{A+B}\right)^{2}I_{1}+\frac{1}{A+B}I_{4}\right]+\right.$
$\displaystyle\left.\frac{\mathcal{C}_{1}(\text{cos}\phi+\mathcal{C}_{2}\text{cos}\phi\text{cos}\Phi-\mathcal{C}_{2}\text{sin}\phi\text{sin}\Phi)}{\Gamma
W^{2}M^{3}}\left(\frac{1}{A+B}I_{0}\right.\right.$
$\displaystyle\left.\left.+\frac{B}{A+B}I_{1}\right)\right\\}+\mathrm{c.c.},$

where $\mathrm{c.c.}$ represents the complex conjugate of the former part,
$A=\frac{\Gamma}{M}+i\left[\left(\frac{W}{M}\right)^{2}-1\right]$,
$B=-i\left(\frac{W}{M}\right)^{2}$ and
$b=1-\left(\frac{W_{\text{min}}}{W}\right)^{2}$. The five integrations
$I_{0}$, $I_{1}$, $I_{2}$, $I_{3}$ and $I_{4}$ are

$\displaystyle I_{0}$ $\displaystyle=$
$\displaystyle\beta(1+\delta)\mathcal{B}(b,\beta,0)-\frac{b\beta}{2}+\log(1-b)\left(\frac{1}{4}\beta^{2}+\frac{\beta}{2}+\frac{3}{8}\beta^{2}b\right)+\frac{1}{16}\beta^{2}\log^{2}(1-b)$
(39) $\displaystyle~{}~{}-\frac{\beta^{2}}{2}b\log
b-\frac{\beta^{2}}{2}\left(\text{Li}_{2}(1-b)-\frac{\pi^{2}}{6}-\text{Li}_{2}(b)\right).$

$\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\beta(1+\delta)\left(\frac{b^{\beta}\,{}_{2}F_{1}(1,\beta,1+\beta,-bB/A)}{A\beta}\right)+\left(\frac{\beta^{2}}{8}+\frac{\beta}{2}\right)\left(\frac{b}{B}-\frac{A}{B^{2}}\log\frac{A+Bb}{A}\right)-$
(40)
$\displaystyle~{}~{}\left(\beta+\frac{3}{4}\beta^{2}\right)\frac{1}{B}\log\frac{A+Bb}{A}+\frac{3}{4}\beta^{2}\left[-\frac{1}{B}\left[\text{Li}_{2}\left(\frac{A+Bb}{A+B}\right)-\text{Li}_{2}\left(\frac{A}{A+B}\right)\right.\right.$
$\displaystyle~{}~{}\left.\left.+\log\frac{B}{A+B}\log\frac{A+Bb}{A}\right]\right]-\frac{3}{8}\beta^{2}\frac{1}{B}\left[-b+(b-1)\log(1-b)+\right.$
$\displaystyle~{}~{}\left.\frac{A}{B}\left[\text{Li}_{2}\left(\frac{A+Bb}{A+B}\right)-\text{Li}_{2}\left(\frac{A}{A+B}\right)+\log\frac{B}{A+B}\log\frac{A+Bb}{A}\right]\right]$
$\displaystyle~{}~{}-\beta^{2}\frac{1}{B}\left[-\frac{\pi^{2}}{6}+\frac{1}{2}\log^{2}\left(1+\frac{B}{A}b\right)+\text{Li}_{2}\left(\frac{A}{A+Bb}\right)-\log\frac{B}{A}\log\left(1+\frac{B}{A}b\right)\right]$
$\displaystyle~{}~{}+\frac{1}{2}\beta^{2}\frac{1}{B}\left[b(\log
b-1)-\frac{A}{B}\left[-\frac{\pi^{2}}{6}+\frac{1}{2}\log^{2}\left(1+\frac{B}{A}b\right)+\text{Li}_{2}\left(\frac{A}{A+Bb}\right)-\right.\right.$
$\displaystyle~{}~{}\left.\left.\log\frac{B}{A}\log\left(1+\frac{B}{A}b\right)\right]\right]-\frac{1}{2}\beta^{2}\frac{1}{A}\left[-\text{Li}_{2}(b)+\left[\text{Li}_{2}\left(\frac{A+Bb}{A+B}\right)-\right.\right.$
$\displaystyle~{}~{}\left.\left.\text{Li}_{2}\left(\frac{A}{A+B}\right)+\log\frac{B}{A+B}\log\frac{A+Bb}{A}\right]\right],$

$\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\beta(1+\delta)\frac{b^{\beta}}{\beta}+\frac{\beta^{2}}{2}\text{Li}_{2}(b)+b^{2}\left(\frac{\beta}{4}+\frac{1}{32}\beta^{2}\right)+b\left(-\beta-\frac{5}{16}\beta^{2}\right)$
(41)
$\displaystyle~{}~{}+\log(1-b)\left(-\frac{9}{16}\beta^{2}+\frac{3}{4}\beta^{2}b-\frac{3}{16}\beta^{2}b^{2}\right)+\log
b\left(\frac{\beta^{2}b^{2}}{4}-\beta^{2}b\right),$ $\displaystyle I_{3}$
$\displaystyle=$
$\displaystyle\beta(1+\delta)\mathcal{B}(b,\beta,-2)+\left(\frac{\beta^{2}}{8}+\frac{\beta}{2}\right)\frac{-1+(1-b)^{3}+b(b+3-3b)}{2(1-b)^{3}}$
(42)
$\displaystyle~{}~{}-\frac{1}{2}\left(\beta+\frac{3}{4}\beta^{2}\right)\left(\frac{1}{(1-b)^{2}}-1\right)-\frac{\beta^{2}}{8}\frac{1-(1-b)^{2}+2\log(1-b)}{4(1-b)^{2}}$
$\displaystyle~{}~{}-\frac{3}{4}\beta^{2}\log(1-b)-\frac{\beta^{2}}{4}\frac{b(b-1-(b-2)\log
b)}{(1-b)^{2}}+\frac{\beta^{2}}{2}\frac{b\log b}{b-1}$
$\displaystyle~{}~{}-\frac{1}{2}\beta^{2}\left(-\text{Li}_{2}(b)-\frac{1}{2}\log^{2}(1-b)\right)-\frac{\beta}{8}\frac{b+\text{log}(1-b)}{1-b},$

and

$\displaystyle I_{4}$ $\displaystyle=$
$\displaystyle\beta(1+\delta)\mathcal{B}(b,\beta,-1)+\frac{\beta^{2}}{2}(\text{Li}_{2}(b)-\text{Li}_{2}(1-b)+\frac{\pi^{2}}{6})+\frac{b}{1-b}\left(-\frac{3}{4}\beta^{2}-\frac{\beta}{2}\right)$
(43)
$\displaystyle~{}~{}+\log(1-b)\left(\frac{\beta}{2}-\frac{3}{8}\beta^{2}\right)-\frac{b\log
b}{1-b}\frac{\beta^{2}}{2}+\log^{2}(1-b)\frac{\beta^{2}}{16}-\frac{\beta^{2}}{8}\frac{\log(1-b)}{1-b}.$

The special functions used in the above formulas are

$\text{Li}_{2}=-\int_{0}^{z}\frac{\log(1-u)}{u}\text{d}u,$ (44)
$\mathcal{B}(x,a,b)=\int_{0}^{x}t^{a-1}(1-t)^{b-1}\text{d}t,$ (45)

and

$\,{}_{2}F_{1}(a,b,c,x)=\sum_{n=0}^{\infty}\frac{(a)_{n}(b)_{n}}{(c)_{n}}\frac{x^{n}}{n!},$
(46)

where $\text{Li}_{2}$ is the Spence function, $\mathcal{B}$ is the beta
function and ${}_{2}F_{1}$ is the Gaussian hypergeometric function with the
Pochhammer symbol $(a)_{n}=\Gamma(a+n)/\Gamma(a)$ and $\Gamma(x)$ is the Gamma
function.

## References

  * Brambilla  _et al._ [2011] N. Brambilla _et al._ , Heavy quarkonium: progress, puzzles, and opportunities, Eur. Phys. J. C 71, 1534 (2011).
  * Wang  _et al._ [2004] P. Wang, C. Z. Yuan, X. H. Mo, and D. H. Zhang, The interference between virtual photon and $1^{--}$ charmonium in $e^{+}e^{-}$ experiment, Phys. Lett. B 593, 89 (2004).
  * Wang  _et al._ [2006] P. Wang, X. H. Mo, and C. Z. Yuan, Measurement of the exclusive light hadron decays of the $\psi^{\prime\prime}$ in $e^{+}e^{-}$ experiments, Int. J. Mod. Phys. A 21, 5163 (2006).
  * Guo and Yuan [2022] Y. P. Guo and C. Z. Yuan, Impact of the interference between the resonance and continuum amplitudes on vector quarkonia decay branching fraction measurements, Phys. Rev. D 105, 114001 (2022).
  * Actis  _et al._ [2010] S. Actis _et al._ , Quest for precision in hadronic cross sections at low energy: Monte carlo tools vs. experimental data, Eur. Phys. J. C 66, 585 (2010).
  * Dong  _et al._ [2018] X. K. Dong, L. L. Wang, and C. Z. Yuan, Derived Born cross sections of $e^{+}e^{-}$ annihilation into open charm mesons from CLEO-c measurements, Chin. Phys. C 40, 043002 (2018).
  * Kuraev and Fadin [1985] E. A. Kuraev and V. S. Fadin, On radiative corrections to the cross section for single photon annihilation of an $e^{+}e^{-}$ pair at high energy, Sov. J. Nucl. Phys. 41, 466 (1985).
  * Sun  _et al._ [2021] W. Y. Sun, T. Liu, M. Q. Jing, B. Zhong, and W. M. Song, An iterative weighting method to apply isr correction to $e^{+}e^{-}$ hadronic cross-section measurements, Frontiers of Physics 16, 64501 (2021).
  * Gribanov and Popov [2021] S. S. Gribanov and A. S. Popov, A new method for obtaining a born cross section using visible cross section data from $e^{+}e^{-}$ colliders, J. High Energ. Phys. 11, 203 (2021).
  * Wang  _et al._ [2023] Y. D. Wang, Y. N. Wang, and P. Wang, Analytical formula for the cross section of hadron production from $e^{+}e^{-}$ collisions around the narrow charmouinum resonances, arXiv 10.48550/arXiv.2311.13292 (2023).
  * Ablikim  _et al._ [2019] M. Ablikim _et al._ (BESIII Collaboration), Measurement of the phase between strong and electromagnetic amplitudes of $J/\psi$ decays, Physics Letter B. 791, 375 (2019).
  * Cahn [1987] R. N. Cahn, Analytic forms for the $e^{+}e^{-}$ annihilation cross section near the X including initial-state radiation, Phys. Rev. D 36, 2666 (1987).
  * Chen  _et al._ [1990] F. Z. Chen, P. Wang, J. M. Wu, and Y. S. Zhu, Radiative correction of the resonant cross section in $e^{+}e^{-}$ collision, Chin. Phys. C 14, 585 (1990).
  * Zhou  _et al._ [2017] X. Y. Zhou, Y. D. Wang, and L. G. Xia, Analytic forms for cross sections of di-lepton production from $e^{+}e^{-}$ collisions around the $J/\psi$ resonance, Chinese Phys. C 41, 083001 (2017).
  * Workman  _et al._ [2022] R. L. Workman _et al._ (Particle Data Group), Review of particle physics, Prog. Theor. Exp. Phys. 2022, 083C01 (2022).
  * Anashin  _et al._ [2012] V. V. Anashin _et al._ (KEDR collaboration), Measurement of $\Gamma_{ee}(J/\psi)\cdot\mathcal{B}(J/\psi\to e^{+}e^{-})$ and $\Gamma_{ee}(J/\psi)\cdot\mathcal{B}(J/\psi\to\mu^{+}\mu^{-})$, Phys. Lett. B 711, 280 (2012).
  * Blondel  _et al._ [2019] A. Blondel, J. Gluza, S. Jadach, P. Janot, and T. Riemann, eds.,  _Theory for the FCC-ee: Report on the 11th FCC-ee Workshop Theory and Experiments_, CERN Yellow Reports: Monographs, Vol. 3/2020 (CERN, Geneva, 2019).
  * Ablikim  _et al._ [2021] M. Ablikim _et al._ (BESIII Collaboration), Number of $J/\psi$ events at BESIII, Chin. Phys.C 46, 074001 (2021).
  * Ablikim  _et al._ [2024] M. Ablikim _et al._ (BESIII Collaboration), Determination of the number of $\psi(3686)$ events taken at BESIII, arXiv: 10.48550/arXiv.2403.06766 (2024).
  * Kou  _et al._ [2019] E. Kou _et al._ (Belle II Collaboration), The Belle II Physics Book, Prog. Theor. Exp. Phys. 2019, 123C01 (2019).
  * Peng  _et al._ [2020] H. P. Peng, Y. H. Zhen, and X. R. Zhou, Super tau-charm facility of china, Physics 49, 8 (2020).
  * Yuan and Karliner [2021] C. Z. Yuan and M. Karliner, Cornucopia of Antineutrons and Hyperons from a Super $J/\psi$ Factory for Next-Generation Nuclear and Particle Physics High-Precision Experiments, Phys. Rev. Lett. 127, 012003 (2021).