format: text

text:
# Detected jump-error correcting quantum codes, quantum error designs and quantum computation

## I Introduction

The discovery of powerful, fast quantum algorithms Shor1 ; Grover ; Simon  launched new efforts to implement such quantum algorithms in real physical systems. Quantum algorithms simultaneously exploit two characteristic features of quantum theory. Namely, the fundamental phenomenon of quantum interference and the fact that for distinguishable quantum systems the dimension of the Hilbert space increases exponentially with the number of systems. Therefore, to implement a quantum algorithm in a real quantum system we must be able to create and manipulate arbitrary superpositions of quantum states and to preserve quantum coherence during a computation. Unfortunately, quantum coherence is very fragile. Typically, any coupling to an environment leads to decoherence so that quantum mechanical superpositions are rapidly destroyed.

The urgent need to develop efficient methods to protect quantum coherence has led to the study of very general classes of quantum error-correcting codes Shor2 ; Knill1 ; Ekert ; Knill2 . The main idea is to restrict the dynamics of a quantum algorithm to a subspace of the Hilbert space, in which errors can be identified uniquely by suitable measurements and where the error operators can be inverted by unitary operations. Typically, this is achieved by an encoding of the logical information and by a suitable choice of quantum gates.

For some special cases it is also possible to design a passive error-correcting quantum code Lidar1 ; Zanardi ; Lidar2 . Such a passive quantum code relies on a subspace of the Hilbert space which is not affected by any errors at all. In this situation the unitary recovery operation is the identity operation so that an active correction of the errors is not necessary.

In principle, any type of error can be corrected by these strategies as long as enough physical qubits are available to achieve the required redundancy and one can make a large number of control measurements and perform the rapid recovery operations. However, in view of current-day experimental possibilities Sackett  it is generally difficult to achieve both requirements. Therefore it is desirable to develop alternative error-correcting strategies which possibly correct a restricted class of errors only, but which tend to minimize both redundancy and the number of recovery operations. Recently, the first steps in this direction have been taken by defining a new class of one detected jump-error correcting quantum codes which are capable of stabilizing distinguishable qubits against spontaneous decay processes into statistically independent reservoirs Alber01 . These codes are constructed by embedding an active error-correcting code in a passive code space and by exploiting information available on error positions. This embedding procedure leads to a significant reduction of redundancy and the number of control measurements and recovery operations.

In this paper the physical principles underlying detected jump-error correcting quantum codes are explored and generalized, motivated by the practical need for quantum error-correcting codes which minimize both redundancy and the number of recovery operations. Based on these physical principles an upper bound is established on the number of logical states of a general embedded detected jump-error correcting quantum code.

From this bound it is apparent that the recently discovered one detected jump-error correcting quantum codes have minimal redundancy. Based on this family of optimal one detected jump-error correcting quantum codes, we establish links with the general notions of combinatorial design theory Beth1 . For this purpose the new concept of a spontaneous emission error design is introduced. This is a powerful tool for constructing multiple detected jump-error correcting quantum codes capable of stabilizing distinguishable qubits against spontaneous decay processes. As an example, we present a new embedded three detected jump-error correcting quantum code.

This paper is organized as follows. In Sec. II basic physical aspects concerning the spontaneous emission of photons by qubit-systems are summarized. In Sec. III the physical principles are explored which lead to the construction of one detected jump-error correcting quantum codes. The conditions for general detected jump-error correcting quantum codes are given in Sec. IV. The links with combinatorial design theory are established in Sec. V. Finally, in Sec. VI numerical examples are presented which exhibit basic stability properties of the optimal one detected jump-error correcting quantum codes.

## II Spontaneous decay of distinguishable qubits and photon detection

In this section we summarize basic facts about the dynamical description of a quantum system interacting with initially unoccupied modes of the electromagnetic field. These considerations are the starting point for the development of optimal strategies of error correction, which we pursue in the subsequent sections.

We consider a model of a quantum computer in which $`N`$ two-level atoms (qubits) interact with external laser pulses which synthesize the quantum gates underlying a quantum algorithm. These $`N`$ qubits are assumed to be arranged in an array with well defined positions $`𝐱_\alpha `$ $`(\alpha =1,\mathrm{},N)`$ (see Fig. 1). In addition, these $`N`$ qubits are assumed to be distinguishable, which requires that their mean nearest neighbor distance is large in comparison with the optical wave lengths involved. Their distinguishability guarantees that the dimension of their associated Hilbert space is $`K_{}=2^N`$ and thus scales exponentially with the number of qubits. In addition, it is assumed that these $`N`$ qubits couple to the unoccupied modes of the electromagnetic field. This coupling causes spontaneous decay processes of the qubits from their excited states $`|1_\alpha `$, to their stable lower lying states $`|0_\alpha `$. Within the Born, Markov, and the rotating wave approximations the resulting dynamics of the reduced density operator $`\rho (t)`$ of this $`N`$-qubit system are described by the master equation Mollow ; Carmichael

$$\frac{d\rho (t)}{dt}=\frac{i}{\mathrm{}}\left(H_{\mathrm{eff}}\rho (t)\rho (t)H_{\mathrm{eff}}^{}\right)+\underset{\alpha =1}{\overset{N}{}}L_\alpha \rho (t)L_\alpha ^{}$$
(1)
with the non-hermitian effective Hamiltonian

$$H_{\mathrm{eff}}=H\frac{i\mathrm{}}{2}\underset{\alpha =1}{\overset{N}{}}L_\alpha ^{}L_\alpha .$$
(2)
Thereby, the coherent dynamics of the $`N`$-qubit system in the absence of the coupling to the vacuum modes of the electromagnetic field are described by the Hamiltonian $`H`$ which incorporates the influence of the external laser pulses. In addition, we assume that the mean distance between the qubits is much larger than the wave lengths of the spontaneously emitted radiation. Therefore, to a good degree of approximation each qubit $`\alpha `$ couples to a different set of modes of the radiation field so that these sets constitute statistically independent reservoirs Carmichael . In Eq. (1) the coupling of qubit $`\alpha `$ to its reservoir and the resulting spontaneous decay process $`|1_\alpha |0_\alpha `$ is characterized by the Lindblad operator

$$L_\alpha =\sqrt{\kappa _\alpha }𝟙_{\beta \alpha }|\mathrm{𝟘}_{\alpha \alpha }\mathrm{𝟙}|$$
(3)
where $`𝟙_{\beta \alpha }`$ denotes the identity on every except the $`\alpha `$-th qubit, and $`\kappa _\alpha `$ is the associated spontaneous decay rate.

Provided that initially the $`N`$-qubit system is in a pure state, say $`|\psi _0`$, a formal solution of the master Eq. (1) is given in Mollow

$`\rho (t)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle _0^t}𝑑t_n{\displaystyle _0^{t_n}}𝑑t_{n1}\mathrm{}{\displaystyle _0^{t_2}}𝑑t_1p(t;t_n\alpha _n,\mathrm{},t_1\alpha _1)`$ (4)
$`|t;t_n\alpha _n,\mathrm{},t_1\alpha _1t;t_n\alpha _n,\mathrm{},t_1\alpha _1|`$
with the pure quantum state

$`|t;t_n\alpha _n,\mathrm{},t_1\alpha _1`$
$`={\displaystyle \frac{e^{iH_{\mathrm{eff}}(tt_n)/\mathrm{}}L_{\alpha _n}\mathrm{}L_{\alpha _1}e^{iH_{\mathrm{eff}}t_1/\mathrm{}}}{p(t|t_n\alpha _n,\mathrm{},t_1\alpha _1)^{1/2}}}|\psi _0`$
and with the probabilities

$`p(t;t_n\alpha _n,\mathrm{},t_1\alpha _1)`$ (5)
$`=\psi _0|e^{iH_{\mathrm{eff}}^{}t_1/\mathrm{}}L_{\alpha _1}^{}\mathrm{}L_{\alpha _n}^{}e^{iH_{\mathrm{eff}}^{}(tt_n)/\mathrm{}}`$
$`e^{iH_{\mathrm{eff}}(tt_n)/\mathrm{}}L_{\alpha _n}\mathrm{}L_{\alpha _1}e^{iH_{\mathrm{eff}}t_1/\mathrm{}}|\psi _0.`$
It can be shown that each pure state $`|t;t_n\alpha _n,\mathrm{},t_1\alpha _1`$ describes the quantum state of the $`N`$-qubit system at time $`t`$ conditioned on the emission of precisely $`n`$ photons at times $`t_1\mathrm{}t_n`$ by qubits $`\alpha _1,\mathrm{},\alpha _n`$ Mollow ; Carmichael ; Carmichael2 ; Zoller . Thus, each of the pure quantum states of Eq. (4) corresponds to a possible measurement record $`(t_1\alpha _1,\mathrm{},t_n\alpha _n)`$ in an experiment in which each qubit is observed continuously by photodetectors. In the subsequent discussion it is important to note that due to the large separation between the qubits ideally this measurement record not only determines the spontaneous decay times $`t_i`$, but also the associated positions $`\alpha _i`$ ($`i=1,\mathrm{},n`$) of the qubits which have been affected by these decay processes. The measurement record $`(t_1\alpha _1,\mathrm{},t_n\alpha _n)`$ is observed with probability $`p(t;t_n\alpha _n,\mathrm{},t_1\alpha _1)`$. According to Eq. (4) the quantum state $`|t;t_n\alpha _n,\mathrm{},t_1\alpha _1`$ resulting from a particular measurement record is determined by two types of effects. First, the time evolution between two successive photon emission events is characterized by the non-hermitian Hamiltonian $`H_{\mathrm{eff}}`$ of Eq. (2). Thus, even in the absence of any spontaneous photon emission process in a given time interval $`[0,t]`$, in general an arbitrary quantum state is modified by the couplings to the vacuum modes of the electromagnetic field. Second, immediately after the spontaneous emission of a photon by qubit $`\alpha `$ the quantum state is modified by the Lindblad operator $`L_\alpha `$. This distinction between the two different types of dynamics is crucial for the development of a useful quantum error-correcting strategy which minimizes the redundancy as well as the number of control measurements and recovery operations.

## III Designing optimal one detected jump-error correcting quantum codes

In this section we introduce the main ideas leading to the construction of optimal one detected jump-error correcting embedded quantum codes which are capable of stabilizing systems of distinguishable qubits against spontaneous decay processes. These quantum codes exploit the physical possibilities which are offered by the continuous observation of the distinguishable qubits by photodetectors.

How can we stabilize an $`N`$-qubit system, such as the one depicted in Fig. 1, against spontaneous decay processes if we can observe the qubits continuously with the help of photodetectors? To achieve this two major tasks have to be tackled. Namely, we have to correct the modifications taking place during successive photon emission events. These modifications are described by the effective Hamiltonian of Eq. (2). And, we have to invert each quantum jump $`|1_\alpha |0_\alpha `$ caused by the spontaneous emission of a photon by qubit $`\alpha `$.

In principle, the errors taking place during two successive photon emission events can be corrected by an active error-correcting quantum code with a sufficiently rapid sequence of control measurements and recovery operations Pellizzari . However for practical applications such an approach is not very attractive. One faces all the problems which are related to the implementation of large numbers of control measurements and recovery operations. It is therefore much more desirable to develop passive error-correcting methods for stabilizing the $`N`$-qubit system against modifications arising from the effective Hamiltonian of Eq. (2).

The main idea is to determine an appropriate subspace $`D`$ of the Hilbert space of the $`N`$ distinguishable qubits in which the perturbing part of the Hamiltonian acts as a multiple of the unit operator, i.e.

$$H_{\mathrm{eff}}|_D=H|_D\frac{i}{2}c𝟙|_𝔻$$
(6)
where $`c`$ denotes a c-number. (In our setting this number is real-valued and positive.) If the (unperturbed) dynamics characterized by the Hamiltonian $`H`$ do not take an initial pure quantum state outside this decoherence free subspace $`D`$, this condition implies that, provided no photon is emitted in the time interval $`[0,t]`$, the quantum state at time $`t`$ is identical with the unperturbed state (compare with Eq. (4)). Thus, if one can find such a sufficiently high dimensional decoherence free subspace, the dynamics taking place between successive spontaneous photon emission events are stabilized perfectly without the need for control measurements and recovery operations. In practice it is desirable to choose the dimension of the decoherence free subspace to be as large as possible.

An important special case occurs when all the $`N`$ qubits have identical spontaneous decay rates, i.e. $`\kappa _\alpha \kappa `$. In this situation it follows that

$$\underset{\alpha =1}{\overset{N}{}}L_\alpha ^{}L_\alpha \kappa \underset{\alpha =1}{\overset{N}{}}|1_{\alpha \alpha }1|$$
(7)
and any subspace formed by basis states involving an equal number, say $`wN`$, of excited qubits is a decoherence free subspace. For a given number of $`N`$ qubits the dimension $`K_D`$ of such a decoherence free subspace is given by $`K_D=\left(\genfrac{}{}{0pt}{}{N}{w}\right)`$ which is maximal if $`w=N/2`$. ($`x`$ denotes the largest integer smaller or equal to $`x`$.)

In general, the first spontaneous emission of a photon will affect the quantum state of the $`N`$-qubit system in an irreversible way. According to Eq. (4) the spontaneous emission of a photon by qubit $`\alpha `$, for example, is described by the application of the Lindblad operator $`L_\alpha `$ which induces a quantum jump $`|1_\alpha |0_\alpha `$. This Lindblad operator is not invertible over the decoherence free subspace $`D`$ so that this quantum jump cannot be corrected. In order to correct for this quantum jump actively we have to restrict the dynamics to a still smaller subspace $`CD`$ in which a unitary operator, say $`U_\alpha `$, can be found having the property

$$U_\alpha L_\alpha |_C=𝟙|_{}.$$
(8)
Therefore, if we still want to take advantage of passive error correction between successive photon emission events we have to construct an active error-correcting quantum code within the relevant decoherence free subspace.

We now construct a one detected jump-error correcting embedded quantum code in the special case of identical spontaneous decay rates considered above. According to the criterion given in Knill1 , the orthogonal basis states $`\{|𝒄_i:i=1,\mathrm{},K\}`$ of a subspace $`C`$ constitute an active error-correcting quantum code with respect to the set of error operators $`\{L_\alpha :\alpha =1,\mathrm{},N\}`$ if and only if

$`𝒄_i|L_\alpha ^{}L_\beta |𝒄_j=\delta _{ij}\mathrm{\Lambda }_{\alpha \beta }`$ (9)
for all possible values of $`i,j`$ and $`\alpha ,\beta `$. Eq. (9) states the necessary and sufficient conditions for the existence of unitary recovery operations which fulfill Eq. (8) for the error operators $`\{L_\alpha :\alpha =1,\mathrm{},N\}`$. In the physical setting this criterion states that: (i) different orthogonal quantum states remain orthogonal under the action of error operators; and (ii) all basis states are affected by a given pair of errors $`L_\alpha `$ and $`L_\beta `$ in a similar way. The latter condition necessarily implies that the scalar products between the states $`L_\alpha |𝒄_i`$ and $`L_\beta |𝒄_i`$ are state-independent. It is plausible that a larger set of error operators leads to a more restrictive set of conditions of the type of Eq. (9). Furthermore, we also expect that more restrictive conditions lead to a higher redundancy of an active quantum code.

As an example, consider the situation where continuous observation of the $`N`$-qubit system by photodetectors does not reveal which qubit $`\alpha `$ has emitted the registered photon. This implies that the error operators which could induce a spontaneous decay process are in the set $`\{L_\alpha :\alpha =1,\mathrm{},N\}`$. It has been shown by Plenio et al. Plenio  that when the error positions are unknown, eight physical qubits are needed to encode two orthogonal logical states by an embedded quantum code. This should be compared with the optimal active one-error correcting code using five qubits 5bit .

Thus, the advantage offered by using an embedded quantum code, capable of passively correcting errors between successive photon emission events leads to a significant increase of redundancy in comparison to purely active methods. However, this disadvantage can be overcome if besides knowing the error time, information about the error position $`\alpha `$ is also available. In principle, this information can be obtained from continuous observation of the $`N`$-qubit system by photodetectors as long as the mean distance between adjacent qubits is large in comparison with the wave length of the spontaneously emitted radiation. For this purpose it is important that each photon which is emitted by one of the qubits can be detected.

How can we construct a one detected jump-error correcting embedded quantum code which exploits information about the error position in an optimal way so that its redundancy is minimized? Let us concentrate again on our previously introduced example of identical spontaneous decay rates. In this setting we have a decoherence free subspace which involves $`w`$ excited qubits. This stabilizes the dynamics between successive photon emission events passively. For example, in the simple case of $`N=4`$ and $`w=2`$, the orthogonal basis states $`|𝒄_i`$ of the decoherence free subspace are given by $`\{|1100,|0011,|1010,|0101,|1001,|0110\}`$. Within this six-dimensional decoherence free subspace a possible active quantum code, capable of correcting one quantum jump at a time at a known position $`\alpha `$, is determined by the following three orthonormal code words:

$`|𝒄_0`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}[|1100+e^{i\phi }|0011],`$
$`|𝒄_1`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}[|01010+e^{i\phi }|1010],`$
$`|𝒄_2`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}[|1001+e^{i\phi }|0110].`$ (10)
These code words are formed by all possible different pairs of basis states of the decoherence free subspace which can be constructed by interchanging states $`|0`$ and $`|1`$. These complementary pairings are also characterized by equally weighted probability amplitudes which involve an arbitrary phase $`\phi `$. It can be easily checked that this embedded quantum code fulfills the criterion of Eq. (9) for any of the Lindblad operators $`L_\alpha `$ of Eq. (3) separately, but not for two different Lindblad operators with $`\alpha \beta `$. Thus, provided that the error time and the error position are known this embedded quantum code can correct one quantum jump at a time. In addition, errors between successive quantum jumps are corrected passively. In this way the number of control measurements and recovery operations is reduced and so is the redundancy. This is apparent by comparing our quantum code with the code proposed by Plenio et al. Plenio , which requires eight physical qubits for the encoding of two orthogonal quantum states.

The above construction based on complementary pairings can be generalized to an even number $`N`$ of physical qubits; giving an infinite family of one detected jump-error correcting embedded quantum codes. The construction is the following. We start from the highest dimensional decoherence free subspace which involves all quantum states in which half of the qubits are excited to states $`|1_\alpha `$. Subsequently, we construct an active quantum code within this decoherence free subspace by complementary pairings. The dimension $`K_\mathrm{J}`$ of the resulting embedded quantum code is given by

$`K_\mathrm{J}`$ $`=`$ $`(1/2)\left({\displaystyle \genfrac{}{}{0pt}{}{N}{N/2}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{N1}{N/21}}\right).`$ (11)
It is a remarkable fact that, provided one wants to correct errors passively between successive quantum jumps, it is not possible for an even number of physical qubits $`N`$ to reduce the redundancy of such an embedded jump code (see Sec. V and Ref. DCC ). Therefore, our family of embedded quantum codes has the desirable feature that it minimizes redundancy and that it requires only a small number of recovery operations. For instance, in a given time interval, say $`[0,t]`$, the mean number $`n`$ of required recovery operations is of the order of $`n2t/(N\kappa )`$.

To stabilize any quantum algorithm against spontaneous decay processes using an embedded one detected jump-error correcting quantum code three requirements have to be met. First, one has to be able to register the time $`t`$ and position $`\alpha `$ of each spontaneous decay event which takes place during the performance of the quantum algorithm. As indicated schematically in Fig. 1 this can be achieved by continuous observation of the $`N`$-qubit system with photodetectors. In principle, an identification of the perturbed qubit $`\alpha `$ is possible provided the mean nearest neighbor spacing of the qubits is large in comparison with the wave lengths of the radiation emitted spontaneously. However, in practice the error position $`\alpha `$ might not be determined so easily due to imperfect detection efficiencies of the photodetectors. Therefore, in actual applications shelving techniques Dehmelt ; Toschek  might be useful which amplify each spontaneously emitted photon to such an extent that it can be detected with an efficiency arbitrarily close to unity.

Second, we have to ensure that each spontaneous decay event is corrected immediately by application of an appropriated unitary transformation which inverts the effect of the Lindblad operator $`L_\alpha `$. In practice, this inversion has to be performed on a time scale which is small in comparison with the natural time scale of the quantum algorithm and with the mean spontaneous decay time. Third, one has to ensure that the sequence of quantum gates which constitute the quantum algorithm does not leave the code space at any time. This can be done by encoding the logical information within the code space, and to develop a universal set of quantum gates which leaves this code space invariant.

Ideally these quantum gates are implemented by suitable Hamiltonians. This ensures that the code space is left invariant even during the application of one of these universal quantum gates. Such universal sets of Hamiltonian-induced quantum gates have already been developed for decoherence free subspaces $`D`$ of the kind discussed above Lidar1 ; Zanardi ; Lidar2 . But in general, unitary gates based on swapping Hamiltonians need not be universal on the embedded quantum code, or the swapping Hamiltonians do not leave the embedded quantum code invariant. The solution of this intricate and yet unsolved problem is beyond the scope of the present work. However, some preliminary results have already been obtained recently Fortschr ; QIT ; Khodjasteh .

So far we have shown that any Lindblad operator $`L_\alpha `$ of the form of Eq. (3) can be inverted by our one detected jump-error correcting quantum codes. We provide an example of a unitary transformation which achieves this inversion in the case of the one detected jump-error correcting quantum code involving four physical qubits. A possible sequence of quantum gates capable of inverting a spontaneous decay process affecting qubit $`\alpha =2`$, for example, is depicted in Fig. 2. This example demonstrates the basic fact that it is indeed possible to perform a unitary inversion of the Lindblad operator $`L_2`$ provided Eq. (9) is fulfilled for $`\alpha =\beta =2`$.

## IV General Detected Jump-Error Correcting Quantum Codes

To define general detected jump-error correcting quantum codes, we introduce some notation. For a set of positions $`E=\{\alpha _1,\alpha _2,\mathrm{},\alpha _n\}`$, we denote by $`𝒥_E`$ the operator

$$𝒥_E:=𝒥_{\{\alpha _1,\alpha _2,\mathrm{},\alpha _n\}}:=L_{\alpha _1}L_{\alpha _2}\mathrm{}L_{\alpha _n}.$$

The associated error times $`t_1,\mathrm{},t_n`$ are no longer mentioned explicitly, but it is understood that they are known. Note that the operators $`L_{\alpha _1},\mathrm{},L_{\alpha _n}`$ commute, because the $`\alpha _i`$ are pairwise different. Since by Eq. (5), the errors which involve two equal indices, say $`\alpha _1=\alpha _2`$, cannot occur.

As discussed in Sec. III, for $`H0`$ all the states which are superpositions of states with a constant number of excited qubits are common eigenstates of the non-unitary effective time evolution (6) between quantum jumps. A subspace $`C`$ of such a decoherence free subspace with orthonormal basis $`\{|𝒄_1,\mathrm{},|𝒄_K\}`$ is called a $`d`$-detected jump-error correcting quantum code, and is denoted by $`C=(N,K,d)_w`$ if the following condition holds for sets $`E`$ of jump positions with at most $`d`$ elements and for all basis states $`|𝒄_i`$ and $`|𝒄_j`$:

$$𝒄_i|𝒥_E^{}𝒥_E|𝒄_j=\delta _{ij}\lambda (E).$$
(12)
The notation $`C=(N,K,d)_w`$ is motivated by classical coding theory. Similarly, the notation $`d`$-JC$`(N,w,K)`$ is motivated by notations from design theory (see Alber01 ).

The validity of this statement follows from the general conditions on quantum error-correcting codes (cf. Knill1 ). Since we know on which positions the jump operator $`𝒥_E`$ acts, only products of the form $`𝒥_E^{}𝒥_E`$ have to be considered.

There is a natural connection with combinatorics. For a basis state $`|x_1\mathrm{}x_N`$ of $`N`$ qubits, the positions which are in state $`|1`$ define a subset of $`\{1,\mathrm{},N\}`$. Furthermore, a collection of such subsets corresponds to an equally weighted superposition of basis states.

Let $`^{(1)},\mathrm{},^{(K)}`$ be $`K`$ disjoint sets of subsets $`X\{1,\mathrm{},N\}`$, where each subset $`X`$ contains $`w`$ elements. Identifying the set $`X`$ and the binary word $`𝒙`$ where $`𝒙_i=1`$ if $`iX`$ and $`𝒙_i=0`$ otherwise, we define the states

$$|𝒄_i=\frac{1}{\sqrt{|^{(i)}|}}\underset{X^{(i)}}{}|𝒙.$$
(13)
Where $`|^{(i)}|`$ denotes the number of elements (cardinality) of the set $`^{(i)}`$. These orthonormal states span a $`d`$-detected jump-error correcting quantum code $`C=(N,K,d)_w`$ provided that for all sets of jump positions $`E`$ with no more than $`d`$ elements and all sets $`^{(i)}`$ the following condition holds:

$$\frac{|X^{(i)}:EX\}|}{|^{(i)}|}=\lambda (E).$$
(14)
We note that the disjointness of the sets $`^{(i)}`$ implies condition (12) for $`ij`$. Rewriting the operator $`𝒥_E^{}𝒥_E`$ as $`𝒥_E^{}𝒥_E=_{XE}|𝒙𝒙|`$ shows that for $`i=j`$ and the states (13) the expectation value $`𝒄_i|𝒥_E^{}𝒥_E|𝒄_i`$ equals the expression in (14).

## V Block Designs and Jump Codes

In this section we will show that $`d`$-detected jump-error correcting quantum codes are naturally connected with $`d`$-designs. These are combinatorial structures which have been extensively studied for many decades, cf. Beth1 . To denote this class of combinatorial structures we introduce some notations using the language of finite incidence structures. Let $`V`$ be a set of $`N`$ elements, called points, say $`V=\{1,\mathrm{},N\}`$ where $`N>0`$ is an integer and $`|V|=N`$. For $`0wN`$ the class of $`w`$-subsets of $`V`$ containing $`w`$ elements will be denoted by $`\left(\genfrac{}{}{0pt}{}{V}{w}\right)`$. In a suggestive way its cardinality, i.e.

$$\left|\left(\genfrac{}{}{0pt}{}{V}{w}\right)\right|=\left(\genfrac{}{}{0pt}{}{|V|}{w}\right)=\left(\genfrac{}{}{0pt}{}{N}{w}\right),$$

is just the binomial coefficient $`\left(\genfrac{}{}{0pt}{}{N}{w}\right).`$ An incidence structure in $`V`$ is specified by a distinguished class $``$ of subsets of $`V`$. The elements $`B`$ of $``$ are called blocks (or sometimes lines) of the incidence structure. If $`\left(\genfrac{}{}{0pt}{}{V}{w}\right)`$, we say that $``$ has constant block size $`w`$. As an example, any undirected graph is an incidence structure of block size two, if we choose $`V`$ as the set of points of the graph and the points which are directly connected by an edge as blocks.

For any point $`\alpha V`$ the class $`_\alpha :=\{B:\alpha B\}`$ denotes the class of blocks containing the point $`\alpha `$ (or: “the lines through $`\alpha `$”). If $`|_\alpha |=r`$ is constant for all $`\alpha V`$ the incidence structure $``$ is called regular. For a graph, $`|_\alpha |`$ is the degree $`d(\alpha )`$ of the vertex $`\alpha `$, i.e., the number of edges on which $`\alpha `$ lies. The incidence structure $``$ —as well as the graph itself—is called regular, if $`d(\alpha )`$ is constant. If there exists a constant $`\lambda _d`$ such that for all $`E=\{\alpha _1,\mathrm{},\alpha _d\}\left(\genfrac{}{}{0pt}{}{V}{d}\right)`$ the class $`_E:=(((_{\alpha _1})_{\alpha _2})\mathrm{})_{\alpha _d}`$ of all blocks containing the set $`d`$ of points has size $`|_E|=\lambda _d`$, the incidence structure is called $`d`$-regular. A $`d`$-design is a $`d`$-regular incidence structure with constant block size $`w`$. It is denoted by $`d`$-$`(N,w,\lambda _d)`$-design or as $`S_{\lambda _d}(d,w,N)`$ Beth1 . Graphs which correspond to $`2`$-designs are depicted in Fig. 3.

The preceding discussion leads to the notion of $`d`$-spontaneous emission error designs which we denote by $`d`$-$`\mathrm{SEED}(N,w,K)`$ DCC . The essential property of these combinatorial objects is the local multiplicity $`\lambda (E)`$ of a subset $`E`$ of $`\{1,\mathrm{},N\}`$ containing at most $`d`$ elements which is defined by

$$\frac{|_E^{(i)}|}{|^{(i)}|}=\lambda (E),$$
(15)
where $`^{(1)},\mathrm{},^{(K)}`$ are disjoint subsets of $`\left(\genfrac{}{}{0pt}{}{V}{w}\right)`$. Any $`d`$-$`\mathrm{SEED}(N,w,K)`$ produces a $`d`$-detected jump-error correcting quantum code $`C=(N,K,d)_w`$ using the encoding defined in Eq. (13).

We conclude this section by constructing of a three-jump correcting code. The permutation group $`G=(\mathrm{1\hspace{0.17em}2})(\mathrm{3\hspace{0.17em}4}),(\mathrm{1\hspace{0.17em}4})(\mathrm{2\hspace{0.17em}3}),(\mathrm{5\hspace{0.17em}6})(\mathrm{7\hspace{0.17em}8}),(\mathrm{5\hspace{0.17em}8})(\mathrm{6\hspace{0.17em}7}),(\mathrm{1\hspace{0.17em}2\hspace{0.17em}3})(\mathrm{5\hspace{0.17em}6\hspace{0.17em}7})`$ of order $`48`$ acts on the $`4`$-element subsets of $`\{1,2,\mathrm{},8\}`$. The orbits under $`G`$ of the sets $`\{1,2,5,6\}`$, $`\{1,3,5,6\}`$, and $`\{1,4,5,6\}`$ are mutually disjoint. Direct calculation shows that they fulfill the local multiplicity condition (15). Hence the sets $`^{(1)}:=\{1,2,5,6\}^G`$, $`^{(2)}:=\{1,3,5,6\}^G`$, and $`^{(3)}:=\{1,4,5,6\}^G`$ define an $`(8,3,3)_4`$ jump code. The corresponding (not normalized) basis states are given by:

$`|c_1`$ $`=`$ $`|00110011+|00111100+|01010101+|01011010`$
$`+|01100110+|01101001+|10010110+|10011001`$
$`+|10100101+|10101010+|11000011+|11001100`$
$`|c_2`$ $`=`$ $`|00110110+|00111001+|01010011+|01011100`$
$`+|01100101+|01101010+|10010101+|10011010`$
$`+|10100011+|10101100+|11000110+|11001001`$
$`|c_3`$ $`=`$ $`|00110101+|00111010+|01010110+|01011001`$
$`+|01100011+|01101100+|10010011+|10011100`$
$`+|10100110+|10101001+|11000101+|11001010`$
Further examples of $`d`$-$`\mathrm{SEED}`$s are discussed in DCC . In that paper, there are also general bounds on the parameters of jump codes derived. In particular, the dimension $`K`$ of a $`d`$-detected jump-error correcting $`(N,K,d)_w`$ code is bounded above by

$$K\mathrm{min}\{\left(\genfrac{}{}{0pt}{}{Nd}{wd}\right),\left(\genfrac{}{}{0pt}{}{Nd}{w}\right)\}.$$
(16)
For completeness, we repeat the main ideas of the proof.

The dimension of the space spanned by the basis states of $`N`$ qubits where $`w`$ of them are in the excited state is $`\left(\genfrac{}{}{0pt}{}{N}{w}\right)`$. This implies the bound $`K\left(\genfrac{}{}{0pt}{}{N}{w}\right)`$. A jump on $`j`$ positions reduces the number of excitations to $`wj`$. After the jump, the $`j`$ positions where the jump occurred are zero. There are $`\left(\genfrac{}{}{0pt}{}{Nj}{wj}\right)`$ such basis states. A jump must not reduce the dimension of the code, hence $`K\left(\genfrac{}{}{0pt}{}{Nj}{wj}\right)`$. For $`0jd`$ possible quantum jumps the lowest upper bound is achieved for $`j=d`$, as

$$\left(\genfrac{}{}{0pt}{}{N}{w}\right)=\frac{N(N1)\mathrm{}(Nd+1)}{w(w1)\mathrm{}(wd+1)}\left(\genfrac{}{}{0pt}{}{Nd}{wd}\right).$$

To obtain the second upper bound in Eq. (16), we note that starting with a jump code $`C=(N,K,d)_w`$, applying $`\sigma _x`$ to all qubits yields a jump code $`\overline{C}=(N,K,d)_{Nw}`$. Note that, $`\sigma _x`$ interchanges ground and excited state, hence the code $`\overline{C}`$ lies in the decoherence free subspace with $`Nw`$ excitations. Also, the linear space spanned by the operators $`𝒥_E^{}𝒥_E`$ for all subsets $`E`$ with no more than $`d`$ elements is invariant under conjugation by $`\sigma _x`$ on all qubits. This holds as $`\sigma _xL_\alpha ^{}L_\alpha \sigma _x=|0_{\alpha \alpha }0|=𝟙|\mathrm{𝟙}_{\alpha \alpha }\mathrm{𝟙}|`$. Thus, for the code $`\overline{C}`$ we obtain the bound $`K\left(\genfrac{}{}{0pt}{}{Nd}{Nwd}\right)=\left(\genfrac{}{}{0pt}{}{Nd}{w}\right)`$. If there is no restriction on the number $`w`$ of excited states, choosing $`w=N/2`$ maximizes the upper bound of Eq. (16), i.e.

$$K\left(\genfrac{}{}{0pt}{}{Nd}{N/2d}\right).$$
(17)
As mentioned above, the upper bound of Eq. (17) is achieved for $`d=1`$ and for an even number $`N`$ of qubits. A table of lower bounds (obtained by constructions from $`d`$-$`\mathrm{SEED}`$s) and upper bounds for small values of $`N`$ and $`d`$ is provided in DCC .

## VI Stability properties of the one detected jump-error correcting quantum codes

The detected jump-error correcting quantum codes constructed in the previous sections can stabilize quantum algorithms provided three conditions are fulfilled. First, the decay rates of the qubits are equal. Second, the time and position of each quantum jump is detected with hundred percent efficiency. Third, the appropriate unitary recovery operations are applied perfectly and instantaneously immediately following the detection of an error. Experimentally these requirements can only be approximated. Therefore, the natural question arises how non-ideal conditions influence the robustness of the embedded error-correcting quantum codes.

In this section several types of imperfections affecting the ideal performance of the one detected jump-error correcting embedded quantum codes of Sec. III are studied numerically. For this purpose we investigate the stabilization of a quantum memory and of a simple Hamiltonian dynamics against spontaneous decay processes. The effective two-level Rabi Hamiltonian considered, i.e.

$$H=i\mathrm{}\mathrm{\Omega }(|x_0v||vx_0|),$$
(18)
can be viewed as modeling the quantum dynamics of the ideal Grover search algorithm Grover  in the limit of a large number of qubits. Thereby, the Rabi frequency $`\mathrm{\Omega }`$ can be related to the characteristic time $`\tau _{op}`$ required for performing an oracle-operation and to the number of qubits $`N`$ according to $`\mathrm{\Omega }=2/(\tau _{op}\sqrt{N})`$ Mussinger . Here, $`|v`$ denotes an equally weighted superposition of all (orthonormal) code words which may be used as an initial state in Grover’s quantum search algorithm. The final state we are searching for is denoted $`|x_0`$. For this choice of $`|v`$ the states $`|v`$ and $`|x_0`$ are not orthogonal because $`v|x_0=1/\sqrt{N}`$. However, if the number of qubits $`N`$ becomes large their overlap tends to zero. According to the Hamiltonian Eq. (18) and consistent with Grover’s quantum algorithm Grover , after the interaction time $`\tau =\pi /(2\mathrm{\Omega })`$, the initial state $`|v`$ is transformed to the final state $`|x_0`$.

### VI.1 Imperfect detection of error positions

Let us first of all consider situation where the jump position can be detected with a given non-zero error rate only. Such an imperfection might occur if, for example, a photon emitted by a particular trapped ion is detected by the photodetector associated with a different ion (see with Fig. 4).

The probability to detect the emitted photon at the correct position is denoted by $`p`$. The probability that an emitted photon is detected falsely by the next nearest neighbor is given by $`q\times p`$. Analogously, the probability of detecting the photon by the $`n`$-th nearest neighbor is $`q^n\times p`$ with the normalization condition $`p\times (1+q+q^2+\mathrm{})1`$.

The influence of this type of imperfection on a quantum memory, i.e. $`H0`$ in Eq.(1), is depicted in Fig. 5. A state $`|\psi (0)`$ of the jump code $`C=(4,3,1)_2`$ is propagated according to Eq. (4) with a Monte-Carlo simulation of the quantum-trajectories Zoller . If a jump is detected, the appropriate recovery operation is applied. In the case of a correct detection, the quantum state of the memory is recovered perfectly. In the case of a false detection, the quantum state of the memory leaves the code space. Therefore, in this simulation the full Hilbert space of four physical qubits has to be taken into account. As a measure of fidelity the squared absolute value of the overlap between the state $`|\psi (\tau )`$ after a time $`\tau =\pi /(2\kappa )`$ and the initial state $`|\psi (0)`$ of the memory is plotted as a function of the error parameter $`q`$.

### VI.2 Unequal spontaneous decay rates

For $`\kappa _\alpha \kappa `$ the orthonormal basis states $`|𝒄_i`$ of the code $`CD`$ are degenerate eigenvectors of the operator $`_\alpha L_\alpha ^{}L_\alpha `$ appearing in the effective Hamiltonian of Eq. (2). This property ensures that these states form a passive code for the effective time evolution between successive quantum jumps. The existence of such degenerate eigenstates of the operator $`_\alpha L_\alpha ^{}L_\alpha `$ relies on the assumption that the decay rates of all qubits are equal. Although this physical situation can be realized in a laboratory, it is of interest to investigate what happens if this condition of equal decay rates is violated. In this latter case our code does not correct errors between successive quantum jumps passively.

For this purpose let us consider the Rabi-hamiltonian of Eq. (18) which describes the ideal quantum dynamics. In addition, we assume that the decay rates of the physical qubits are selected randomly according to a Gaussian distribution whose mean value $`\overline{\kappa }`$ is equal to the characteristic Rabi frequency $`\mathrm{\Omega }`$.

To study the resulting time evolution we choose the $`(6,10,1)_3`$ code based on complementary pairings (see Sec. III). In Fig. 6 the fidelity of the quantum state $`\rho (\tau )`$ is depicted as a function of the variance $`\mathrm{\Delta }\kappa `$ of the Gaussian distribution. The fidelity is defined as the overlap $`x_0|\rho (\tau )|x_0`$ between the (mixed) system state $`\rho (\tau )`$ at time $`\tau =\pi /(2\mathrm{\Omega })`$ and the desired state $`|x_0x_0|`$ which would result from the ideal dynamics at this particular interaction time. In this numerical simulation the master equation (1) was integrated up to time $`\tau =\pi /(2\mathrm{\Omega })`$, whereas each jump operator $`L_\alpha `$ was replaced by a sequence consisting of $`L_\alpha `$ and an immediately applied unitary recovery operation $`U_\alpha `$ (see Eq.(8)). In this simulation it was assumed that the recovery operations are performed perfectly.

It is apparent that the code stabilizes the quantum dynamics successfully, despite the fact that the $`(6,10,1)_3`$ code is not a perfect one detected jump-error correcting quantum code for this situation.

### VI.3 Time delay between error detection and recovery operation

Immediately after the detection of a spontaneous emission event the qubits are described by a quantum state belonging to a subspace involving one excited qubit less than the original code space. This subspace also constitutes a passive error-correcting code. Therefore, a time delay between the detection and the application of a recovery operation does not lead to an additional error caused by the effective time evolution $`H_{\mathrm{eff}}`$ provided the ideal quantum dynamics characterized by $`H`$ is not affected. Nevertheless, this time delay must be short in comparison with the mean time between two successive spontaneous emission events. Otherwise, a second spontaneous emission may map the state of the system onto another subspace, from which a recovery is no longer possible. Fig. 7 demonstrates that, as long as the delay between detection and correction is not too large compared with the mean decay time, error correction is still possible.

### VI.4 Dead times of the photodetectors

Another important condition for correct implementation of a detected jump-error correcting quantum code is the ability to observe the environment of each qubit continuously. However, immediately after the detection of a spontaneous emission event, typically the detector is not able to respond to another photon. During the latent response time of a photodetector a second spontaneous emission event can take place which may destroy quantum coherence. In Fig. 8 the dependence of the fidelity $`|\psi (\tau )|x_0|^2`$ on the response time of a photodetector is depicted for various decay rates $`\kappa `$.

It is apparent from Fig. 8 that the detected jump-error correcting quantum code can stabilize an algorithm as long as the response time of the photodetectors is small in comparison with the average time between successive spontaneous emission events.

## VII Conclusions

We have studied quantum error-correcting codes that exploit additional information about the locations of the errors. This information is obtained by continuously monitoring the system. Errors caused by the resulting non-unitary dynamics are corrected passively by embedding the error-correcting code in a decoherence free subspace. To construct such codes we have established connections to design theory. The numerical simulations presented demonstrate that the jump codes discussed can stabilize quantum systems even in cases of imperfect detections and recovery operation.

## VIII Acknowledgments

This work is supported by the DFG (SPP ‘Quanteninformationsverarbeitung’) and by IST-1999-13021 and IST-2001-38869 of the European Commission.


text:
# Production of a pseudo-scalar Higgs boson at hadron colliders at next-to-next-to-leading order

## 1 Introduction

In its minimal version, the Standard Model of particle interactions requires a complex scalar weak-isospin doublet to spontaneously break the electro-weak gauge symmetry. Three of its four original degrees of freedom transform into the longitudinal modes of the $`W^\pm `$ and $`Z`$ gauge bosons while the fourth manifests itself as a neutral scalar field, the Higgs boson. Yukawa-type couplings between the Higgs field and the fermions generate gauge invariant masses for the fermions.

There is no strong theoretical reason to believe that the minimal version of the Standard Model is indeed realized in nature. Extended models give rise to several Higgs particles which can be electrically charged or neutral and which can be even or odd under CP inversion.

The most appealing extension of the Standard Model, at the moment, is the Minimal Supersymmetric Standard Model (MSSM), which basically doubles the field content of the Standard Model and requires two Higgs doublets for the generation of the particle masses. This results in five physical Higgs bosons: three neutral ($`h,H,A`$) and two charged ($`H^\pm `$). $`h`$ and $`H`$ are CP-even while $`A`$ is CP-odd (for a review see Ref. ).

To date, no Higgs boson has been observed in nature, in spite of great efforts that have been made at particle accelerators, especially at the Large Electron Positron collider (LEP) at CERN. The null results lead to lower limits on the mass of possible Higgs bosons. For a scalar Higgs boson in the framework of the minimal Standard Model, this limit is $`M_H114.4`$ GeV at $`95\%`$ CL. Assuming the MSSM, the limit for a neutral scalar Higgs goes down to about 91 GeV. The $`95\%`$ CL limit on a CP-odd Higgs boson, the subject of this paper, is around 92 GeV .

The Large Hadron Collider (LHC) at CERN is scheduled to commence taking data in the year 2007. It will be a proton-proton collider with a center of mass energy of $`\sqrt{s}=14`$ TeV and has been designed specifically to search for Higgs bosons. For all relevant values of the Higgs boson mass and most kinds of Higgs bosons in various models, the gluon-gluon production mode will be of the greatest importance both for discovery and for measuring the Higgs boson mass. For CP-even Higgs bosons, the theoretical prediction is now fairly well under control, having been calculated to next-to-next-to-leading order (NNLO) in the strong coupling constant .

In this paper, we use the techniques of Ref.  to evaluate the production rate for a pseudo-scalar Higgs boson to the same accuracy, i.e., to next-to-next-to-leading order in QCD. We work in the heavy-top limit, using an effective lagrangian for the interaction of the pseudo-scalar Higgs boson with the gluons. As in the case of the scalar Higgs boson, this does not restrict the validity of the results to the Higgs-mass region well below $`2M_t`$ if one factors in the full top mass dependence at leading order .

In the effective lagrangian approach, the massive top-quark loop that mediates the coupling between the Higgs and the constituents of the initial state hadrons, reduces to effective vertices with known coefficient functions. What remains to be computed at NNLO are $`21`$ processes up to two loops, $`22`$ up to one loop, and $`23`$ at tree level, where all internal particles are massless, and all external particles are taken on-shell ($`p_i^2=0`$ for quarks and gluons, $`q^2=M_A^2`$ for the Higgs).

We note that we do not consider contributions from virtual particles in extended theories in this paper. We also do not consider the effect of $`b`$ quark loops at this order, which can be important for example in the MSSM, when the coupling is enhanced by a large value of $`\mathrm{tan}\beta `$. For light quarks, like the $`b`$, the effective lagrangian cannot be formulated and one must perform a true three-loop calculation at this order. That calculation is beyond the current state of the art.

NLO corrections to this process were evaluated in Ref.  and found to be very similar in size to the NLO effects for scalar Higgs production. We find that this is also true at NNLO: The K-factors at NNLO for the scalar and the pseudo-scalar are comparable. This means in turn that the production rate for pseudo-scalar Higgs bosons is fairly well under control. Though still sizable, the NNLO term in the perturbative series is significantly smaller than the NLO term and still higher order effects are presumably negligible. Accordingly, the unphysical dependence of the cross section on the renormalization and factorization scales is reasonably small, allowing for a prediction of the total rate with errors of the order of, or below, the expected experimental precision.

The paper is organized as follows: In the second section, we describe our theoretical framework, including the effective lagrangian, its Wilson coefficients and the renormalization of its operators, and our prescription for handling Levi-Civita tensors and $`\gamma _5`$. In the third section, we briefly discuss the matrix element calculations and the calculation of the partonic cross sections. In section 4, we present our results for the partonic cross sections. The result for pseudo-scalar production is very similar to that for scalar production and the expression for the difference between the two is quite compact. The full expression is presented in the appendix. In section 5 we compute the hadronic cross section by folding in the parton distributions and finally we present our conclusions.

## 2 Theoretical setup

In principle, the whole theoretical background is laid out very clearly in Ref. . Nevertheless, let us review the necessary ingredients for our calculation.

We assume that the pseudo-scalar Higgs boson, $`A`$, couples only to top quarks, $`t`$. The interaction vertex is given by

$$\begin{array}{c}\hfill _{At\overline{t}}=ig_t\frac{A}{v}M_t\overline{t}\gamma _5t,\end{array}$$
(1)
where $`g_t`$ is a coupling constant that depends on the specific theory under consideration. In the MSSM, for example, one has $`g_t=\mathrm{cot}\beta `$, where $`\mathrm{tan}\beta `$ is the standard ratio of the vacuum expectation values of the two Higgs fields.

The coupling of $`A`$ to two gluons is then mediated by a top quark loop. In the heavy-top limit, this interaction can be described by an effective lagrangian :

$$\begin{array}{cc}\hfill _{Agg}& =g_t\frac{A}{v}\left[\stackrel{~}{C}_1^\mathrm{B}\stackrel{~}{𝒪}_1^\mathrm{B}+\stackrel{~}{C}_2^\mathrm{B}\stackrel{~}{𝒪}_2^\mathrm{B}\right],\hfill \\ \hfill \stackrel{~}{𝒪}_1^\mathrm{B}& =G_{\mu \nu }^a\stackrel{~}{G}^{a,\mu \nu },\stackrel{~}{𝒪}_2^\mathrm{B}=_\mu \left(\underset{q}{}\overline{q}\gamma ^\mu \gamma _5q\right).\hfill \end{array}$$
(2)
$`G_{\mu \nu }^a`$ is the gluon field strength tensor, and $`\stackrel{~}{G}_{\mu \nu }^a`$ is its dual:

$$\begin{array}{c}\hfill \stackrel{~}{G}_{\mu \nu }^a=ϵ_{\mu \nu \alpha \beta }G^{a,\alpha \beta }.\end{array}$$
(3)
All quantities in these equations are to be understood as bare quantities in the effective theory of five massless flavors. Thus the sum over quarks ($`_q`$) in Eq. (2) does not include the top quark. $`\stackrel{~}{C}_1^\mathrm{B}`$ and $`\stackrel{~}{C}_2^\mathrm{B}`$ are coefficient functions that can be evaluated perturbatively. One may define renormalized operators and coefficient functions as follows:

$$\begin{array}{c}\hfill \stackrel{~}{𝒪}_i=\underset{j=1}{\overset{2}{}}Z_{ij}\stackrel{~}{𝒪}_j^\mathrm{B},\stackrel{~}{C}_i=\underset{j=1}{\overset{2}{}}(Z^1)_{ji}\stackrel{~}{C}_j^\mathrm{B},i=1,2.\end{array}$$
(4)
The renormalization matrix $`Z`$ is given by

$$\begin{array}{cc}\hfill Z_{11}& =Z_{\alpha _s},Z_{12}=\frac{4}{ϵ}\frac{\alpha _s}{\pi }+𝒪(\alpha _s^2),\hfill \\ \hfill Z_{21}& =0,Z_{22}=Z_{\mathrm{MS}}^sZ_5^s,\hfill \end{array}$$
(5)
with

$$\begin{array}{cc}\hfill Z_{\alpha _s}& =1\frac{\alpha _s}{\pi }\frac{\beta _0}{ϵ}+\left(\frac{\alpha _s}{\pi }\right)^2\left(\frac{\beta _0^2}{ϵ^2}\frac{\beta _1}{2ϵ}\right)+𝒪(\alpha _s^3),\hfill \\ \hfill \beta _0& =\frac{1}{4}\left(11\frac{2}{3}n_f\right),\beta _1=\frac{1}{16}\left(102\frac{38}{3}n_f\right),\hfill \\ \hfill Z_{\mathrm{MS}}^s& =1+\left(\frac{\alpha _s}{\pi }\right)^2\frac{1}{ϵ}\left[\frac{11}{6}+\frac{5}{36}n_f\right]+𝒪(\alpha _s^3),\hfill \\ \hfill Z_5^s& =1\frac{\alpha _s}{\pi }\frac{4}{3}+𝒪(\alpha _s^2),\hfill \end{array}$$
(6)
where $`n_f`$ is the number of light (i.e. massless in our approach) quark flavors. In our numerical analysis below, we will always assume $`n_f=5`$. $`Z_{\alpha _s}`$ and $`Z_{\mathrm{MS}}^s`$ are the renormalization constants of the strong coupling and the singlet axial current, respectively. $`Z_5^s`$ is a finite renormalization constant that will be discussed below.

The coefficient functions are process independent and have been evaluated in the context of pseudo-scalar Higgs decay to NNLO :

$$\begin{array}{cc}\hfill \stackrel{~}{C}_1(\alpha _s)& \frac{\alpha _s}{\pi }\frac{1}{16}+𝒪(\alpha _s^4),\hfill \\ \hfill \stackrel{~}{C}_2(\alpha _s)& =\left(\frac{\alpha _s}{\pi }\right)^2\left(\frac{1}{8}\frac{1}{4}\mathrm{ln}\frac{\mu ^2}{M_t^2}\right)+𝒪(\alpha _s^3),\hfill \end{array}$$
(7)
where $`M_t`$ is the pole mass of the top quark.

The matrix elements to be evaluated have the form

$$\begin{array}{cc}\hfill \stackrel{~}{C}_1\stackrel{~}{𝒪}_1& +\stackrel{~}{C}_2\stackrel{~}{𝒪}_2=\stackrel{~}{C}_1Z_{11}\stackrel{~}{𝒪}_1^\mathrm{B}+(\stackrel{~}{C}_1Z_{12}+\stackrel{~}{C}_2Z_{22})\stackrel{~}{𝒪}_2^\mathrm{B},\hfill \end{array}$$
(8)
with $`\stackrel{~}{𝒪}_n^\mathrm{B}ab|\stackrel{~}{𝒪}_n^\mathrm{B}|XH`$, where $`a`$ and $`b`$ label the two partons in the initial state, and $`X`$ denotes an arbitrary number of partons in the final state. We need to evaluate the square of this matrix element up to $`𝒪(\alpha _s^4)`$.

$`\stackrel{~}{𝒪}_1^\mathrm{B}`$ generates vertices which couple two or three gluons to the pseudo-scalar Higgs (the $`ggggA`$-vertex vanishes due to the Jacobi-identity of the structure functions of SU(3)). Thus, the diagrams related to $`\stackrel{~}{𝒪}_1^\mathrm{B}`$ are the same as in the scalar Higgs case (cf. Ref. ). A sample of typical diagrams is shown in Fig. 1. Higher order contributions are obtained by dressing these diagrams with additional gluons and quarks, either as virtual particles, or as real particles in the final state. The actual results for these diagrams are different from the ones in the scalar case, of course, due to the different Feynman rules corresponding to $`G_{\mu \nu }G^{\mu \nu }`$ and $`G_{\mu \nu }\stackrel{~}{G}^{\mu \nu }`$.

As for $`\stackrel{~}{𝒪}_2^\mathrm{B}`$, since its prefactor in Eq. (8) is of order $`\alpha _s^2`$, it is only required to one order less than $`\stackrel{~}{𝒪}_1^\mathrm{B}`$. Diagrams contributing to $`\stackrel{~}{𝒪}_2^\mathrm{B}`$ for the different subprocesses are shown in Fig. 2. They appear in interference terms with the corresponding diagrams of $`\stackrel{~}{𝒪}_1^\mathrm{B}`$ (see Fig. 1). At NNLO, the square of diagram $`(c)`$ is the only potential contribution to the process arising solely through $`\stackrel{~}{𝒪}_2^\mathrm{B}`$. It vanishes, however, in the limit that the light quarks are massless.

When evaluating the diagrams, the presence of the manifestly 4-dimensional Levi-Civita tensor in $`\stackrel{~}{𝒪}_1^\mathrm{B}`$ as well as the $`\gamma _5`$ in $`\stackrel{~}{𝒪}_2^\mathrm{B}`$ requires special care. We strictly follow the strategy outlined in Ref. . This means that first we replace

$$\begin{array}{cc}\hfill \text{/}q\gamma _5& \frac{i}{3!}q^\alpha ϵ_{\alpha \beta \gamma \delta }\gamma ^\beta \gamma ^\gamma \gamma ^\delta .\hfill \end{array}$$
(9)
Since we need to evaluate squared amplitudes, there will always be exactly two Levi-Civita tensors in our expressions. This allows us to express them in terms of metric tensors, which can be interpreted as $`d`$-dimensional objects:

$$\begin{array}{c}\hfill ϵ_{\alpha \beta \gamma \delta }ϵ^{\overline{\alpha }\overline{\beta }\overline{\gamma }\overline{\delta }}=g_\alpha ^{[\overline{\alpha }}g_\beta ^{\overline{\beta }}g_\gamma ^{\overline{\gamma }}g_\delta ^{\overline{\delta }]}=g_\alpha ^{\overline{\alpha }}g_\beta ^{\overline{\beta }}g_\gamma ^{\overline{\gamma }}g_\delta ^{\overline{\delta }}+g_\alpha ^{\overline{\beta }}g_\beta ^{\overline{\alpha }}g_\gamma ^{\overline{\gamma }}g_\delta ^{\overline{\delta }}\mathrm{},\end{array}$$
(10)
where, as indicated, the brackets around the upper indices indicate that their positions are to be anti-symmetrized while the lower indices remain fixed. All further calculations (i.e., integrations, contraction of indices and spinor algebra) can be performed in $`d`$ dimensions. This procedure is slightly different from that of Ref. , but we have checked that it gives the same result for the pseudo-scalar Higgs decay up to $`𝒪(\alpha _s^4)`$.111We thank M. Steinhauser for providing us with unpublished intermediate results concerning Ref. . Furthermore, we have computed the two-loop virtual corrections for our process in both approaches and obtained identical results.

The finite renormalization constant $`Z_5^s`$, introduced in Eq. (6), is determined by requiring that the one-loop character of the operator relation of the axial anomaly is preserved also at higher orders:

$$\begin{array}{c}\hfill \stackrel{~}{𝒪}_2=\frac{\alpha _s}{\pi }\frac{n_f}{8}\stackrel{~}{𝒪}_1.\end{array}$$
(11)
Note that $`\stackrel{~}{𝒪}_1`$ and $`\stackrel{~}{𝒪}_2`$ denote renormalized operators here, as defined in Eq. (4).

## 3 Methods of evaluation

### 3.1 Virtual corrections

Potential sub-processes are $`ggA`$ and $`q\overline{q}A`$. It turns out, however, that the latter does not contribute at NNLO. The diagrams needed are two-loop vertices with one massive and two massless external legs.

For their evaluation, we use the technique of Ref.  that has already been applied successfully to the evaluation of the virtual correction to scalar Higgs production . This means that the two-loop vertex diagrams are first mapped onto three-loop propagator diagrams by interpreting the massless external lines (gluons or quarks) as part of an additional loop. The resulting integrals, including tensor structures, can be treated by means of the integration-by-parts algorithm of Ref. . In particular, we can use the FORM  program MINCER  as the basis of an algebraic program that reduces all integrals encountered to a set of master integrals . The analytic expressions for the latter have been known for a long time . Note that a generalized version of the method of Ref.  has been constructed in Ref. .

### 3.2 Single real emission

The radiation of one additional parton has to be evaluated up to one-loop level. The contributing processes are $`ggAg`$, $`gqAq`$, $`g\overline{q}A\overline{q}`$, and $`q\overline{q}Ag`$. The one-loop matrix elements have been evaluated to all orders in the dimensional regularization parameter $`ϵ=(4d)/2`$. After interfering with the tree-level amplitude, the squared matrix element can be integrated over single-emission phase space to obtain the contribution to the partonic cross section in closed analytic form. The interference of bare operators $`\stackrel{~}{𝒪}_1^\mathrm{B}`$ and $`\stackrel{~}{𝒪}_2^\mathrm{B}`$ is of order $`ϵ`$. Nonetheless, operator $`\stackrel{~}{𝒪}_2^\mathrm{B}`$ contributes to the single real emission cross section through operator mixing since $`Z_{12}`$ is of order $`1/ϵ`$. The Feynman rules, loop integrals and phase space integration have all been implemented in FORM programs.

### 3.3 Double real emission

We need the tree-level expressions for the processes $`ggAgg`$, $`ggAq\overline{q}`$, $`gqAgq`$, $`q\overline{q}Agg`$, $`q\overline{q}Aq\overline{q}`$, and $`qqAqq`$ (and the corresponding charged conjugated processes). The squared matrix elements can be evaluated straightforwardly in $`d=42ϵ`$ space-time dimensions with the help of FORM. The result is a rather large expression of several thousand terms which must then be integrated over phase space.

The phase space for double real emission is quite complicated and we perform the integration in the method of Ref. . That is, the matrix element and the phase space are expanded in terms of $`(1x)`$, where $`x=M_A^2/\widehat{s}`$. The result is a power series expansion in $`(1x)`$ and $`\mathrm{ln}(1x)`$ (at NNLO, the highest power of the logarithm is $`\mathrm{ln}^3(1x)`$) for the double real emission contribution to the hadronic cross section. If one were to compute all terms in the series, this would be an exact result.

In fact, a truncated series is sufficient for obtaining the NNLO cross section to very high numerical precision. In this case, one obtains a cancellation of infrared singularities by expanding the other contributions to the partonic cross section (single real emission, mass factorization; the dependence of the virtual terms on $`x`$ is simply $`\delta (1x)`$) to the same order in $`(1x)`$ as the double real emission term. In Ref. , we computed scalar Higgs boson production to order $`(1x)^{16}`$ and arrived at a prediction for the hadronic cross section that is phenomenologically equivalent to one based on the closed analytic form of the partonic cross section that has recently become available . This comes as no great surprise. The functions which contribute to the closed analytic form can all be expanded in terms of $`(1x)`$. For example, the dilogarithm can be represented as

$$\begin{array}{c}\hfill \mathrm{Li}_2(x)=\frac{\pi ^2}{6}\underset{n=1}{\overset{\mathrm{}}{}}\frac{(1x)^n}{n}\left(\frac{1}{n}\mathrm{ln}(1x)\right).\end{array}$$
(12)
Furthermore, steeply falling parton distributions ensure that the threshold region dominates and that convergence in $`(1x)`$ is quite rapid.

Although the truncated series leads to a physical result that is by all means equivalent to the exact expression, the approach of expanding in $`(1x)`$ can be taken farther. If the expansion can be evaluated up to sufficiently high order in $`(1x)`$, one can actually invert the series and obtain the partonic cross section in closed analytic form. This is due to the fact that only a limited number of functions appear in this closed form representation: logarithms, dilogarithms, and trilogarithms of various arguments, multiplied by $`1/x`$, $`1/(1+x)`$and $`(1x)^n,n=1,0,1,2,3`$. Taking the functions which appear in the result for the Drell-Yan cross section  and these possible factors, one finds that carrying out the expansion to order $`(1x)^{96}`$ should suffice to allow inversion of the series. In fact, we have carried out the expansion to order $`(1x)^{100}`$ so that we could over-determine the system. As a check, this procedure has also been carried out for the scalar Higgs boson , and complete agreement with the result of Ref.  was found.

It is interesting to observe that the task of inverting the series is much simpler when one examines only the difference between the scalar and pseudo-scalar cases since they are already very similar at the level of squared amplitudes. The difference may be inverted with less than twenty terms.

## 4 Partonic results

As noted above, the partonic cross sections for scalar and pseudo-scalar Higgs production are very similar, so that many terms cancel in the difference between the two. Thus, the partonic cross section for pseudo-scalar Higgs production can be expressed conveniently in terms of the known scalar Higgs boson cross section (with modified normalization) plus a remainder. For this purpose, we write

$$\begin{array}{cc}\hfill \widehat{\sigma }_{ab\mathrm{\Phi }}& =\sigma _\mathrm{\Phi }^0\mathrm{\Delta }_{ab\mathrm{\Phi }},\mathrm{\Phi }\{H,A\},a,b\{g,q,\overline{q}\},\hfill \end{array}$$
(13)
where $`\widehat{\sigma }_{ab\mathrm{\Phi }}`$ is the cross section for the process $`ab\mathrm{\Phi }+X`$. $`a`$ and $`b`$ label the partons in the initial state, $`\mathrm{\Phi }`$ means either a scalar ($`H`$) or pseudo-scalar ($`A`$) Higgs boson, and $`X`$ denotes any number of quarks or gluons in the final state. In the normalization factors, we keep the full top mass dependence:

$$\begin{array}{cc}\hfill \sigma _H^0& =\frac{\pi \sqrt{2}G_\mathrm{F}}{256}\left(\frac{\alpha _s}{\pi }\right)^2\left|\tau _H\left[1+(1\tau _H)f(\tau _H)\right]\right|^2,\hfill \\ \hfill \sigma _A^0& =\frac{\pi \sqrt{2}G_\mathrm{F}}{256}\left(\frac{\alpha _s}{\pi }\right)^2\left|g_t\tau _Af(\tau _A)\right|^2,\tau _\mathrm{\Phi }=\frac{4M_t^2}{M_\mathrm{\Phi }^2},\hfill \end{array}$$
(14)
where $`G_\mathrm{F}1.664\times 10^5\mathrm{GeV}^2`$, and $`g_t`$ has been introduced in Eq. (1). The one-loop function appearing in Eq. (14) is defined by

$$\begin{array}{cc}\hfill f(\tau )& =\{\begin{array}{cc}\mathrm{arcsin}^2\frac{1}{\sqrt{\tau }},\hfill & \tau 1,\hfill \\ \frac{1}{4}\left[\mathrm{ln}\frac{1+\sqrt{1\tau }}{1\sqrt{1\tau }}i\pi \right],\hfill & \tau <1.\hfill \end{array}\hfill \end{array}$$
(15)
For completeness, we note that in the heavy-top limit, the normalization factors approach the values

$$\begin{array}{cc}\hfill \sigma _H^0& \stackrel{M_t\mathrm{}}{}\frac{\pi \sqrt{2}G_\mathrm{F}}{576},\sigma _A^0\stackrel{M_t\mathrm{}}{}\frac{\pi \sqrt{2}G_\mathrm{F}}{256}|g_t|^2.\hfill \end{array}$$
(16)
The kinematic terms are written as a perturbative expansion:

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{ab\mathrm{\Phi }}(x)& =\delta _{ag}\delta _{bg}\delta (1x)+\frac{\alpha _s}{\pi }\mathrm{\Delta }_{ab\mathrm{\Phi }}^{(1)}(x)+\left(\frac{\alpha _s}{\pi }\right)^2\mathrm{\Delta }_{ab\mathrm{\Phi }}^{(2)}(x)+𝒪(\alpha _s^3).\hfill \end{array}$$
(17)
The results for scalar Higgs production up to NNLO ($`\mathrm{\Phi }=H`$) can be found in Refs. .

With the help of these expressions, the corresponding results for the pseudo-scalar Higgs production cross section can be written in a rather compact form. At NLO, the result has been evaluated some time ago :

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{abA}^{(1)}(x)& =\mathrm{\Delta }_{abH}^{(1)}(x)+\frac{1}{2}\delta _{ag}\delta _{bg}\delta (1x),a,b\{g,q,\overline{q}\}.\hfill \end{array}$$
(18)
The main result of this paper are the NNLO terms in $`\mathrm{\Delta }_{abA}`$. We find:

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{ggA}^{(2)}& (x)=\mathrm{\Delta }_{ggH}^{(2)}(x)+\left[\frac{1939}{144}\frac{19}{8}l_{Ht}+3\zeta _2\right]\delta (1x)+6𝒟_1(x)\hfill \\ & (12x6x^2+6x^3)\mathrm{ln}(1x)9x\mathrm{ln}^2(x)\hfill \\ & +\frac{3}{2}\frac{(10x13x^2+4x^32x^4)}{1x}\mathrm{ln}(x)+\frac{(154189x+24x^2+11x^3)}{4}\hfill \\ & +n_f[(\frac{13}{16}\frac{2}{3}l_{Ht}+2\delta _2)\delta (1x)+\frac{2}{3}x\mathrm{ln}^2(x)+x\mathrm{ln}(x)\hfill \\ & \frac{(111x+10x^2)}{6}],\text{with}\delta _2=\frac{1}{4}+\frac{1}{2}l_{Ht},\hfill \\ \hfill \mathrm{\Delta }_{gqA}^{(2)}& (x)=\mathrm{\Delta }_{gqH}^{(2)}(x)+\frac{(44x+2x^2)}{3}\mathrm{ln}(1x)\frac{28}{9}x\mathrm{ln}^2(x)\hfill \\ & +\frac{(22+30xx^2)}{3}\mathrm{ln}(x)+\frac{(337382x+51x^2)}{18},\hfill \\ \hfill \mathrm{\Delta }_{qqA}^{(2)}& (x)=\mathrm{\Delta }_{qqH}^{(2)}(x)\frac{64}{27}x\mathrm{ln}^2(x)+\frac{16}{27}(6+11x)\mathrm{ln}(x)+\frac{8}{27}(3740x+3x^2),\hfill \\ \hfill \mathrm{\Delta }_{q\overline{q}A}^{(2)}& (x)=\mathrm{\Delta }_{q\overline{q}H}^{(2)}(x)+\frac{32}{27}x\mathrm{ln}^2(x)+\frac{32}{27}(3+8x)\mathrm{ln}(x)\hfill \\ & +\frac{16}{27}(11x9x^2x^3)+n_f\left[\frac{32}{27}x\mathrm{ln}(x)\frac{16}{27}(1x^2)\right],\hfill \\ \hfill \mathrm{\Delta }_{qq^{}A}^{(2)}& (x)=\mathrm{\Delta }_{qq^{}H}^{(2)}(x)\frac{16}{9}x\mathrm{ln}^2(x)+\frac{16}{9}(2+3x)\mathrm{ln}(x)+\frac{8}{9}(1112x+x^2),\hfill \end{array}$$
(19)
where $`l_{Ht}\mathrm{ln}(M_A^2/M_t^2)`$, $`\zeta _n\zeta (n)`$ is Riemann’s $`\zeta `$ function ($`\zeta _2=\pi ^2/6=1.64483\mathrm{}`$, $`\zeta _3=1.20206\mathrm{}`$) and $`𝒟_n(x)[\mathrm{ln}^n(1x)/(1x)]_+`$. Of course one must change $`M_HM_A`$ wherever it appears in $`\mathrm{\Delta }_{abH}`$. For the sake of completeness, we list the full result for $`\mathrm{\Delta }_{abA}`$ in App. A. The difference $`\mathrm{\Delta }_{abA}\mathrm{\Delta }_{abH}`$ can be expanded readily in terms of $`(1x)`$ in order to bring it to a form consistent with the results of Ref. .

The only contribution that originates from the presence of the (renormalized!) operator $`\stackrel{~}{𝒪}_2`$ is the term $`\delta _2`$ in the equations above. We obtain it by computing the interference of diagram $`(a)`$ in Fig. 1 with diagram $`(a)`$ of Fig. 2. Alternatively, it can be derived from Eq. (11) by using the LO result for $`\mathrm{\Delta }_{ggA}`$:

$$\begin{array}{cc}\hfill (\stackrel{~}{C}_1\stackrel{~}{C}_2)^1n_f\delta _2\left(\frac{\alpha _s}{\pi }\right)^4\delta (1x)& =\frac{\alpha _s}{\pi }\frac{n_f}{8}\left[\stackrel{~}{C}_1^2\left(\frac{\alpha _s}{\pi }\right)^2\delta (1x)\right],\hfill \end{array}$$
(20)
which leads to the same result for $`\delta _2`$ as above. This provides a welcome check on the normalization of the contribution from $`\stackrel{~}{𝒪}_2`$.

## 5 Hadronic results

In complete analogy to scalar Higgs production, one has to convolute the partonic rate with the proper parton distribution functions, in order to arrive at a physical prediction for the hadronic production rate:

$$\begin{array}{cc}\hfill \sigma _{h_1h_2}(s)& =\underset{a,b}{}_0^1𝑑\tau _\tau ^1\frac{dx_a}{x_a}\left[f_{a/h_1}(x_a)f_{b/h_2}(\tau /x_a)\right]\widehat{\sigma }_{ab}(\widehat{s}=s\tau ).\hfill \end{array}$$
(21)
A consistent NNLO hadronic result requires not only a NNLO partonic cross section, but also PDFs that have been evaluated at NNLO. Strictly speaking, such a set of PDFs is not yet available, because the NNLO evolution equation is not fully known. Nevertheless, we use an approximate set of PDF, based on approximations to the evolution equation  derived from the available moments of the structure functions . At lower order, we use the corresponding lower order PDF.

Fig. 3 shows the total cross section as a function of the pseudo-scalar Higgs boson mass, $`M_A`$, at $`(a)`$ the LHC, and $`(b)`$ Tevatron Run II. Fig. 4 shows the curves for the LHC over a larger mass range. The renormalization and factorization scales have been set to $`M_A`$. The cusp in Fig. 4 is a leading order effect caused by the $`t\overline{t}`$ threshold, cf. Eq. (14). Since we do not want to specify a particular extension of the Standard Model, we choose the coupling of the pseudo-scalar Higgs to the top quarks such that $`g_t=1`$ (cf. Eq. (1)). The total cross section scales with $`|g_t|^2`$, so that the actual numbers can be easily obtained from the ones shown in the figure. As already noted, the behavior of the corrections is very similar to the ones for scalar Higgs production: The NNLO corrections are significantly smaller than the NLO corrections, indicating a nicely converging result with uncertainties well under control.

This is further affirmed in Figs. 5 and 6, which show the variation of the cross section with the renormalization and factorization scales for a fixed Higgs mass $`M_A=115`$ GeV. The scales are varied between $`M_A/4`$ and $`4M_A`$: in sub-panels $`(a)`$, the renormalization scale and the factorization scale are identified and varied simultaneously; in sub-panels $`(b)`$, the renormalization scale $`\mu _R`$ is identified with $`M_A`$, and $`\mu _F`$ is varied, while in sub-panels $`(c)`$, $`\mu _F=M_A`$ is fixed and $`\mu _R`$ is varied. One observes a clear reduction of the scale dependence at NNLO with respect to NLO. Compared to the LO curve, there is a significant relative improvement in the scale variation, while the absolute ranges of variation at LO and NNLO are comparable. It is also clear that the scale dependence is dominated by renormalization scale dependence; factorization scale dependence is quite small. Using the rather conservative range of $`M_A/4<\mu <4M_A`$, one arrives at an uncertainty in $`\sigma `$ from scale variation of about $`\pm 30\%`$ ($`\pm 65\%`$) at LO, $`\pm 30\%`$ ($`\pm 50\%`$) NLO, and less than $`\pm 20\%`$ ($`\pm 25\%`$) at NNLO for the LHC (numbers in brackets for Tevatron Run II). Varying $`\mu `$ between $`M_A/2`$ and $`2M_A`$ results in a variation of $`\sigma `$ of $`\pm 20\%`$ ($`\pm 40\%`$) at LO, $`\pm 15\%`$ ($`\pm 25\%`$) at NLO, and less than $`\pm 10\%`$ ($`\pm 15\%`$) at NNLO.

## 6 Conclusions

The hadronic cross section for the production of a pseudo-scalar Higgs boson has been calculated at NNLO in QCD. We find that the corrections are similar to the scalar case, both in their magnitude and in their uncertainty due to scale dependence. While this uncertainty is still rather large, it seems that the NNLO calculation yields a reliable prediction for the total rate.

#### Note Added:

As we completed this manuscript, we became aware of a similar paper by Anastasiou and Melnikov. We have compared results for the partonic cross sections and find complete agreement.

#### Acknowledgments.

R.V.H. thanks the High Energy Theory group at Brookhaven National Laboratory, where part of this work has been done, for hospitality. The work of W.B.K. was supported by the U.S. Department of Energy under Contract No. DE-AC02-98CH10886.

## Appendix A Analytic results

In this appendix, the analytic expressions for the partonic cross section of pseudo-scalar Higgs production at NNLO are listed. Note that using the formulas of Eqs. (18) and (19), one can transform these expressions into the corresponding ones for scalar Higgs production.

At NLO , the result for the gluon-gluon sub-process is

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{ggA}^{(1)}(x)& =(6+6\zeta _2)\delta (1x)+12𝒟_1(x)(24x12x^2+12x^3)\mathrm{ln}(1x)\hfill \\ & 6\frac{(1x+x^2)^2}{1x}\mathrm{ln}(x)\frac{11}{2}(1x)^3.\hfill \end{array}$$
(22)
For the quark-gluon channel one finds

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{gqA}^{(1)}(x)& =\frac{(88x+4x^2)}{3}\mathrm{ln}(1x)\frac{(44x+2x^2)}{3}\mathrm{ln}(x)\frac{(36x+x^2)}{3},\hfill \end{array}$$
(23)
and for the quark–anti-quark channel

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{q\overline{q}A}^{(1)}(x)=& \frac{32}{27}(1x)^3.\hfill \end{array}$$
(24)
The contributions from different sub-processes at NNLO are written as

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{abA}^{(2)}& =\mathrm{\Delta }_{abA}^A+n_f\mathrm{\Delta }_{abA}^F,\hfill \end{array}$$
(25)
where $`n_f`$ is the number of light (i.e., massless in our approach) quark flavors.

For the sub-process with two gluons in the initial state, we find

$$\begin{array}{cc}& \mathrm{\Delta }_{ggA}^A(x)=\left[\frac{741}{8}+\frac{139}{2}\zeta _2\frac{165}{4}\zeta _3\frac{9}{20}\zeta _2^2\right]\delta (1x)\hfill \\ & \left[\frac{101}{3}33\zeta _2\frac{351}{2}\zeta _3\right]𝒟_0(x)+\left[13990\zeta _2\right]𝒟_1(x)33𝒟_2(x)+72𝒟_3(x)\hfill \\ & (144x72x^2+72x^3)\mathrm{ln}^3(1x)(297381x+348x^2330x^3)\mathrm{ln}^2(1x)\hfill \\ & \frac{9}{2}\frac{(3130x+93x^294x^3+31x^4)}{1x}\mathrm{ln}^2(1x)\mathrm{ln}(x)\hfill \\ & +\left[\frac{(20272735x+2182x^22583x^3)}{4}+(180x90x^2+90x^3)\zeta _2\right]\mathrm{ln}(1x)\hfill \\ & +3\frac{(88211x+312x^2365x^3+187x^4)}{1x}\mathrm{ln}(1x)\mathrm{ln}(x)\hfill \\ & +9\frac{(7+3x+19x^23x^319x^4+9x^5)}{1x^2}\mathrm{ln}(1x)\mathrm{ln}^2(x)\hfill \\ & +36\frac{(16x13x^26x^3+x^4)}{1+x}\mathrm{ln}(1x)\mathrm{Li}_2(1x)\hfill \\ & 18\frac{(1+2x+3x^2+2x^3+x^4)}{1+x}\mathrm{ln}(1x)\mathrm{Li}_2(1x^2)\hfill \\ & \frac{9}{2}\frac{(2438x+8x^2+54x^319x^4+9x^5)}{1x^2}\mathrm{Li}_3(1x)\hfill \\ & \frac{9}{2}\frac{(27+35x+75x^229x^378x^4+6x^5)}{1x^2}\mathrm{Li}_3\left(\frac{(1x)}{x}\right)\hfill \\ & +\frac{9}{8}\frac{(1+2x+3x^2+24x^3+16x^4)}{1+x}\mathrm{Li}_3(1x^2)\hfill \\ & \frac{9}{8}\frac{(1+2x+3x^28x^38x^4)}{1+x}\mathrm{Li}_3\left(\frac{(1x^2)}{x^2}\right)\hfill \\ & \frac{9}{2}\frac{(7+14x+21x^2+8x^3+4x^4)}{1+x}\left[\mathrm{Li}_3\left(\frac{1x}{1+x}\right)\mathrm{Li}_3\left(\frac{1x}{1+x}\right)\right]\hfill \\ & \frac{3}{4}\frac{(317398x87x^2+300x^3121x^4)}{1x}\mathrm{Li}_2(1x)\hfill \\ & +\frac{9}{2}\frac{(11+31x+59x^225x^365x^4+11x^5)}{1x^2}\mathrm{Li}_2(1x)\mathrm{ln}(x)\hfill \\ & +\frac{(42+36x63x^233x^3)}{4}\mathrm{Li}_2(1x^2)+\frac{9}{4}\frac{(5+10x+15x^22x^4)}{1+x}\mathrm{Li}_2(1x^2)\mathrm{ln}(x)\hfill \\ & +\frac{3}{4}\frac{(21+23x+41x^237x^344x^4+4x^5)}{1x^2}\mathrm{ln}^3(x)\hfill \\ & \frac{3}{8}\frac{(154365x+675x^2827x^3+374x^4)}{1x}\mathrm{ln}^2(x)\hfill \\ & \frac{1}{8}\frac{(22135599x+6603x^27003x^3+4342x^4)}{1x}\mathrm{ln}(x)\hfill \\ & +9\frac{(92x+27x^234x^3+9x^4)}{1x}\zeta _2\mathrm{ln}(x)\frac{(1630920611x+23819x^222749x^3)}{48}\hfill \\ & +\frac{3}{4}(319277x+233x^2363x^3)\zeta _2\frac{351}{2}(2xx^2+x^3)\zeta _3\hfill \end{array}$$
(26)
and

$$\begin{array}{cc}& \mathrm{\Delta }_{ggA}^F(x)=\left[\frac{689}{72}+l_{Ht}\frac{5}{3}\zeta _2+\frac{5}{6}\zeta _3\right]\delta (1x)\hfill \\ & +\left[\frac{14}{9}2\zeta _2\right]𝒟_0(x)\frac{10}{3}𝒟_1(x)+2𝒟_2(x)\hfill \\ & +\frac{2}{9}(812x+3x^217x^3)\mathrm{ln}^2(1x)+\frac{8}{3}(x+x^2)\mathrm{ln}^2(1x)\mathrm{ln}(x)\hfill \\ & \frac{(922294x+249x^21570x^3)}{108}\mathrm{ln}(1x)\hfill \\ & \frac{2}{9}\frac{(17+7x+21x^261x^3+25x^4)}{1x}\mathrm{ln}(1x)\mathrm{ln}(x)\hfill \\ & \frac{8}{3}(x+x^2)\mathrm{ln}^2(x)\mathrm{ln}(1x)+\frac{16}{3}(x+x^2)\mathrm{ln}(1x)\mathrm{Li}_2(1x)\hfill \\ & \frac{(2+14x+17x^2)}{6}\mathrm{Li}_3(1x)\frac{(234x31x^2)}{12}\mathrm{Li}_3\left(\frac{(1x)}{x}\right).\hfill \\ & +\frac{1}{36}\frac{(68302x+21x^2+227x^3+4x^4)}{1x}\mathrm{Li}_2(1x)\hfill \\ & +\frac{(250x47x^2)}{12}\mathrm{Li}_2(1x)\mathrm{ln}(x)\hfill \\ & +\frac{(2+6x+9x^2)}{72}\mathrm{ln}^3(x)+\frac{1}{72}\frac{(68+100x+69x^2351x^3+132x^4)}{1x}\mathrm{ln}^2(x)\hfill \\ & +\frac{1}{216}\frac{(1282382x+117x^23041x^3+2384x^4)}{1x}\mathrm{ln}(x)\frac{8}{3}(x+x^2)\zeta _2\mathrm{ln}(x)\hfill \\ & +\frac{(12707606x+1641x^217774x^3)}{1296}\frac{2}{9}(812x+3x^217x^3)\zeta _2\hfill \end{array}$$
(27)
For the quark-gluon channel, we have:

$$\begin{array}{cc}& \mathrm{\Delta }_{gqA}^A(x)=\frac{367}{54}(22x+x^2)\mathrm{ln}^3(1x)\hfill \\ & \frac{(25922278x111x^2288x^3)}{36}\mathrm{ln}^2(1x)\frac{(642+190x+553x^2)}{18}\mathrm{ln}^2(1x)\mathrm{ln}(x)\hfill \\ & +\left[\frac{(2388717388x2538x^2784x^3)}{162}\frac{50}{9}(22x+x^2)\zeta _2\right]\mathrm{ln}(1x)\hfill \\ & +\frac{(16652040x+174x^2400x^3)}{27}\mathrm{ln}(1x)\mathrm{ln}(x)+\frac{4(38+21x+39x^2)}{9}\mathrm{ln}(1x)\mathrm{ln}^2(x)\hfill \\ & \frac{2}{9}(46+298x+139x^2)\mathrm{ln}(1x)\mathrm{Li}_2(1x)2(2+2x+x^2)\mathrm{ln}(1x)\mathrm{Li}_2(1x^2)\hfill \\ & \frac{2}{9}(42142xx^2)\mathrm{Li}_3(1x)\frac{(302+474x+339x^2)}{9}\mathrm{Li}_3\left(\frac{(1x)}{x}\right)\hfill \\ & \frac{(2+2x+x^2)}{2}\mathrm{Li}_3(1x^2)\frac{(2+2x+x^2)}{2}\mathrm{Li}_3\left(\frac{(1x^2)}{x^2}\right)\hfill \\ & 4(2+2x+x^2)\left[\mathrm{Li}_3\left(\frac{1x}{1+x}\right)\mathrm{Li}_3\left(\frac{1x}{1+x}\right)\right]\hfill \\ & \frac{(979144x215x^2+52x^3)}{18}\mathrm{Li}_2(1x)+\frac{(142+374x+245x^2)}{9}\mathrm{Li}_2(1x)\mathrm{ln}(x)\hfill \\ & +\frac{(166+222x+33x^2+4x^3)}{54}\mathrm{Li}_2(1x^2)+2(2+2x+x^2)\mathrm{ln}(x)\mathrm{Li}_2(1x^2)\hfill \\ & +\frac{(133+202x+115x^2)}{27}\mathrm{ln}^3(x)\frac{(8371296x+234x^2226x^3)}{54}\mathrm{ln}^2(x)\hfill \\ & \left[\frac{(2204232040x5847x^22464x^3)}{324}\frac{(194+222x+213x^2)}{9}\zeta _2\right]\mathrm{ln}(x)\hfill \\ & \frac{(173719156324x12687x^26148x^3)}{1944}+\frac{(1071710x130x^2144x^3)}{18}\zeta _2\hfill \\ & +\frac{311}{18}(22x+x^2)\zeta _3\hfill \end{array}$$
(28)
and

$$\begin{array}{cc}& \mathrm{\Delta }_{gqA}^F(x)=\frac{(22x+x^2)}{18}\mathrm{ln}^2(1x)\hfill \\ & \frac{(1316x+9x^2)}{9}\mathrm{ln}(1x)\frac{(44x+2x^2)}{9}\mathrm{ln}(1x)\mathrm{ln}(x)\hfill \\ & +\frac{(22x+x^2)}{9}\mathrm{ln}^2(x)+\frac{(2938x+19x^2)}{27}\mathrm{ln}(x)+\frac{(265418x+179x^2)}{162}.\hfill \end{array}$$
(29)
For the scattering of two identical quarks, we find

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{qqA}^A& (x)=\frac{32}{9}(32xx^2)\mathrm{ln}^2(1x)\frac{16}{9}(4+4x+x^2)\mathrm{ln}^2(1x)\mathrm{ln}(x)\hfill \\ & +\frac{4}{3}(1712x5x^2)\mathrm{ln}(1x)+\frac{8}{9}(128x5x^2)\mathrm{ln}(1x)\mathrm{ln}(x)\hfill \\ & +\frac{8}{9}(4+4x+x^2)\mathrm{ln}(1x)\mathrm{ln}^2(x)\frac{32}{9}(4+4x+x^2)\mathrm{ln}(1x)\mathrm{Li}_2(1x)\hfill \\ & \frac{8}{27}(22x+x^2)\mathrm{Li}_3(1x)\frac{8}{27}(50+46x+13x^2)\mathrm{Li}_3\left(\frac{(1x)}{x}\right)\hfill \\ & \frac{8}{9}(64xx^2)\mathrm{Li}_2(1x)+\frac{16}{27}(19+17x+5x^2)\mathrm{Li}_2(1x)\mathrm{ln}(x)\hfill \\ & +\frac{8}{81}(19+17x+5x^2)\mathrm{ln}^3(x)\frac{4}{27}(1810x9x^2)\mathrm{ln}^2(x)\hfill \\ & \frac{2}{27}(129212x69x^2)\mathrm{ln}(x)+\frac{16}{9}(4+4x+x^2)\zeta _2\mathrm{ln}(x)\hfill \\ & +\frac{4}{27}(86+53x+33x^2)+\frac{32}{9}(32xx^2)\zeta _2\hfill \end{array}$$
(30)
and

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{qqA}^F& (x)=0.\hfill \end{array}$$
(31)
For the scattering of a quark–(anti-)quark pair of distinct flavor, we find

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{qq^{}A}^A& (x)=\frac{32}{9}(32xx^2)\mathrm{ln}^2(1x)\frac{16}{9}(4+4x+x^2)\mathrm{ln}^2(1x)\mathrm{ln}(x)\hfill \\ & +\frac{4}{3}(1712x5x^2)\mathrm{ln}(1x)+\frac{8}{9}(128x5x^2)\mathrm{ln}(1x)\mathrm{ln}(x)\hfill \\ & +\frac{8}{9}(4+4x+x^2)\mathrm{ln}(1x)\mathrm{ln}^2(x)\frac{32}{9}(4+4x+x^2)\mathrm{ln}(1x)\mathrm{Li}_2(1x)\hfill \\ & \frac{32}{9}(4+4x+x^2)\mathrm{Li}_3\left(\frac{(1x)}{x}\right)\hfill \\ & \frac{8}{9}(64xx^2)\mathrm{Li}_2(1x)+\frac{8}{3}(4+4x+x^2)\mathrm{Li}_2(1x)\mathrm{ln}(x)\hfill \\ & +\frac{4}{9}(4+4x+x^2)\mathrm{ln}^3(x)\frac{8}{9}(32xx^2)\mathrm{ln}^2(x)\hfill \\ & \frac{2}{9}(4368x29x^2)\mathrm{ln}(x)+\frac{16}{9}(4+4x+x^2)\zeta _2\mathrm{ln}(x)\hfill \\ & \frac{2}{9}(6146x15x^2)+\frac{32}{9}(32xx^2)\zeta _2\hfill \end{array}$$
(32)
and

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{qq^{}A}^F& (x)=0.\hfill \end{array}$$
(33)
For quark–anti-quark scattering (same quark flavor), we have

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{q\overline{q}A}^A& (x)=\frac{32}{81}(14+21x48x^2+13x^3)\mathrm{ln}^2(1x)\frac{16}{9}(4+4x+x^2)\mathrm{ln}^2(1x)\mathrm{ln}(x)\hfill \\ & +\frac{4}{81}(75+892x1351x^2+384x^3)\mathrm{ln}(1x)\hfill \\ & +\frac{8}{81}(76+72x189x^2+64x^3)\mathrm{ln}(1x)\mathrm{ln}(x)\hfill \\ & +\frac{8}{9}(4+4x+x^2)\mathrm{ln}(1x)\mathrm{ln}^2(x)\frac{32}{9}(4+4x+x^2)\mathrm{ln}(1x)\mathrm{Li}_2(1x)\hfill \\ & \frac{40}{27}(2+2x+x^2)\mathrm{Li}_3(1x)\frac{8}{9}(18+18x+5x^2)\mathrm{Li}_3\left(\frac{(1x)}{x}\right)\hfill \\ & +\frac{10}{27}(2+2x+x^2)\mathrm{Li}_3(1x^2)+\frac{2}{9}(2+2x+x^2)\mathrm{Li}_3\left(\frac{(1x^2)}{x^2}\right)\hfill \\ & +\frac{8}{27}(2+2x+x^2)\left[\mathrm{Li}_3\left(\frac{1x}{1+x}\right)\mathrm{Li}_3\left(\frac{1x}{1+x}\right)\right]\hfill \\ & \frac{8}{81}(12+30x93x^226x^3)\mathrm{Li}_2(1x)+\frac{16}{9}(7+7x+2x^2)\mathrm{Li}_2(1x)\mathrm{ln}(x)\hfill \\ & \frac{8}{27}(2x+x^2+6x^3)\mathrm{Li}_2(1x^2)\frac{4}{9}(2+2x+x^2)\mathrm{Li}_2(1x^2)\mathrm{ln}(x)\hfill \\ & +\frac{8}{27}(5+5x+x^2)\mathrm{ln}^3(x)\frac{8}{81}(27+27x81x^2+44x^3)\mathrm{ln}^2(x)\hfill \\ & +\frac{2}{81}(91064x+2111x^2768x^3)\mathrm{ln}(x)+\frac{16}{9}(4+4x+x^2)\zeta _2\mathrm{ln}(x)\hfill \\ & +\frac{20}{81}(101462x+520x^2159x^3)+\frac{16}{81}(11+93x147x^2+43x^3)\zeta _2\hfill \end{array}$$
(34)
and

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{q\overline{q}A}^F& (x)=\frac{32(1x)^3}{81}\mathrm{ln}(1x)\frac{16(39x+12x^24x^3)}{81}\mathrm{ln}(x)\hfill \\ & \frac{8(41111x+111x^241x^3)}{243}.\hfill \end{array}$$
(35)
In the expressions above, we have identified the renormalization and the factorization scale with the Higgs boson mass, $`\mu _R=\mu _F=M_A`$. The dependence on these scales is logarithmic and can readily be reconstructed by employing scale invariance of the total partonic cross section.