text:
# Total suppression of a large spin tunneling barrier in quantum adiabatic computation

## I Introduction

Recently a novel paradigm was suggested for the design of quantum algorithms for solving combinatorial search and optimization problems based on quantum adiabatic evolution Farhi . In the quantum adiabatic evolution algorithm (QAA) a quantum state closely follows a ground state of a specially designed slowly varying in time control Hamiltonian. At the initial moment of time the control Hamiltonian has a simple form with the known ground state that is easy to prepare; at the final moment of time it coincides with the “problem” Hamiltonian $`H_P`$ whose ground state encodes the solution of the classical optimization problem in question. It can also be chosen to reflect the bit-structure and cost spectrum of the problem. For example,

$`H_P={\displaystyle \underset{𝐳}{}}E_𝐳|𝐳𝐳|`$ (1)
$`|𝐳=|z_1_1|z_2_2\mathrm{}|z_n_n.`$
Here $`E_𝐳`$ is a cost function defined on a set of $`2^n`$ binary strings $`𝐳=\{z_1,\mathrm{},z_n\}`$ $`z_j=0,1`$, each containing $`n`$ bits. The summation in (1) is over $`2^n`$ states $`|𝐳`$ forming the computational basis of a quantum computer with $`n`$ qubits. State $`|z_j_j`$ of the $`j`$-th qubit is an eigenstate of the Pauli matrix $`\widehat{\sigma }_z`$ with eigenvalue $`12z_j=\pm 1)`$. If at the end of the QAA the quantum state is sufficiently close to the ground state of $`H_P`$ then the solution to the optimization problem can be retrieved by measurement.

Running of the algorithm for several NP-complete problems has been simulated on a classical computer using a large number of randomly generated problem instances that are believed to be computationally hard for classical algorithms FarhiSc ; FarhiSat ; FarhiCli . Results of these numerical simulations for relatively small size of the problem instances ( $`n`$ 20) suggest a quadratic scaling law of the run time of the quantum adiabatic algorithm with $`n`$. Furthermore, it was shown in Vazirani02  that the previous query complexity argument that lead to the exponential lower bound for unstructured search Bennett  cannot be used to rule out the polynomial time solution of NP-complete Satisfiability problem by QAA.

On the other hand, a set of examples of the 3-Satisfiability problem has been recently constructed Vazirani01 ; annealing ; Vazirani02  to test analytically the power of QAA in the situations where the optimization problem in question has multiple well-separated local minima.

In these examples the cost function $`E_𝐳`$ depends on a bit-string $`𝐳`$ with $`n`$ bits, $`𝐳=\{z_1,z_2,\mathrm{},z_n\}`$, only via the Hamming weight of the string, $`w_𝐳=(z_1+z_2+\mathrm{}+z_n)`$, so that $`E_𝐳=f\left(w_𝐳\right)`$. The function $`f(w)`$ is multi-modal, it has a local minimum separated from the global minimum by the barrier of an order-$`n`$ width in $`w`$. For that reason classical local search like simulating annealing provably fails to find a globally optimal solution in time polynomial in $`n`$. It annealing ; Vazirani02  a “standard” QAA was applied to this problem in which a quantum evolution begins in a uniform superposition state $`\frac{1}{\sqrt{2^n}}_𝐳|𝐳`$, ends in a target (solution) state $`|𝐳_t`$, and the control Hamiltonian is a linear interpolation in time between the initial and final Hamiltonians.

For the above examples it was shown annealing ; Vazirani02  that the system can be trapped during the QAA in the local minimum of a cost function for a time that grows exponentially in the problem size $`n`$. It was also shown annealing  that an exponential delay time in QAA can be computed in terms of a quantum-mechanical tunneling for an auxiliary large spin system.

It can also be inferred from Vazirani02 ; annealing  that QAA will have an exponential complexity even if the cost function $`E_𝐳`$ no longer depends strictly on a Hamming weight but the deviation only occurs for states $`|𝐳`$ that have exponentially small (in $`n`$) overlap with the adiabatic ground state wavefunction $`|\varphi _0(t)`$ at all times during the algorithm execution.

The above example has a significance more than just being a particular simplified case of the binary optimization problem with symmetized cost. Indeed one can argue that it shows one of the mechanisms for setting “locality traps” in the 3-Satisfiability problem Vazirani\_talk . But most importantly, this example demonstrates that exponential complexity of QAA results from a collective phenomenon in which transitions between the bit-configurations with low-lying energies can only occur by the simultaneous flipping of large clusters containing order-n bits. In many cases these transitions can be analyzed as a tunneling of spin variables. A similar phenomenon related to the tunneling of magnetization was recently observed in the large-spin molecular nanomagnets Wernsdorfer .

However low-energy collective behavior is also well known in spin glass models, many of which are in one-to-one correspondence with random NP-complete problems Anderson . In particular, an important ingredient of the “replica symmetry breaking” picture of an infinite-range spin glass by Parizi Parizi  is that there are collective spin exitations that are of the order of the system size $`n`$ whose energy is $`(𝒪(1)`$, i.e., it does not grow with the size of the system. A similar picture may be applicable to random Satisfiability problems Monasson .

Therefore in connection to the above example, it is important to understand how to design a polynomial time QAA without a prior knowledge of a particular form of the cost function $`f\left(w\right)`$ (possibly multi-modal or even randomly sampled), so that a tunneling barrier between the local minima will be totally suppressed. This is a focus of the present paper footnote1 .

In Sec. III we present a theory of quantum adiabatic evolution of a large spin system introducing a control Hamiltonian that guarantees polynomial time complexity for QAA in the symmetrized 3-Satisfiability example mentioned. In Sec. IV we provide an estimate of the minimum gap between the two lowest eigenvalues of the control Hamiltonian that determines the complexity of QAA. Sec. V contains concluding remarks.

## II  Quantum Adiabatic Evolution Algorithm

In a standard QAA Farhi  one specifies the time-dependent Hamiltonian $`H(t)=\stackrel{~}{H}(t/T)`$

$$\stackrel{~}{H}(\tau )=(1\tau )H_D+\tau H_P,$$
(2)
where $`\tau =t/T`$ is dimensionless “time”. This Hamiltonian guides the quantum evolution of the state vector $`|\psi (t)`$ according to the Schrödinger equation $`i|\psi (t)t=H(t)|\psi (t)`$ from $`t=0`$ to $`t=T`$, the run time of the algorithm (we let $`\mathrm{}=1`$). $`H_P`$ is the “problem” Hamiltonian given in (1). $`H_D`$ is a “driver” Hamiltonian, that is designed to cause the transitions between the eigenstates of $`H_P`$. In this algorithm one prepares the initial state of the system $`|\psi (0)`$ to be a ground state of $`\stackrel{~}{H}(0)=H_D`$. It is typically constructed assuming $`\mathrm{𝑛𝑜}`$ knowledge of the solution of the classical optimization problem and related ground state of $`H_P`$. In the simplest case

$$H_D=C\underset{j=0}{\overset{n1}{}}\widehat{\sigma }_x^j,|\psi (0)=2^{n/2}\underset{𝐳}{}|𝐳,$$
(3)
where $`\sigma _x^j`$ is a Pauli matrix for $`j`$-th qubit and $`C>0`$ is some scaling constant. Consider instantaneous eigenstates $`|\varphi _k(\tau )`$ of $`\stackrel{~}{H}(\tau )`$ with energies $`E_k(\tau )`$ arranged in nondecreasing order at any value of $`\tau (0,1)`$

$$\stackrel{~}{H}|\varphi _k=E_k|\varphi _k,k=0,1,\mathrm{},2^n1.$$
(4)
Provided the value of $`T`$ is large enough and there is a finite gap for all $`t(0,T)`$ between the ground and exited state energies, $`\mathrm{\Delta }E(\tau )=E_1(\tau )E_0(\tau )>0`$, quantum evolution is adiabatic and the state of the system $`|\psi (t)`$ stays close to an instantaneous ground state, $`|\varphi _0(t/T)`$ (up to a phase factor). Because $`H(T)=H_P`$ the final state $`|\psi (T)`$ is close to the ground state $`|\varphi _0(\tau =1)`$ of the problem Hamiltonian. Therefore a measurement performed on the quantum computer at $`t=T(\tau =1)`$ will find one of the solutions of combinatorial optimization problem with large probability. Quantum transition away from the adiabatic ground state occurs most likely in the vicinity of the point $`\tau \tau _c`$ where the energy gap $`\mathrm{\Delta }E(\tau )`$ reaches its minimum (avoided-crossing region). The probability of the transition, $`1|\psi (t)|\varphi _0(t/T)|_{t=T}^2`$, is small provided that

$$T\frac{|\varphi _1|\stackrel{~}{H}_\tau |\varphi _0|_{\tau =\tau _c}}{\mathrm{\Delta }E_{\mathrm{min}}^2},\mathrm{\Delta }E_{\mathrm{min}}=\underset{0\tau 1}{\mathrm{min}}\left[E_1(\tau )E_0(\tau )\right],$$
(5)
($`\stackrel{~}{H}_\tau d\stackrel{~}{H}/d\tau `$). The fraction in (5) gives an estimate for the required runtime of the algorithm and the task is to find its asymptotic behavior in the limit of large $`n1`$. The numerator in (5) is less than the largest eigenvalue of $`\stackrel{~}{H}_\tau =H_PH_D`$, typically polynomial in $`n`$ Farhi . However, $`\mathrm{\Delta }E_{\mathrm{min}}`$ can scale down exponentially with $`n`$ and in such cases the runtime of quantum adiabatic algorithm will grow exponentially fast with the size of the input $`n`$.

### II.1 Binary optimization problems with symmetrized cost function

Consider a cost function $`E_𝐳`$ in the following form:

$$E_𝐳=f\left(w_𝐳\right),w_𝐳=\underset{j=1}{\overset{n}{}}z_j.$$
(6)
This cost is symmetric with respect to the permutation of bits and $`w_𝐳`$ is a Hamming weight of a string $`𝐳`$. A particular example of this problem related to 3-Satisfiability was introduced in Vazirani01 ; annealing ; Vazirani02  (the discussion in this subsection closely follows annealing , Sec. 4). In this example

$`E_𝐳={\displaystyle \underset{i<j<k}{}}c(z_i+z_j+z_k),`$
$`c(m)=(1\delta _{m,0})\left(1+\delta _{m,1}(q1)\right),m=0,1,2,3`$ (7)
where $`\delta _{k,l}`$ is a Kronecker delta. For this particular case function $`f(w)`$ in (6) takes the following form:

$$f(w)=\frac{q}{2}w(nw)(nw1)+\frac{1}{2}w(w1)(nw)+\frac{1}{6}w(w1)(w2).$$
(8)
where $`q`$ is an integer greater than or equal to 3. In the leading order in $`n1`$ one can write:

$$f(w)=\left(\frac{n}{3}\right)^3g\left(\frac{w}{n}\right)+𝒪(n^2),$$
(9)
where

$$h(u)=4qu(1u)^2+4u^2(1u)+\frac{4}{3}u^3.$$
(10)
$`g(u)`$ is a non-monotonic function with global minimum at $`u=0`$ corresponding to $`z_1,z_2,\mathrm{},z_n=0`$. It also has a local minimum at $`u=1`$ corresponding to $`z_1,z_2,\mathrm{},z_n=1`$ ( cf. Ref. annealing , Fig. 1).

In QAA the symmetrized cost function (8) corresponds to the following problem Hamiltonian $`H_P`$ of total spin-$`n/2`$ system (cf. Eq. (1) )

$$H_P=f\left(\frac{n}{2}\widehat{S}_z\right),$$
(11)
where

$$\widehat{S}_z|𝐳=S_z)|𝐳,S_z=\frac{n}{2}w_𝐳.$$
(12)
Here $`\widehat{S}_z`$ is operator of $`z`$-projection of a total spin $`\frac{n}{2}`$ of a system of $`n`$ individual spins $`1/2`$. We used an obvious connection between the values of the Hamming weight function $`w_𝐳`$ and corresponding eigenvalues $`S_z`$ of the operator $`\widehat{S}_z`$. In what following we use “hat” notation for the operators of a total spin.

It was shown annealing  that if one uses $`H_D`$ in a form corresponding to (3)

$$H_D=\left(\genfrac{}{}{0pt}{}{n1}{2}\right)\left(\frac{n}{2}\widehat{S}_x\right).$$
(13)
then the minimum gap $`\mathrm{\Delta }E_{\mathrm{min}}`$ scales down exponentially with $`n`$ implying the exponential complexity of QAA for this problem.

## III  Adiabatic Evolution of a Large Spin

### III.1 Extended vs localized initial states

In this paper we show that the main reason for the exponentially small minimum gap for the problem (7) with $`H_D`$ given in (13) is the fact that the ground state of the driving Hamiltonian is localized in the space of the $`z`$-projections of a total spin. We construct an example of the driving Hamiltonian with extended ground state and show that in this case the evolution time of QAA is polynomial in the number of qubits.

In the following we adopt the notation for the total spin $`\widehat{𝐒}=\left\{\widehat{S}_j\right\}`$, with $`j=x,y,z`$, where $`\widehat{S}_j`$ are the projections of the total spin operator on the $`j`$-th axis. The total spin operator’s components equal $`\widehat{S}_j=_{k=1}^n\widehat{S}_j^{\left(k\right)}`$ and are symmetric in all one-qubit spin operators $`\left\{\widehat{S}_j^{\left(k\right)}\right\}`$ for $`k=1,\mathrm{}n`$. We use $`S_j`$ for the eigenvalues of the operator $`\widehat{S}_j`$.

We now introduce a new driver Hamiltonian

$$H_D=\widehat{S}_x^2.$$
(14)
Consider, for example, its ground state in the case when the total number of qubits $`n=2s`$ is even and therefore the total spin $`l=n/2=s`$ is an integer. Then the ground state is a state with $`S_x=0`$ in a basis where $`x`$ is chosen as a quantization axis. Making use of Wigner’s rotation matrix $`d_{m,m^{}}^l\left(\theta \right)`$ LANDAU1 , one easily projects this state onto the computational (problem) basis with the $`z`$ quantization axis and $`m=S_z`$. This gives us the ground state wave function in the problem basis as $`\mathrm{\Psi }_0\left(m\right)=d_{m,0}^l\left(\pi /2\right)`$. In the limit of large spin (which is the case of interest for us), the wave function is given by

$$\mathrm{\Psi }_0\left(m\right)\frac{\left(1\right)^k}{\sqrt{2\pi }}\frac{1}{\left(l^2m^2\right)^{1/4}}\delta _{lm,2k},l=\frac{n}{2}(n\mathrm{is}\mathrm{even}).$$
(15)
Note that $`\mathrm{\Psi }_0\left(m\right)=0`$ for the odd difference $`lm=2k+1`$ and is given by (15) for $`lm=2k`$. From (15), it follows that the ground state of $`S_x^2`$ is delocalized in the $`m`$\- space and spread over the whole interval, $`lml`$.

On the other hand, the ground state of the operator (13) used in annealing  as a driver, is localized in the $`m`$-space. Indeed, that ground state is a state with $`S_x=l`$ in the $`x`$ basis. This gives us the ground state wave function in the problem basis, $`\mathrm{\Psi }_0^{}\left(m\right)=d_{m,l}^l\left(\pi /2\right)`$. In the limit of large spin one has

$$\mathrm{\Psi }_0^{}\left(m\right)\frac{1}{\sqrt{\pi l}}\mathrm{exp}\left(\frac{m^2}{l}\right),$$
(16)
which is clearly localized on the scale $`m^2l`$. The same conclusion holds for the case when $`n`$ is odd and $`l`$ take half-integer values.

As we will show below, the adiabatic evolution of the delocalized states is fundamentally different from the localized ones. In particular, the localized ground states in general result in macroscopic tunneling. This gives the exponentially small in $`n`$ ground state energy gaps and consequently the exponentially large complexity of QAA. Using the delocalized states one can avoid macroscopic tunneling. We argue that in a general situation when the information about the ground state is not used for constructing the driver Hamiltonian, the driver with extended ground state should lead to polynomial complexity of adiabatic algorithms independently of the specific form of the problem Hamiltonian provided that it is expressed as a function (in general, nonlinear) of the total spin operators, $`\widehat{S}_j`$ (cf. Eq. (11)).

### III.2 WKB approximation for the large spin

To be specific, we refer to the same problem Hamiltonian as in Eq. (11). The full Hamiltonian takes the following form in the limit of large $`n`$

$`H\left(\tau \right)=\left(1\tau \right)n\widehat{S}_x^2+\tau \left({\displaystyle \frac{n}{2}}\right)^3h\left(\widehat{u}\right),`$ (17)
$`\widehat{u}={\displaystyle \frac{1}{2}}\left(1{\displaystyle \frac{\widehat{S}_z}{l}}\right),l={\displaystyle \frac{n}{2}},`$
where $`\tau =t/T`$ ; $`0\tau 1`$ and function $`h(u)`$ is given in (10).

In order to get a simple physical picture, we will refer to the WKB-type approach commonly used in the theory of quantum spin tunneling in magnetics (QTM) CHUD1 , GARG1 , GARG2 . This approach is applicable for large spins $`(l>>1)`$, which is the case of interest for us. We choose $`x`$ as a quantization axis and following the standard procedure obtain the effective quasi-classical Hamiltonian in polar coordinates $`\{\theta ,\phi \}`$ with $`\theta [0,\pi ]`$ and $`\phi [0,2\pi ]`$. In doing this, we make use of the following relations valid in the limit $`(l>>1)`$

$$S_z^n\left(l\mathrm{sin}\theta \mathrm{sin}\phi \right)^n.$$

As it was shown in GARG2 , in the quasi-classical limit we have

$$\widehat{S}_x^2\frac{^2}{\phi ^2}\widehat{p}^2,\theta \frac{\pi }{2},$$
(18)
which means that the motion of the quasi-classical spin can be described as a 1D motion of a massive particle on a unit radius ring in the appropriate effective potential $`V(\varphi )`$. We also have

$$u=\frac{1}{2}\left(1\mathrm{sin}\phi \right).$$
(19)
Substituting this into (17), we finally obtain the effective quasi-classical Hamiltonian in the form

$`H\left(\tau \right)=n\left(1\tau \right)\widehat{p}^2+\tau \left({\displaystyle \frac{n}{2}}\right)^3V(\phi ),`$
$`V(\phi )=h\left({\displaystyle \frac{1\mathrm{sin}\phi }{2}}\right).`$ (20)
One should note that the driver $`H_D=n\widehat{S}_x^2`$ does not change the effective potential caused by $`H_P`$ and only introduces the effective kinetic energy into the problem. As we show below, this driver has extended eigenstates in the space of the $`z`$-projections $`m`$ of a total spin, as opposed to the case $`H_D=n^2S_x`$ considered in annealing , where the eigenstates are localized in $`m`$. As we will see below, this leads to the absence of tunneling and polynomial gap in our case as opposed to the exponentially small tunneling amplitude arising in annealing  . The adiabatic wave functions $`\mathrm{\Psi }_k\left(\phi \right)`$ satisfy

$$\left[\left(1\tau \right)\frac{}{\phi ^2}+\tau \frac{n^2}{8}V(\phi )\right]\mathrm{\Psi }_k\left(\phi \right)=\stackrel{~}{E}_k\left(\tau \right)\mathrm{\Psi }_k\left(\phi \right),$$
(21)
where $`\stackrel{~}{E}_k\left(\tau \right)=E_k\left(\tau \right)/n`$ is a rescaled energy. The first term in the l.h.s. of (21) corresponds to the driving Hamiltonian $`H_D`$ and presents the kinetic energy of the particle moving in the effective potential due to the problem Hamiltonian $`H_P.`$ At initial moment $`\tau =0`$, the total Hamiltonian reduces to the driver $`H_D`$ which describes a free massive particle on a ring and its eigenstates are well known LANDAU1 . In this case, the Schrödinger equation (21) gives exact wavefunctions and spectrum. Namely, we have

$`\mathrm{\Psi }_k\left(\phi \right)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}\mathrm{exp}\left(ik\phi \right),`$ (22)
$`\stackrel{~}{E}_k\left(0\right)`$ $`=`$ $`k^2,`$
where $`k`$ is an integer number $`k=0,\pm 1,\pm 2,\mathrm{}`$ for even $`n`$ when the total spin $`l=n/2`$ is integer, and a half-integer number $`k=\pm 1/2,\pm 3/2,\mathrm{}`$ for an odd number of qubits $`n`$, when $`l=[n/2]+1/2`$. One should note that in case of integer spin $`l`$ , all eigenstates are twofold degenerate except for the ground state corresponding to $`k=0`$ and for the half-integer spin, all states including the ground state are twofold degenerate. This is a particular case of the Kramers’ degeneracy LANDAU1 , GARG1 , which occurs due to the symmetry with respect to the spin-flip transformation. Note that since from (19) it follows that $`m=l\mathrm{sin}\phi `$, the states (22) are indeed delocalized in the $`m`$-space.

From (21), it follows that the problem Hamiltonian is of the same order of magnitude as the driving one only for sufficiently short times $`\tau n^21`$. This means that the problem Hamiltonian is of the order of the level separation of the kinetic term only for sufficiently small times. Since the eigenstates of the driving Hamiltonian are extended (delocalized) in $`m`$, this implies the following qualitative picture of adiabatic evolution of the ground energy level with $`\tau `$ . At sufficiently small times $`\tau n^21`$, the term due to $`H_P`$ can be considered as a perturbation to the driver term and the energy levels are not strongly distorted. The eigenstates are delocalized in the $`m`$ space. As $`\tau `$ increases, the ground state is affected by the perturbation and the ground state energy increases.

### III.3 Minimum gap analysis

If the ground state is not degenerate (this is the case when the total number of qubits $`n=2s`$ is even), the gap may have a non-monotonic behavior in $`\tau `$ in the range $`\tau n^21`$. Qualitatively, this can be described as follows. For sufficiently small times $`\tau n^21`$, the ground state energy is increasing in $`\tau `$ due to the diagonal matrix elements of the problem Hamiltonian until compensated by the level repulsion from the first excited state. After this, the ground state is pushed down and gradually approaches the ground state of $`H_P`$ (which is $`E_g=0`$ in the present case). Since the upper levels are not as strongly affected by the perturbation as the ground state and since the number of upper levels is very large, $`𝒪(n)`$, the system of upper levels behaves as a rigid one. For this reason the strong interaction between the ground and first excited states occurs in the range of energies $`\mathrm{\Delta }E`$ corresponding to the the level separation in the driver Hamiltonian and is linear in $`n`$. Therefore

$$\mathrm{\Delta }E=𝒪\left(n^2\mathrm{\Delta }E_P\right),$$
(23)
where $`\mathrm{\Delta }E_P`$ denotes the separation between the two lowest eigenvalues of $`H_P`$. Clearly, the energy scale of the minimal gap is polynomial in the number of qubits $`n`$ provided that $`\mathrm{\Delta }E_P`$ is polynomial.

As we have discussed above, in case of an odd number of qubits, all initial eigenstates of $`H_D`$ have a twofold degeneracy (Kramers’ degeneracy). In this case, the degeneracy is removed by the effective potential (Zeeman splitting) and there is no avoided-crossing between the ground state and the first excited states of the total Hamiltonian $`H(\tau )`$. One can easily show that in this case the level separation $`E_1(\tau )E_0(\tau )`$ grows monotonically in $`\tau `$ LANDAU1  (it is linear in $`\eta `$ for $`\eta 1`$).

Apart from the features occurring in the range of small $`\tau n^2`$, the evolution of the ground state on the large time scale $`\tau 1`$ is the same for both even and odd number of qubits. Most importantly, the estimate (23) holds globally in both cases. As we will see below, this is confirmed by numerical simulations.

## IV Minimal Gap Estimate

As we discussed above, the ground state energy may have a non-monotonic time dependence on $`\tau `$ when the total number of qubits is even and the ground state of $`H_D`$ is not degenerate. Because the eigenvalues of the driver Hamiltonian, $`nk^2`$, grow rapidly with the quantum number $`k`$ it is possible to obtain a minimum gap estimate analyzing how $`H_P`$ affects the two lowest eigenvalues of the driver. Making use of (2), we obtain the adiabatic gap as a function of time in the range $`\tau n^21`$

$$E_1(\tau )E_0(\tau )\sqrt{\left[\left(1\tau \right)n\tau \left(H_P\right)_{e_1,e_2}\right]^2+8\tau ^2\left|\left(H_P\right)_{g,e}\right|^2},$$
(24)
Here we explicitely take into account that the first excited level of $`H_D`$ is twofold degenerate. Ssubscripts $`g`$, $`e_{1,2}`$ above denote the ground state of $`H_D`$ and the two lowest exited states, respectively. Matrix elements of $`H_P`$ on these states satisfy the following relations:

$$\left(H_P\right)_{g,g}=\left(H_P\right)_{e,e},\left(H_P\right)_{g,e_1}=\left(H_P\right)_{g,e_2}\left(H_P\right)_{g,e}.$$
(25)
From (24), we obtain an estimate for the time when the minimal gap is achieved

$$\tau _c=\frac{n\left[\left(H_P\right)_{e_1,e_2}+n\right]}{\left|\left(H_P\right)_{g,e}\right|^2+\left[\left(H_P\right)_{e_1,e_2}+n\right]^2},$$
(26)
In the limit $`n1`$ we can write

$$\eta _c=\tau _cn^2\frac{n^3\left(H_P\right)_{e_1,e_2}}{\left|\left(H_P\right)_{g,e}\right|^2+\left[\left(H_P\right)_{e_1,e_2}\right]^2}.$$
(27)
Since matrix elements of $`H_Pn^3`$, it can be easily verified that the scaled quantity $`\eta _c`$ does not depend on $`n`$. The matrix elements can be calculated either in the quasiclassical basis given by (22) or exactly using (17) (see Appendix). For the Hamiltonian (17), the quasiclassical matrix elements yield

$$\left[\left(H_P\right)_{e_1,e_2}\right]^2=\frac{9}{1024}n^6.$$
(28)
(for $`q=3`$). Therefore, we have $`\eta _{\mathrm{min}}=\frac{64}{9}7.1`$. This is in a qualitative agreement with $`\eta _{\mathrm{min}}5.0`$ obtained in our numerical simulations. Substituting into (24), we obtain the estimate for the minimal gap

$$\mathrm{\Delta }E_{\mathrm{min}}=\mathrm{\Delta }E\left(\tau _c\right)=n\left(\frac{2}{3}\right)^{1/2}.$$
(29)
The corresponding value of the slope $`\mathrm{\Delta }E_{\mathrm{min}}/n=\sqrt{\frac{2}{3}}0.82`$ is again in a qualitative agreement with the value $`0.86`$ obtained in the numerical simulations (Fig.5).

### IV.1 Numerical analysis

We also performed numerical simulations of the adiabatic spectrum with Hamiltonian (17). In Fig. 2, we plot the ground state and the first excited state energies as functions of the dimensionless ”time” $`\eta =\tau n^2`$ for even values of $`n`$. According to the above discussion, in this case the ground state is not degenerate and the energy level differnce exhibits a non-monotonic behavior (cf. Fig. 3) leading to the polynomial minimal gap $`\mathrm{\Delta }E_{\mathrm{min}}n`$ at $`\eta 1`$.

In Fig. 4 we plot an eigenvalue spectrum of the Hamiltonian $`H(\tau )`$ at the avoided crossing point, $`\tau =\tau _c`$ for an even value of $`n`$. For $`k1`$ $`E_k=𝒪(n^3k^2)`$ which corresponds to the scaling analysis in Eqs. (31),(32). Zeeman splitting of a doublet of the two lowest exited energy levels in Fig. 4 is of the order of $`n^3\tau _c=𝒪(n)`$.

In Fig.3 we plot the dependence of minimal gap $`\mathrm{\Delta }E_{\mathrm{min}}`$ vs $`n`$ for even values of $`n`$ . The insert to this plot corresponds to the dependence of the avoided crossing point $`\tau _c`$ vs $`n`$.

Using the simple quasi-classical picture presented above it is not difficult to compute the minimum gap in the case of the driver Hamiltonian (13) in terms of the appropriate tunneling exponent and recover the answer given in annealing .

## V  Conclusion

We show that macroscopic tunneling in QAA with the symmetrized cost function can be totally suppressed using a driver Hamiltonian with the ground state extended in the space of a total spin projection onto the direction of computational basis. This leads to a polynomial time complexity of QAA. We give a simple form of a driver Hamiltonian (14) that has the aforementioned property and therefore makes the algorithm polynomial.

We developed a simple intuitive picture of this phenomenon based on WKB approximation for large spins. It follows from the analysis that the results will hold for any form of $`f(w)`$ that is sufficiently smooth on the scale $`w1`$. We conjecture that this phenomenon holds even if the cost function $`E_𝐳`$ is symmetrical with respect to the permutation of bits but is not expressed only throught the Hamming weight of a string.

We also argue using a picture of QAA as a quantum local search Vazirani02 ; ST  that suppression of tunneling barriers with the operator $`S_x^2`$ changes the global neighborhood properties for QAA in a very profound way that can have an effect on the algorithm complexity for a larger class of cost functions, $`E_𝐳=f(w_𝐳)+\mathrm{\Delta }f_𝐳`$ where $`\mathrm{\Delta }f_𝐳`$ breaks the symmetry between the bits. In particular, it can be large for those states $`|𝐳^{}`$ that have exponentially small overlap $`|𝐳^{}|\mathrm{\Psi }_0(\tau )|`$ with adiabatic ground states at all times (cf. Sec. 7.3 of Vazirani02 ).

A possible generalization of the above analysis is related to random optimization problems with frustration, such as NP-hard problems and corresponding spin glass models. Exponential complexity of quantum adiabatic evolution algorithms for these problems is not necessarily related to tunneling but rather to the quantum diffusion phenomenon associated with the rapid falloff of correlations in the bit-structure with growing size of neighborhood around a given string ST . The role of the tunneling and collective phenomena involving the low cost configurations in the performance of the quantum adiabatic evolution algorithms in random NP-hard is yet to be analyzed.

## VI Acknowledgments

We wish to thank Edward Farhi, Samuel Gutmann, Andrew Childs ( MIT), and Umesh Vazirani ( UC Berkley) for stimulating discussions. We also wish to thank Alex Burin (Northwestern University) for turning our attention to Refs. GARG1 ; GARG2 . This research was supported by NASA Intelligent Systems Revolutionary Computing Algorithms program (project No:749-40).

## VII  Note Added

Shortly after this work was completed, we learned about the work of E.Farhi, J. Goldstone, S. Gutmann farhi\_paths  wherein another approach was suggested to achieve polynomial complexity of the quantum adiabatic evolution algorithm for the same optimization problem with the symmetrized cost function (8). The algorithm proposed in Ref. farhi\_paths  is quite different from our algorithm. We believe that the proposal in farhi\_paths  for random generation of interpolating “paths” $`H(t)`$ in different trials of the algorithm is a very promising tool. It can perhaps be modified, by including intermediate measurents farhimeas , to become an efficient adaptive algorithm for random NP-hard optimization problems.

However we believe the particular method of achieving polynomial complexity presented in farhi\_paths  is less robust than ours in the problems with $`n`$ qubits where the dynamics of the total spin $`n/2`$ is a key. The difference between the two algorithms is that in our approach a simple universal form of the driver Hamiltonian guarantees that the minimum gap scales polynomially in the problem size for a broad class of symmetrized cost functions $`f(w)`$. Also, by construction our method does not require a specific knowledge of the solution, or even a specific form of the cost function $`f(w)`$.

We also believe that simple quasi-classical approach presented above enables one to study analytically a “volume” in the space of possible interpolating paths farhi\_paths  that find a solution in polynomial time. Detailed study of this issue will be presented elsewhere.

## VIII Appendix: Minimum gap estimate using matrix representation

The eigenstates and eigenvectors of the Hamiltonian (17) can be obtained solving from the matrix representation. Choosing $`x`$ as a quantization axis, we have $`S_z=\left(i\right)\left(S_+S_{}\right)`$ and

$`\left(S_x^2\right)_{m,m}`$ $`=`$ $`m^2,m=l,l+1,\mathrm{},l,`$
$`\left(S_+\right)_{m,m1}`$ $`=`$ $`\left(S_{}\right)_{m1,m}={\displaystyle \frac{1}{2}}\sqrt{\left[l\left(l+1\right)m\left(m1\right)\right]}.`$
The full Hamiltonian (17) is given by

$`H\left(\tau \right)`$ $`=`$ $`\left(1\tau \right)nS_x^2+\tau \left({\displaystyle \frac{n}{2}}\right)^3h\left(S_z\right),`$
$`h\left(S_z\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\left(q+{\displaystyle \frac{4}{3}}\right)+{\displaystyle \frac{\left(2q\right)}{l}}S_z{\displaystyle \frac{q}{l^2}}S_z^2+{\displaystyle \frac{\left(q2/3\right)}{l^3}}S_z^3\right],`$
and has the following matrix elements

$`\left(H\right)_{m,m}`$ $`=`$ $`\left(1\tau \right)nm^2+\tau \left({\displaystyle \frac{n}{2}}\right)^3{\displaystyle \frac{1}{2}}\left(q+{\displaystyle \frac{4}{3}}q\mathrm{\Lambda }_0\right),`$
$`\left(H\right)_{m,m1}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\tau \left({\displaystyle \frac{n}{2}}\right)^3\left[\left(2q\right)\mathrm{\Lambda }_1+\left(q2/3\right)\mathrm{\Lambda }_2\right],`$ (30)
$`\left(H\right)_{m,m2}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\tau \left({\displaystyle \frac{n}{2}}\right)^3\left(q\right)\mathrm{\Lambda }_3,`$
$`\left(H\right)_{m,m3}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\tau \left({\displaystyle \frac{n}{2}}\right)^3\left(q2/3\right)\mathrm{\Lambda }_4,`$
with

$`\mathrm{\Lambda }_0(l,m)`$ $`=`$ $`{\displaystyle \frac{1}{2l^2}}\left[l\left(l+1\right)m^2\right],`$
$`\mathrm{\Lambda }_1(l,m)`$ $`=`$ $`{\displaystyle \frac{\left(i\right)}{2l}}\sqrt{l\left(l+1\right)m\left(m1\right)},`$
$`\mathrm{\Lambda }_2(l,m)`$ $`=`$ $`{\displaystyle \frac{3i}{4l^3}}\sqrt{l\left(l+1\right)m\left(m1\right)}\left[l\left(l+1\right)m\left(m1\right)1\right],`$
$`\mathrm{\Lambda }_3(l,m)`$ $`=`$ $`{\displaystyle \frac{1}{4l^2}}\sqrt{\left[l\left(l+1\right)m\left(m1\right)\right]\left[l\left(l+1\right)\left(m1\right)\left(m2\right)\right]},`$
$`\mathrm{\Lambda }_4(l,m)`$ $`=`$ $`{\displaystyle \frac{i}{2l}}\sqrt{\left[l\left(l+1\right)\left(m2\right)\left(m3\right)\right]}\mathrm{\Lambda }_3(l,m).`$
One should note that the coefficients $`\left\{\mathrm{\Lambda }_k\right\}`$ do not scale with $`l=2n`$, i.e. $`\mathrm{\Lambda }_k=𝒪(1)`$ for $`n1`$. From (30), it follows that the eigenvalue problem $`H\left(\tau \right)\mathrm{\Psi }_k=E_k\left(\tau \right)\mathrm{\Psi }_k`$ is rescaled for $`n11`$ in terms of dimensionless variables $`\eta =\tau n^2`$ and $`\stackrel{~}{E}_k\left(\tau \right)=E_k\left(\tau \right)/n`$ as

$$\stackrel{~}{H}\left(\tau \right)\mathrm{\Psi }_k=\stackrel{~}{E}_k\left(\tau \right)\mathrm{\Psi }_k,$$
(31)
with

$`\left(\stackrel{~}{H}\right)_{m,m}`$ $`=`$ $`m^2+{\displaystyle \frac{\eta }{16}}\left(q+{\displaystyle \frac{4}{3}}q\mathrm{\Lambda }_0\right),`$
$`\left(\stackrel{~}{H}\right)_{m,m1}`$ $`=`$ $`{\displaystyle \frac{\eta }{16}}\left[\left(2q\right)\mathrm{\Lambda }_1+\left(q2/3\right)\mathrm{\Lambda }_2\right],`$ (32)
$`\left(\stackrel{~}{H}\right)_{m,m2}`$ $`=`$ $`{\displaystyle \frac{\eta }{16}}\left(q\right)\mathrm{\Lambda }_3,`$
$`\left(\stackrel{~}{H}\right)_{m,m3}`$ $`=`$ $`{\displaystyle \frac{\eta }{16}}\left(q2/3\right)\mathrm{\Lambda }_4,`$
implying that matrix elements of $`\stackrel{~}{H}`$ do not scale with $`n`$ and neither do the eigenvalues of $`\stackrel{~}{H}`$. It is strighforward to obtain an estimate for the minimum gap from the above equation using a Brillouin-Wigner perturbation theory in parameter $`\eta `$. If one uses a 3-level truncation scheme, ($`m=0,\pm 1`$), and sets $`q=3`$ then a quasi-classical estimate (29) can be recovered to the leading order in $`n`$.


format: text

text:
# AN UNFOLDING METHOD FOR HIGH ENERGY PHYSICS EXPERIMENTS

## 1 THE UNFOLDING PROBLEM

A standard task in high energy physics experiments is the measurement of a distribution $`f(x)`$ of some kinematical quantity $`x`$. With an ideal detector one could measure the quantity $`x`$ in every event and could obtain $`f(x)`$ by a simple histogram of the quantity $`x`$. With *real* detectors the determination of $`f(x)`$ is complicated by three effects:

* Limited acceptance: The probability to observe a given event, the *detector acceptance*, is less than 1. The acceptance depends on the kinematical variable $`x`$.
* Transformation: Instead of the quantity $`x`$ a different, but related quantity $`y`$ is measured. The transformation from $`x`$ to $`y`$ can be caused by the non-linear response of a detector component.
* Finite resolution: The measured quantity $`y`$ is smeared out due to the finite resolution (or limited measurement accuracy) of the detector. Thus there is only a statistical relation between the true kinematical variable $`x`$ and the measured quantity $`y`$.

The really difficult effect in the data correction for experimental effects, or data transformation from $`y`$ to $`x`$ is the finite resolution, causing a *smearing* of the measured quantities. Mathematically the relation between the distribution $`f(x)`$ of the true variable $`x`$, to be determined in an experiment, and the measured distribution $`g(y)`$ of the measured quantity $`y`$ is given by the integral equation,

$$\overline{)g(y)=A(y,x)f(x)\text{d}x},$$
(1)
called a Fredholm integral equation of the first kind. In practice often a known (measured or simulated) background contribution $`b(y)`$ has to be added to the right-hand side of equation (1); this contribution is ignored in this paper. The resolution function $`A(y,x)`$ represents the effect of the detector. For a given value $`x=x_0`$ the function $`A(y,x_0)`$ describes the response of the detector in the variable $`y`$ for that fixed value $`x_0`$. The problem to determine the distribution $`f(x)`$ from measured distributions $`g(y)`$ is called *unfolding*; it is called an inverse problem. Unfolding of course requires the knowledge of the resolution function $`A(y,x)`$, i.e. all the effects of limited acceptance, transformation and finite resolution.

In addition to the imperfections of the detector, there may be further effects between $`x`$ and $`y`$, which are *outside of the experimental control*, even with an ideal detector. One example are radiative effects, which in experiments are often corrected afterwards *(radiative corrections)*, but are in their effect similar to detector effects. If the true kinematical quantity is defined at the *parton level*, further effects from the fragmentation process of partons to the (observable) hadrons influence the measured quantity $`y`$. All these effects are of statistical nature.

A typical example for these effects is shown in Figure 1 taken from a Monte Carlo simulation of all three effects. By unfolding an estimate of the original distribution has to be determined from the distorted measured distribution. Details on the Monte Carlo simulation are given later in section 3, where the unfolding in this example is discussed in detail.

For the numerical solution of equation (1) the distributions have to represented by a finite set of parameters. One posssibility is to represent the distributions by histograms, and the resolution function by a matrix. Equation (1) can then be represented by the matrix equation

$$\overline{)𝒚=𝑨𝒙},$$
(2)
which has to be solved for the vector $`𝒙`$, given the vector $`𝒚`$ (data histogram). The vector $`𝒚`$ with $`n`$ elements represents a histogram of the measured quantity $`y`$, and the distribution $`f(x)`$ is represented by a histogram of the vector $`𝒙`$ with $`m`$ elements. The variables $`y`$ and $`x`$ may be multidimensional, and the multidimensional histograms can be mapped to $`n`$-bin ($`𝒙`$) and $`m`$-bin histograms ($`𝒚`$), respectively. The transition from $`𝒙`$ to $`𝒚`$ is described by the $`n`$-by-$`m`$ matrix $`𝑨`$. The element $`a_{ij}`$ is related to the probability to observe an entry in histogram bin $`i`$ of the histogram $`𝒚`$, if the true value $`x`$ is from histogram bin $`j`$ of the histogram $`𝒙`$. Problems with standard solutions are discussed in the next section.

In high energy physics experiments the problems is even more difficult than in other fields. Often the statistics of the measurement is not high and every $`𝒚`$-bin content (which is distributed due to the Poisson distribution around the expected value) has a large statistical fluctuation. Furthermore the resolution function $`𝑨(x,y)`$ (or the matrix $`𝑨`$) is not known analytically, but is represented by a data set from Monte Carlo simulation of the process, based on some assumed distribution $`f_{\text{MC}}(x)`$,

$$g_{\text{MC}}(y)=A(y,x)f_{\text{MC}}(x)\text{d}x,\text{(Monte Carlo simulation)}$$
(3)
and is also statistically limited. Standard methods for the solution of integral equations or linear equations can not be used in this case.

A simple method like the so-called *bin-by-bin correction* may be meaningful if the measurements $`y`$ are very close to the true values $`x`$. Real *unfolding* methods, taking all the correlations into account, are essential if there are larger effects of *transformation* and *finite resolution*. A solution $`𝒙`$ has to be found, with small deviations between the elements of $`𝑨`$$`𝒙`$ and the elements of the actually measured histogram $`\widehat{𝒚}`$. In the maximum likelihood method a function $`F(𝒙)`$ is constructed as the negative log of the Likelihood function, which describes the statistical relations between data and result:

$$F(𝒙)=\mathrm{log}L(𝒙,𝒚,𝑨)$$
(4)
and the minimum of $`F(𝒙)`$ is determined. Wildly fluctuating results $`𝒙`$ are due to large (negative) correlations between adjacent bins and are not acceptable. The approach to get a more reasonable solution is to impose a measure of the smoothness on the result $`𝒙`$; this method is called regularization. This technique was proposed independently by Phillips  and by Thikhonov . For a function $`f(x)`$ the integrated square of the second derivative

$$C(f)=\left(\frac{\text{d}^2f}{\text{d}x^2}\right)^2\text{d}x$$
(5)
is often used in the regularization which in the linearized version of the problem can be expressed by a quadratic form $`C(𝒙)=𝒙^T𝑪𝒙`$ with a positive-semidefinite matrix $`𝑪`$ (derivatives are replaced by finite differences). Equation (4) is then modified to the form

$$F(𝒙)=\mathrm{log}L(𝒙,𝒚,𝑨)+\tau C(𝒙)$$
(6)
with a factor $`\tau `$ called regularization parameter.

The result of the minimization of the modified function $`F(𝒙)`$ of equation (6) will show smaller fluctuations than the result obtained from equation (4) and may be more useful to compare the measurement with theoretical predictions. However it is clear that unavoidably the regularization introduces a bias. The magnitude of the bias depends on the value of regularization parameter $`\tau `$. A very large value would result in a *linear* function $`f(x)`$ or distribution $`𝒙`$, respectively. It is clear that the method requires an a-priori knowledge about a smooth behaviour of $`f(x)`$. The function $`f_{\text{MC}}(x)`$ used in the Monte Carlo simulation of equation (3) is often very close to the final result $`f(x)`$, i.e. the ratio is rather smooth. This suggests to express $`f(x)`$ in the form $`f(x)=f_{\text{MC}}(x)\times f_{\text{mult}}(x)`$ and to rewrite equation (1) in the form

$$g(y)=\left[A(y,x)f_{\text{MC}}(x)\right]f_{\text{mult}}(x)\text{d}x.$$
(7)
For the discretized form the function $`f_{\text{MC}}(x)`$ can be absorbed in a redefinition of matrix $`𝑨`$ and the vector $`𝒙`$ is interpreted as discretization of the hopefully *smooth* function $`f_{\text{mult}}(x)`$. With this redefinition the equation (2) can remain unchanged. The program RUN  for regularized unfolding is available since almost two decades and has been used in many experiments; early applications are  and . It is based on the reinterpretation of matrix $`𝑨`$ and $`𝒙`$, as described above, and includes a method for the determination of the regularization parameter $`\tau `$ based on the available degrees of freedom. In the method described later in this paper some details are treated differently.

## 2 UNFOLDING AS AN ILL-POSED PROBLEM

The problems inherent to unfolding are discussed in a simple special case, assuming a resolution matrix $`𝑨`$ with some smearing of data into neighbour bins. Assuming a true vector $`𝒙`$ the product $`𝒚=𝑨𝒙`$ describes the distribution expected due to the migration effect. With the same dimensions for the vectors $`𝒙`$ and $`𝒚`$ the matrix $`𝑨`$ is a square matrix and in the example later in this section the following simple symmetric form is assumed for the matrix $`𝑨`$, which depends on a single parameter $`\epsilon `$ ($`\epsilon `$ = migration parameter); for a 5-by-5 matrix the form is

$$𝑨=\left(\begin{array}{ccccc}1\epsilon & \epsilon & 0& 0& 0\\ \epsilon & 12\epsilon & \epsilon & 0& 0\\ 0& \epsilon & 12\epsilon & \epsilon & 0\\ 0& 0& \epsilon & 12\epsilon & \epsilon \\ 0& 0& 0& \epsilon & 1\epsilon \end{array}\right).$$
(8)
A direct solution for $`𝒙`$, given a measurement $`\widehat{𝒚}`$ differing from the expectation $`𝑨`$$`𝒙`$ with the true vector $`𝒙`$ by statistical fluctuations, is possible with inversion of the matrix $`𝑨`$:

$$\text{estimate}\widehat{𝒙}=𝑨^1\widehat{𝒚}\text{error propagation}𝑽(\widehat{𝒙})=𝑨^1𝑽_y\left(𝑨^1\right)^T.$$

The result has certain good statistical properties, for example it has no bias: $`E\left[\widehat{𝒙}\right]=𝑨^1E\left[𝒚\right]=𝑨^1𝑨E\left[𝒙\right]=𝒙`$. In practice the result is however satisfactory only for a matrix $`𝑨`$ with dominating diagonal; the result looks terrible if the matrix $`𝑨`$ describes a large migration to neighbour bins. Consequently the problem is called an ill-posed problem. In the following the solution of the equation $`𝒚=𝑨𝒙`$ using an orthogonal decomposition is discussed; this will allow some insight into the unfolding problem.

The symmetric matrix $`𝑨`$ is expressed by

$$𝑨=𝑼𝑫𝑼^T$$
(9)
with a transformation matrix $`𝑼`$ with property $`𝑼^T𝑼=\mathrm{𝟏}`$, and a diagonal matrix $`𝑫`$, where the diagonal elements of matrix $`𝑫`$ are the eigenvalues $`\lambda _j`$ of matrix $`𝑨`$ (in the order of decreasing value). The transformation matrix $`𝑼`$ contains the corresponding eigenvectors with the eigenvector $`𝒖_j`$ in the $`j`$-th column. The condition number $`\kappa `$ of a matrix is defined by the ratio of eigenvectors $`\kappa =\lambda _{\text{max}}/\lambda _{\text{min}}`$; the value of $`\kappa `$ is important for the quality of unfolding (see below). For values above $`\epsilon =0.20`$ the condition number $`\kappa `$ is very rapidly increasing.

A transformation of equation $`𝒚=𝑨𝒙`$ to a new basis is done by multiplication with matrix $`𝑼^T`$ (which is a rotation in the $`n`$-dimensional space):

$`𝑼^T|𝒚`$ $`=`$ $`𝑨𝒙=𝑼𝑫𝑼^T𝒙.`$
$`\underset{=𝒄}{\underset{}{𝑼^T𝒚}}`$ $`=`$ $`𝑫\underset{=𝒃}{\underset{}{𝑼^T𝒙}}𝒄=𝑫𝒃.`$
The matrix $`𝑼`$ of eigenvectors $`𝒖_j^T`$ allows to transform the vectors $`𝒙`$ and $`𝒚`$ to vectors $`𝒃=𝑼^T𝒙`$ and $`𝒄=𝑼^T𝒚`$, and to transform these vectors back by $`𝒙=𝑼𝒃`$ and $`𝒚=𝑼𝒄`$. The transformed equation $`𝒄=𝑫𝒃`$ with the diagonal matrix $`𝑫`$ shows, that each of the coefficients $`b_j`$ and $`c_j`$ is transformed independently of any other coefficient by the simple relation $`c_j=\lambda _jb_j`$. This operation does not depend on any assumption of the solution $`𝒙`$, and depends only on the properties of the matrix $`𝑨`$. Folding ($`𝒙𝒚`$) and unfolding ($`𝒚𝒙`$) is multiplication and division by the eigenvalues $`\lambda _j`$, respectively, of the coefficients in the transformed space.

In order to unfold a measured vector $`𝒚`$, the vector is transformed by $`𝒄=𝑼^T𝒚`$ to coefficients $`c_j`$, which have values influenced by statistical fluctuations of the elements of vector $`𝒚`$. In the unfolding the coefficients $`c_j`$ are divided by the eigenvalues $`\lambda _j`$ to obtain $`b_j=c_j/\lambda _j`$; the statistical fluctuation of coefficient $`c_j`$ is magnified for small eigenvalues $`\lambda _j`$ (i.e. $`\lambda _j1`$). Eventually, for very small eigenvalues $`\lambda _j`$, the final result $`𝒙=𝑼𝒃`$ will be dominated by one or by few of the coefficients $`b_j`$ with small eigenvalues and large statistical errors, and the complete result is unsatisfactory.

Example. In a numerical example the matrix $`𝑨`$ has the form of equation (8) with $`n=20`$ and a value of the migration parameter of $`\epsilon =0.22`$. The first eigenvalue is $`\lambda _1=1.0`$, and the last one is $`\lambda _{20}=1/7.9`$, giving a condition number $`\kappa =7.9`$. For $`𝒙`$ the ideal distribution of Figure 2a is assumed; the underlying function is of the form $`x\mathrm{exp}(ax)`$. The decomposition of the matrix $`𝑨`$ according to equation (9) is performed and the coefficients $`b_j`$ and $`c_j`$ are calculated. These coefficients are shown in Figure 3a (with $`b_jc_j`$). In addition this figure shows, calculated by standard error propagation, the almost constant error level of the coefficients, of the folded distribution of Figure 2a with Poisson distributed bin contents. Figure 3a shows, that the coefficients $`b_j`$ of the true distribution decrease rapidly with increasing value $`j`$ of the index of the coefficient, by roughly three orders of magnitude. The coefficients $`c_j`$ of the folded distribution drop even faster, because it is more smooth due to the migration effect. Of course the relation $`b_j/c_j=\lambda _j`$ is valid. The last coefficient $`b_j`$ in Figure 3a is reduced to $`c_j`$ by the inverse of the condition number of the matrix, which is $`\kappa =7.9`$ in this case.

The components of the first eigenvector $`u_1`$ (eigenvalue = 1) are all the same. Thus the coefficients $`b_i`$ and $`c_i`$ are identical, and proportional to the total sum of the measured distribution, not at all influenced by the migration. If visualized by functions, interpolating the components the eigenvector $`u_j`$ (eigenvalue $`\lambda _j`$) has $`j1`$ zeros, and the curvature of the visualized eigenvectors is rapidly increasing with index $`j`$. The components of the last eigenvector $`u_n`$ have alternating sign for the bins; it has a small absolute value and its measured value will have a large relative statistical error. The value of $`b_{20}`$ is obtained by $`b_{20}=7.9c_{20}`$ in unfolding, introducing a large bin-to-bin oscillation into the result of unfolding.

In a simulation Poisson distributed bin contents are assumed in the measurement vector $`𝒚`$. The coefficents for this measured distribution are shown in Figure 3b, together with the level of the statistical error. As expected from the size of the errors all coefficients with an index above about $`j=12`$ are dominated by the statistical error and therefore do not significantly contribute to the information content of the measurement. For indices above $`j=12`$ even the sign of the coeffient can not be determined by the measurement.

Using all the ”measured” coefficients for the unfolding the result of Figure 2b is obtained. This result shows large fluctuations around the expected values shown by the curve. The fluctuations are due to the contributions from indices above $`j=12`$, which represent noise and are magnified in the unfolding because of the large values of their inverse eigenvalues. The result is clearly unsatisfactory.

Because all measured coefficients $`c_j`$ with $`j`$ above a value of 12 are dominated by statistical errors (noise) their use in the unfolding makes no sense. A sharp cut-off after index $`j=12`$ or even after index $`j=10`$ will not remove any useful information from the measurement. The unfolding result using only measured coefficients $`c_j`$ up to $`j=10`$ is shown in Figure 2c; compared to Figure 2b the large fluctuations are suppressed and the results seems to be acceptable. Of course the fine structure of the true distribution expressed by the true coefficients $`b_j`$ with $`j>10`$ is not included in the solution and this may represent a bias. It is an unavoidable bias because these coefficients can not be measured.

The covariance matrix of the result can be calculated by standard error propagation. However it is clear that the covariance matrix is singular and has only rank 10 in this case, because the 20 bins are obtained from 10 measured coefficients (10 degrees of freedom). This property is inherent to the cut-off method and to the regularization method, and was already mentioned in . Such singularity of the covariance matrix can be avoided if the final transformation is to a number of bins identical to the degree of freedoms; only a limited number of bins can be obtained in a measurement with large miration effects.

This method of using a sharp cut-off has to be compared to the regularization method. It has been shown  that the use of a regularization function of the type of equation (5) is equivalent to a *smooth* cut-off; essentially the measured coefficients $`c_j`$ are multiplied by a factor depending on the curvature of the orthogonal contributions (see section 3).111 Sometimes the iterative solution of unfolding problems expressed by the equation $`𝒚=𝑨𝒙`$ is proposed in the literature without explicit regularization, starting from a ”good” initial distribution for $`𝒙`$. Of course equations of this type (with a square matrix) have a unique solution and iterative solutions are slow compared to the direct solution; after a large number of iterations with convergence the same unsatisfactory result as by direct solution will be obtained. However in these proposals only a small number of iterations is recommended. It can be shown that iterative methods can in fact include an implicit regularization : there is a different speed of convergence for the various orthogonal contributions and the small contributions with a small eigenvalue will converge very slowly. Thus after a few iterations the (large) coefficients with large eigenvalues are already accurate; the remaining coefficients are still almost unchanged and thus, for a stop after few iterations, their values are still close to the starting values. There is of course some subjectivity in stopping ”early” enough.

## 3 THE PROPOSED UNFOLDING METHOD

The proposed method is similar to the method used in RUN; the differences are emphasized in this section. It is expected that the proposed modifications results in more stable solutions. The proposed method requires large dimension parameters in the resolution matrix $`𝑨`$. Like in RUN the regularization is determined by the required number of degrees of freedom, which determines the regularization parameter.

Figures in this section refer to the example already mentioned in section 1 In a Monte Carlo calculation of all three effects, limited ($`x`$-dependent) acceptance, non-linear transformation and finite resolution are simulated. Details on the function and the distorting effects are identical to the published examples . In total 100 000 ”events” are simulated for ”data” and for the MC defining matrix $`𝑨`$. The input function $`f_{\text{MC}}(x)`$ (equation (7)) is a constant.

In RUN the discretization for $`f(x)`$ and for $`A(y,x)`$ was done using cubic B-spline functions; the effect is the same as for simple histograms namely the integral equation is transformed to a system of linear equations, however the elements of the vectors are B-spline coefficients instead of bin contents. The advantage is that the parametrized solution is a *smooth* function and the curvature as defined by equation (5) can be exactly written as a quadratic form. However the accurate determination of matrix $`𝑨`$ requires a good Monte Carlo statistic. In RUN statistical fluctuations of the elements of matrix $`𝑨`$ could not be treated.

Simple histograms are instead proposed here; the elements of the vector $`𝒚`$ are bin contents (integer numbers). The curvature of the solution is constructed by finite differences: the second derivative in bin $`j`$ is proportional to $`x_{j1}2x_j+x_{j+1}`$. In a histogram some resolution is lost if bins with a width as large as expected for the final resolution would be used. It is recommended to use initially $`m=2n_{\text{df}}`$ bins for $`𝒙`$ for a final number of degrees of freedom of $`n_{\text{df}}`$. For $`y`$ a larger number of bins $`n`$ ($`>m`$) is recommended, to avoid a loss of resolution. Thus the number $`n\times m`$ of elements is large, and a large sample of Monte Carlo events is required to *fill* matrix $`𝑨`$. The statistical error of the elements $`a_{ij}`$ eventually can not be neglected.

Standard Poisson maximum likelihood fit. Ignoring initially eventual statistical errors of the elements $`a_{ij}`$ the expected number of events in bin $`i`$ of $`𝒚`$ is given by $`y_i=_{j=1}^ma_{ij}x_j`$ . For the expected number $`y_i`$, as given by this expression, the observed values $`\widehat{y}_i`$ follows the Poisson distribution. Optimal estimates for the elements $`x_j`$ are obtained by minimizing the (negative) logarithm of the total likelihood with respect to the elements $`x_j`$ of vector $`𝒙`$, assuming the Poisson distribution:

$$F(𝒙)=\mathrm{ln}(𝒙)=\mathrm{ln}\left[\underset{i=1}{\overset{n}{}}P_{y_i}(\widehat{y}_i)\right]=\underset{i=1}{\overset{n}{}}\left(y_i\widehat{y}_i\mathrm{ln}y_i\right)+\text{const.},$$
(10)
where the constant term containing e.g $`\widehat{y}_i!`$ can be ommited. This expression (10) correctly accounts also for bins with a small number of histogram entries $`\widehat{y}_i`$.

An alternative would be to use the (linear) least squares method with singular value decomposition for the fit. However for small number of entries the use of the Poisson distribution seems to be essential. Furthermore the diagonalization used later in the method is almost equivalent to singular value decomposition (eigenvalues are the squares of the singular values).

Fitting with finite Monte Carlo samples. The problem of statistical fluctuations of the elements $`a_{ij}`$ has been neglected so far. A method to treat the problem within the maximum-likelihood method has been developed by R.Barlow and Chr.Beeston . In this method there is for each source bin $`x_j`$ some (unknown) expected number of events $`A_{ij}`$. For each element $`A_{ij}`$ the corresponding number $`a_{ij}`$ from the Monte Carlo sample is generated by a distribution which is taken to be Poisson too. The nice feature of this method is that a bias which would be introduced by ignoring the statistical character of the values of the elements $`a_{ij}`$ is avoided and the maximum likelihood error is more realistic. A large number of slack variables (one for each bin) is introduced and has to be treated in the optimzation. A new fast numerical solution method has been developed (see ).

Combining bins. The likelihood function is a sum over all bins. Combining almost empty bins does not introduce a systematic error. The total number of elements of the matrix may be large, especially if $`x`$ and/or $`y`$ are multidimensional, and a small number of entries (or even zero) in an element may not be uncommon. The combination of almost empty bins is done with a cluster algorithm, taking into account the distance between bins in one, two or three dimensions.

First option: sharp cut-off of orthogonal contributions. This method is rather similar to the method discussed in section 2. The computational problem is to determine the minimum of $`F(𝒙)`$ (see equation (10)). The standard iterative method is based on the representation for the correction $`\mathrm{\Delta }𝒙`$

$$F(\mathrm{\Delta }𝒙)=\frac{1}{2}\mathrm{\Delta }𝒙^T𝑯\mathrm{\Delta }𝒙+\mathrm{\Delta }𝒙^T𝒈+\mathrm{}$$
(11)
with the Hessian $`𝑯`$ (matrix of second derivatives of $`F(\mathrm{\Delta }𝒙)`$) and the gradient vector $`𝒈`$ (first derivatives of $`F(\mathrm{\Delta }𝒙)`$). A Newton step is then calulated from equation $`𝑯\mathrm{\Delta }𝒙+𝒈=0`$. Convergence is usually fast for good starting values and the covariance matrix is equal to the inverse of the Hessian. The starting values can be calculated by a linear least square fit, based on the approximation of the Poisson distribution by a Gaussian distribution for each bin.

A sharp cut-off as discussed in the example of section 2 requires a diagonalization of the symmetric matrix $`𝑯`$ by $`𝑯=𝑼𝑫𝑼^T`$ with a diagonal matrix $`𝑫`$ and a transformation matrix $`𝑼`$. By a transformation (rotation) in $`𝒙`$-space linear combinations of the $`𝒙`$-components are obtained with a diagonal covariance matrix, with variances of the linear combinations given by the inverse of the eigenvalues of matrix $`𝑫`$. A cut-off is done at some index $`j`$ followed by backtransformation to the $`𝒙`$-space of bin-contents using the transformation matrix $`𝑼`$.

Second option: regularization. In this option the regularization is based on the second derivative of the result according to equation (5), which can be expressed by a quadratic form $`𝒙^T𝑪𝒙`$ with a positive-semidefinite matrix $`𝑪`$. In principle the same procedure is used as in RUN; the mathematical details are given elsewhere . Here a simple explanation is given on the standard mathematical operations222In a publication the method has been described to ”*have certain mathematical complications*”, but it is based only on standard linear algebra of symmetric matrices. used. Regularization is done by adding the term $`\tau 𝒙^T𝑪𝒙`$ to the function $`F(\mathrm{\Delta }𝒙)`$ of equation (11). Exactly as in the first option the Hessian is diagonalized.

$$𝑯=𝑼𝑫𝑼^T𝑯^1=𝑼𝑫^1𝑼^T=\left[𝑼𝑫^{1/2}\right]\left[𝑫^{1/2}𝑼^T\right].$$
(12)
Up to this step everything is identical to the cut-off option. Using transformation matrix $`𝑼𝑫^{1/2}`$ the vector $`𝒙`$ is transformed to linear combinations $`\widehat{𝒙}`$, which are orthogonal, with all variances equal to 1 (unit covariance matrix). Because the covariance matrix is equal to the unit matrix, every additional pure rotation will not change the (unit) covariance matrix. In terms of the transformed vector the regularization term can now be written in the form $`\tau \widehat{𝒙}^TC_U\widehat{𝒙}`$, where $`C_U`$ is the transformed curvature matrix $`𝑪`$. Now another diagonalization can be done of matrix $`C_U`$:

$$\tau 𝒙^T𝑪𝒙\tau \widehat{𝒙}^T𝑪_U\widehat{𝒙}=\tau \widehat{𝒙}^T𝑼_C𝑺𝑼_C^T\widehat{𝒙}$$
(13)
with a diagonal matrix $`𝑺`$ and a rotation matrix $`𝑼_C`$. This diagonalization can be used to define a pure rotation from the linear combination $`\widehat{𝒙}`$ to another linear combination $`\stackrel{~}{𝒙}`$

$$\widehat{𝒙}\stackrel{~}{𝒙}=𝑼_C^T\widehat{𝒙}.$$
(14)
The components of the new vector $`\stackrel{~}{𝒙}`$ still have the unit matrix as covariance matrix. The complete transformation from $`𝒙`$ to $`\stackrel{~}{𝒙}`$ is the effect of the transformation by $`𝑼𝑫^{1/2}`$ and by $`𝑼_C`$. The algebra can be explained in other words: the error ellipsoid related to the Hessian is first rotated to have the axes parallel to the axes of the new system. By a change of the scales the ellipsoid is transformed to a sphere, which will remain a sphere for any further rotation. A last rotation is done to bring the axes into the order of increasing curvature.

Some columns of the complete (product) transformation are shown in Figure 4. All linear combinations obtained have the same precision (standard deviation of the coefficient is one). As seen in the Figure linear combinations with large index $`j`$ are oscillating with large amplitude. The diagonal elements $`S_{jj}`$ are the (statistically independent) contributions of the elements of $`\stackrel{~}{𝒙}`$ to the total curvature. Sorted according to increasing value of $`S_{jj}`$ the value of $`S_{jj}`$ will increase rather fast with increasing index $`j`$. The spectrum of eigenvalues $`S_{jj}`$ is shown in Figure 5. In terms of the linear combinations $`\stackrel{~}{𝒙}`$ regularization is simply given by

$$\left(\stackrel{~}{x}_j\right)_{\text{reg}}=\left(\frac{1}{1+\tau S_{jj}}\right)\left(\stackrel{~}{x}_j\right)_{\text{unreg}}.$$
(15)
and this simple form is the reason for the transformations made before.

Determination of the regularization parameters $`\tau `$. The first factors (small $`j`$) on the right-hand-side of equation (15) will be close to 1; for a value $`\tau =1/S_{kk}`$ the factor will be 1/2 and for indices $`j>k`$ will rapidly decrease towards zero. The sum of all factors can be called the *effective* number of degrees of freedom, and can be used to determine the value of the regularization parameter $`\tau `$ from the required number of degrees of freedom, i.e. the regularization parameter $`\tau `$ is determined from the value of $`n_{\text{df}}`$ in the equation

$$n_{\text{df}}=\underset{j=1}{\overset{m}{}}\left(\frac{1}{1+\tau S_{jj}}\right).$$
(16)
Thus the required number of degrees of freedom has to be specified and determines the degree of regularization. This number can be taken from the spectrum of the coefficients or amplitudes, shown in Figure 5. The insignificant part (large $`j`$) is clearly visible in the spectrum and separated from the significant part (small $`j`$). The selected value of $`n_{df}`$ should be equal to or larger than the number of significant terms. The unregularized amplitudes, which have standard deviation one, are shown by the left bars; amplitudes above index 15 are compatible with one and represent noise. They would however make a large contribution to the solution, because the corresponding column vectors (Figure 4) are large. The regularization effectively damps the amplitude (right bars) around and above index 15, which has been chosen as the degree of freedom here. The significant amplitudes are not affected by the regularization.

The final result of the example (measured distribution in Figure 1) is shown in Figure 6. The left figure shows 30 data points with error bars together with the original (true) distribution; within errors the original distribution is nicely reproduced. The rank of the covariance matrix is about 15, which was chosen as the effective number of degrees of freedom; thus inversion of the covariance matrix, needed e.g. for a least-square fit of a model to the data, is not possible. Although the large number of 30 data points seems to be attractive, the data points should be reduced to 15 data points by combining two bins to one, which then have a full-rank covariance matrix. This set of data points is shown in Figure 6 (right). The broader bins of this set of data points are a consequence of the limited acceptance and finite resolution of the measurement.

## ACKNOWLEDGEMENTS

I would like to thank the organizers of the conference on Advanced Statistical Techniques in Particle Physics for their hospitality and the stimulating atmosphere in Durham.


text:
# 1 Introduction

## 1 Introduction

Recently there have been important new data on the charm structure function (SF) $`F_2^c`$, of the proton from the H1  and ZEUS  Collaborations at HERA, which have probed the small-$`x`$ region down to $`x=8\times 10^4`$ and $`x=2\times 10^4`$, respectively. At these values of $`x`$, the charm contribution to the total proton SF, $`F_2^p`$, is found to be around $`25\%`$, which is a considerably larger fraction than that found by the European Muon Collaboration at CERN  at larger $`x`$, where it was only $`1\%`$ of $`F_2^p`$. Extensive theoretical analyses in recent years have generally served to confirm that the $`F_2^c`$ data can be described through perturbative generation of charm within QCD (see, for example, the review in Ref.  and references therein).

We note, that perhaps more relevant analyses of the HERA data, where the $`x`$ values are quite small, are those based on BFKL dynamics , because the leading $`ln(1/x)`$ contributions are summed. The basic dynamical quantity in BFKL approach is the unintegrated gluon distribution $`\phi _g(x,k_{}^2)`$ 111Hereafter $`p^\mu `$ and $`k^\mu `$ are the hadron and the gluon 4-momentums, respectively, and $`q^\mu `$ is the photon 4-momentum. ($`f_g`$ is the (integrated) gluon distribution multiplied by $`x`$ and $`k_{}`$ is the transverse momentum)

$`f_g(x,Q^2)=f_g(x,Q_0^2)+{\displaystyle _{Q_0^2}^{Q^2}}𝑑k_{}^2\phi _g(x,k_{}^2),`$ (1)
which satisfies the BFKL equation. The integral is divergent at the lower limit and it leads to the necessity to use the difference $`f_g(x,Q^2)f_g(x,Q_0^2)`$ with some nonzero $`Q_0^2`$.

Then, in the BFKL-like approach (hereafter the $`k_t`$-factorization approach  is used) the SFs $`F_{2,L}^c(x,Q^2)`$ are driven at small $`x`$ by gluons and are related in the following way to the unintegrated distribution $`\phi _g(x,k_{}^2)`$:

$`F_{2,L}^c(x,Q^2)={\displaystyle _x^1}{\displaystyle \frac{dz}{z}}{\displaystyle 𝑑k_{}^2C_{2,L}^g(z,Q^2,m_c^2,k_{}^2)\phi _g(x/z,k_{}^2)},`$ (2)
The functions $`C_{2,L}^g(x,Q^2,m_c^2,k_{}^2)`$ may be regarded as the structure functions of the off-shell gluons with virtuality $`k_{}^2`$ (hereafter we call them as hard structure functions222This notation reflects the fact that structure functions $`F_{2,L}^c`$ connect with the functions $`C_{2,L}^g`$ at the same form as cross-sections connect with hard ones (see ).). They are described by the quark box (and crossed box) diagram contribution to the photon-gluon interaction (see Fig. 1 in ).

The purpose of the paper is to present the results for these hard SF $`C_{2,L}^g(x,Q^2,m_c^2,k_{}^2)`$ and to analyze experimental data for $`F_2^c(x,Q^2)`$ by applying Eq. (2) with different sets of unintegrated gluon densities (see Ref. ) and to give predictions for the longitudinal SF $`F_L^c(x_B,Q^2)`$.

It is instructive to note that the results should be similar to those of the photon-photon scattering process. The corresponding QED contributions have been calculated many years ago in Ref.  (see also the beautiful review in Ref. ). Our results have been calculated independently in  (based on approaches of ) and they are in full agreement with . However, we hope that our formulas which are given in a simpler form could be useful for others.

## 2 Hard structure functions

The gluon polarization tensor (hereafter the indices $`\alpha `$ and $`\beta `$ are connected with gluons), which gives the main contribution at high energy limit, has the form:

$`\widehat{P}_{BFKL}^{\alpha \beta }={\displaystyle \frac{k_{}^\alpha k_{}^\beta }{k_{}^2}}={\displaystyle \frac{x^2}{k^2}}p^\alpha p^\beta ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{\stackrel{~}{\beta }^4}}\left[\stackrel{~}{\beta }^2g^{\alpha \beta }12bx^2{\displaystyle \frac{q^\alpha q^\beta }{Q^2}}\right],`$ (3)
where $`\stackrel{~}{\beta }^2=14bx^2,b=k^2/Q^2k_{}^2/Q^2>0.`$

Contracting the corresponding photon projectors, we have:

$`C_2^g(x)={\displaystyle \frac{𝒦}{\stackrel{~}{\beta }^2}}\left[f_{BFKL}^{(1)}+{\displaystyle \frac{3}{2\stackrel{~}{\beta }^2}}f_{BFKL}^{(2)}\right],C_L^g(x)={\displaystyle \frac{𝒦}{\stackrel{~}{\beta }^2}}\left[4bx^2f_{BFKL}^{(1)}+{\displaystyle \frac{(1+2bx^2)}{\stackrel{~}{\beta }^2}}f_{BFKL}^{(2)}\right]`$ (4)
$`f_{BFKL}^{(1)}={\displaystyle \frac{1}{\stackrel{~}{\beta }^4}}\left[\stackrel{~}{\beta }^2f^{(1)}3bx^2\stackrel{~}{f}^{(1)}\right],f_{BFKL}^{(2)}={\displaystyle \frac{1}{\stackrel{~}{\beta }^4}}\left[\stackrel{~}{\beta }^2f^{(2)}3bx^2\stackrel{~}{f}^{(2)}\right],`$
where the normalization factor $`𝒦=e_c^2\alpha _s(Q^2)/(4\pi )x`$, $`a=m_c^2/Q^2`$ , $`e_c=2/3`$ is charm quark charge and

$`f^{(1)}`$ $`=`$ $`2\beta [1(12x(1+b2a)[1x(1+b+2a)])f_1`$ (5)
$`+`$ $`(2ab)(12a)x^2f_2],`$
$`f^{(2)}`$ $`=`$ $`8x\beta [(1(1+b)x)2x(bx(1(1+b)x)(1+b2a)+a\stackrel{~}{\beta }^2)f_1`$ (6)
$`+`$ $`bx^2(1(1+b)x)(2ab)f_2],`$
$`\stackrel{~}{f}^{(1)}`$ $`=`$ $`\beta [{\displaystyle \frac{1x(1+b)}{x}}2(x(1x(1+b))(1+b2a)+a\stackrel{~}{\beta }^2)f_1`$ (7)
$``$ $`x(1x(1+b))(12a)f_2],`$
$`\stackrel{~}{f}^{(2)}`$ $`=`$ $`4\beta (1(1+b)x)^2\left[2(1+2bx^2)f_1bx^2f_2\right],`$ (8)
with

$$\beta ^2=1\frac{4ax}{(1(1+b)x)},f_1=\frac{1}{\stackrel{~}{\beta }\beta }\mathrm{ln}\frac{1+\beta \stackrel{~}{\beta }}{1\beta \stackrel{~}{\beta }},f_2=\frac{4}{1\beta ^2\stackrel{~}{\beta }^2}$$

For the important regimes when $`k^2=0`$, $`m_c^2=0`$ and $`Q^2=0`$, the results are coincided with ones of Refs. . Notice that our results in Eq. (4) should be also agree with those in Ref.  but the direct comparison is quite difficult because the structure of their results is quite cumbersome (see Appendix A in Ref. ). We have found numerical agreement in the case of $`F_2(x,Q^2)`$ for several types of unintegrated gluon distributions (see Fig.4 in ).

## 3 Comparison with $`F_2^c`$ experimental data and predictions for $`F_L^c`$

With the help of the results obtained in the previous Section we have analyzed HERA data for SF $`F_2^c`$ from ZEUS  collaboration.

Notice that in Ref.  the $`k_{}^2`$-integral in the r.h.s. of Eq.(2) has been evaluated using the BFKL results for the Mellin transform of the unintegrated gluon distribution and the Wilson coefficient functions have been calculated for the full perturbative series at asymptotically small $`x`$ values. Since we would like to analyze experimental data for the SF $`F_2^c`$, we have an interest to obtain results at quite broad range of small $`x`$ values. For the reason we need in a parameterization of unintegrated gluon distribution function.

We consider two different parametrizations for the unintegrated gluon distribution (see ): the Ryskin-Shabelski (RS) one  and the Blumlein one .

In Fig. 1 we show the SF $`F_2^c`$ as a function $`x`$ for different values of $`Q^2`$ in comparison with ZEUS  experimental data. We see that at large $`Q^2`$ ($`Q^210`$ GeV2) the SF $`F_2^c`$ obtained in the $`k_T`$ factorization approach is higher than the SF obtained in the standard parton model with the GRV  gluon density at the LO approximation (see curve 1) and has a more rapid growth in comparison with the standard parton model results, especially at $`Q^2130`$ GeV2 . At $`Q^210`$ GeV2 the predictions from perturbative QCD (in GRV approach) and those based on the $`k_T`$ factorization approach are very similar 333This fact is due to the quite large value of $`Q_0^2=4`$ GeV2 chosen here. and show a disagreement with data below $`Q^2=7`$ GeV2 444A similar disagreement with data at $`Q^22`$ GeV2 has been observed for the complete structure function $`F_2`$ (see, for example, the discussion in Ref.  and reference therein). We note that the insertion of higher-twist corrections in the framework of usual perturbative QCD improves the agreement with data (see Ref. ) at quite low values of $`Q^2`$.. Unfortunately the available experimental data do not permit yet to distinguish the $`k_T`$ factorization effects from those due to boundary conditions .

The results for the SF $`F_L^c`$ obtained in perturbative QCD and from the $`k_T`$ factorization approach are quite similar to the $`F_2^c`$ case discussed above. The ratio $`R^c=F_L^c/F_2^c`$ is shown in Fig. 2. We see that $`R^c0.1÷0.3`$ in a wide region of $`Q^2`$. The estimation of $`R^c`$ is very close to the results for $`R=F_L/(F_2F_L)`$ ratio (see Refs. -). We would like to note that these values of $`R^c`$ contradict the estimation obtained in Refs. . The effect of $`R^c`$ on the corresponding differential cross-section should be considered in the extraction of $`F_2^c`$ from future more precise measurements.

For the ratio $`R^c`$ we found quite flat $`x`$-behavior at low $`x`$ in the low $`Q^2`$ region (see Fig. 2), where approaches based on perturbative QCD and on $`k_T`$ factorization give similar predictions (see Fig 1). It is in agreement with the corresponding behaviour of the ratio $`R=F_L/(F_2F_L)`$ (see Ref. ) at quite large values of $`\mathrm{\Delta }_P`$ 555The behaviour is in agreement with previous studies . Note that at small values of $`\mathrm{\Delta }_P`$, i.e. when $`x^{\mathrm{\Delta }_P}Const`$, the ratio $`R`$ has got strong negative NLO corrections (see Refs. ) and tends to zero at $`x0`$ after some resummation done in . ($`\mathrm{\Delta }_P>0.20.3`$). The low $`x`$ rise of $`R^c`$ at high $`Q^2`$ disagrees with early calculations  in the framework of perturbative QCD. It could be due to the small $`x`$ resummation, which is important at high $`Q^2`$ (see Fig 1). We plan to study this effect in future.

## 4 Conclusions

We have presented the results for the perturbative parts of the SFs $`F_2^c`$ and $`F_L^c`$ for a gluon target having nonzero momentum squared, in the process of photon-gluon fusion.

We have applied the results in the framework of $`k_T`$ factorization approach to the analysis of present data for the charm contribution to $`F_2`$ and we have given the predictions for $`F_L^c`$. The analysis has been performed with two parameterizations of unintegrated gluon distributions for comparison. We have found good agreement of our results with experimental HERA data for $`F_2^c`$, except at low $`Q^2`$ ($`Q^27`$ GeV2) 666It must be noted that the cross section of inelastic $`c\overline{c}`$ and $`b\overline{b}`$pair photoproduction at HERA are described by the Blumlein parametrization at a smaller value of $`Q_0^2`$ ($`Q_0^2=1`$ GeV2) (see ).. We have also obtained quite large contribution of the SF $`F_L^c`$ at low $`x`$ and high $`Q^2`$ ($`Q^230`$ GeV2).

We would like to note the good agreement between our results for $`F_2^c`$ and the ones obtained in Ref.  by Monte-Carlo studies. Moreover, we have also good agreement with fits of H1 and ZEUS data for $`F_2^c`$ (see recent reviews in Ref.  and references therein) based on perturbative QCD calculations. But unlike to these fits, our analysis uses universal unintegrated gluon distribution, which gives in the simplest way the main contribution to the cross-section in the high-energy limit.

It could be also very useful to evaluate the complete $`F_2`$ itself and the derivatives of $`F_2`$ with respect to the logarithms of $`1/x`$ and $`Q^2`$ with our expressions using the unintegrated gluons.

The consideration of the SF $`F_2`$ in the framework of the leading-twist approximation of perturbative QCD (i.e. for “pure” perturbative QCD) leads to very good agreement (see Ref.  and references therein) with HERA data at low $`x`$ and $`Q^21.5`$ GeV2. The agreement improves at lower $`Q^2`$ when higher twist terms are taken into account . As it has been studied in Refs. , the SF $`F_2`$ at low $`Q^2`$ is sensitive to the small-$`x`$ behavior of quark distributions. Thus, our future analysis of $`F_2`$ in a broader $`Q^2`$ range in the framework of $`k_T`$ factorization should require the incorporation of parametrizations for unintegrated quark densities, introduced recently (see Ref.  and references therein).

The SF $`F_L`$ is very interested too (see ). The structure function $`F_L`$ depends strongly on the gluon distribution (see, for example, Ref. ), which in turn is determined  by the derivative $`dF_2/d\mathrm{ln}Q^2`$. Thus, in the framework of perturbative QCD at low $`x`$ the relation between $`F_L`$, $`F_2`$ and $`dF_2/d\mathrm{ln}Q^2`$ could be violated by non-perturbative contributions, which were expected to be important in the $`F_L`$ case (see Refs. ). The application in  of present analysis to $`F_L`$ gives a “non pure” perturbative QCD predictions for the structure function which are in good agreement with data  and with the “pure” perturbative results of Ref. , that determines a small value of nonperturbative corrections to SF $`F_L`$ at low $`x`$ and quite large $`Q^2`$ values.

Acknowledgments. Authors would like to express their sincerely thanks to the Organizing Committee of the XVIth International Workshop “High Energy Physics and Quantum Field Theory” for the kind invitation, the financial support at such remarkable Conferences, and for fruitful discussions.

A.V.K. was supported in part by Alexander von Humboldt fellowship and INTAS grant N366. G.P. acknowledges the support of Galician research funds (PGIDT00 PX20615PR) and Spanish CICYT (FPA2002-01161). N.P.Z. also acknowledge the support of Royal Swedish Academy of Sciences.