text:
# DESY 01-015 ISSN 0418-9833 hep-ph/0103018 February 2001 Pair production of neutral Higgs bosons at the CERN Large Hadron Collider

## 1 Introduction

The search for Higgs bosons will be among the prime tasks of the CERN Large Hadron Collider (LHC) . While the standard model (SM) of elementary-particle physics contains one complex Higgs doublet, from which one neutral $`CP`$-even Higgs boson $`H`$ emerges in the physical particle spectrum after the spontaneous breakdown of the electroweak symmetry, the Higgs sector of the minimal supersymmetric extension of the SM (MSSM) consists of a two-Higgs-doublet model (2HDM) and accommodates five physical Higgs bosons: the neutral $`CP`$-even $`h^0`$ and $`H^0`$ bosons, the neutral $`CP`$-odd $`A^0`$ boson, and the charged $`H^\pm `$-boson pair. At the tree level, the MSSM Higgs sector has two free parameters, which are usually taken to be the mass $`m_{A^0}`$ of the $`A^0`$ boson and the ratio $`\mathrm{tan}\beta =v_2/v_1`$ of the vacuum expectation values of the two Higgs doublets.

In the SM, the trilinear and quartic Higgs-boson couplings are proportional to the square of the Higgs-boson mass, and they are uniquely fixed once the latter is known. In the MSSM, the various Higgs-boson self-couplings are determined by the gauge couplings multiplied with trigonometric factors that depend on $`\alpha `$, the mixing angle that rotates the weak $`CP`$-even Higgs eigenstates into the mass eigenstates $`h^0`$ and $`H^0`$, and $`\beta `$  \[see Eq. (A.1)\].

At the LHC, the trilinear Higgs-boson couplings may be probed by studying the inclusive hadroproduction of Higgs-boson pairs. These couplings enter the stage via the Feynman diagrams where a virtual neutral Higgs boson is produced in the $`s`$ channel through $`q\overline{q}`$ annihilation or $`gg`$ fusion and in turn decays into a pair of neutral or charged Higgs bosons. Specifically, the partonic subprocesses include $`q\overline{q}HH`$ and $`ggHH`$ in the SM, and $`q\overline{q}\varphi _1\varphi _2,H^+H^{}`$ and $`gg\varphi _1\varphi _2,H^+H^{}`$ in the MSSM, where $`\varphi _i=h^0,H^0,A^0`$. At the tree level, there are two mechanisms of $`q\overline{q}`$ annihilation. On the one hand, it can proceed via a $`Z`$ boson (Drell-Yan process) if an appropriate Higgs-Higgs-$`Z`$ coupling exists. On the other hand, the Higgs bosons can be radiated off the $`q`$-quark line if the relevant Yukawa couplings are sufficiently strong. In the SM, $`q\overline{q}`$ annihilation is greatly suppressed due to the absence of a $`HHZ`$ coupling and the smallness of the $`Hq\overline{q}`$ couplings for the active quarks contained inside the proton, $`q=u,d,s,c,b`$. In the MSSM, however, we have $`h^0A^0Z`$ and $`H^0A^0Z`$ couplings at the tree level \[see Eq. (A.2)\], and the $`\varphi _ib\overline{b}`$ couplings are generally strong if $`\mathrm{tan}\beta `$ is large \[see Eq. (A.3)\]. In the SM, $`gg`$ fusion is mediated via heavy-quark loops. In the MSSM, there are additional contributions from squark loops. The SM case was studied in Ref. . As for $`H^+H^{}`$ pair production, $`q\overline{q}`$ annihilation was investigated in Refs. , the quark loop contribution to $`gg`$ fusion in Refs. , and the squark loop one in Refs. . As for the pair production of neutral Higgs bosons in the MSSM, $`q\overline{q}`$ annihilation via a $`Z`$ boson was analyzed in Ref.  and the quark and squark loop contributions to $`gg`$ fusion in Refs.  and , respectively. In Ref. , also the QCD corrections to $`q\overline{q}`$ annihilation via a $`Z`$ boson and to $`gg`$ fusion via an infinitely heavy top quark were considered. The processes of associated production of a neutral Higgs-boson pair with a dijet (in addition to the remnant jets), an intermediate boson, or another neutral Higgs boson were found to have cross sections that are greatly suppressed, by more than an order of magnitude, compared to those of the respective processes without those additional final-state particles .

The purpose of this paper is to reanalyze the pair production of neutral Higgs bosons in the MSSM, both via $`q\overline{q}`$ annihilation and $`gg`$ fusion. In the case of $`q\overline{q}`$ annihilation, we also allow for $`q=b`$. This makes it necessary to include a new class of diagrams involving $`b`$-quark Yukawa couplings, which are depicted in Fig. 1(a) and the second and third lines of Fig. 1(b). These come in addition to the Drell-Yan diagram, shown in the first line of Fig. 1(b), which is already present for $`q=u,d,s,c`$. For final states with an even number of $`CP`$-even Higgs bosons, $`\varphi _1\varphi _2=h^0h^0,h^0H^0,H^0H^0,A^0A^0`$, $`b\overline{b}`$ annihilation is the only production mechanism at the tree level. As we shall see, its cross section is comparable to — and in certain areas of the MSSM parameter space even in excess of — the one of $`gg`$ fusion. For final states with an odd number of $`CP`$-even Higgs bosons, $`\varphi _1\varphi _2=h^0A^0,H^0A^0`$, $`b\overline{b}`$ annihilation significantly enhances the Drell-Yan contribution for large values of $`\mathrm{tan}\beta `$. As for $`gg`$ fusion, we reproduce the analytical and numerical results for the quark loop contributions of Ref. . Apart from obvious typographical errors, we also agree with the formulas for the squark loop contributions listed in Ref. . However, we find their numerical size to be considerably smaller than what was found in Ref. .

As for $`b\overline{b}`$ annihilation, it should be noted that the treatment of bottom as an active flavor inside the colliding hadrons leads to an effective description, which comprises contributions from the higher-order subprocesses $`gb\varphi _1\varphi _2b`$, $`g\overline{b}\varphi _1\varphi _2\overline{b}`$, and $`gg\varphi _1\varphi _2b\overline{b}`$. If all these subprocesses are to be explicitly included along with $`b\overline{b}\varphi _1\varphi _2`$, then it is necessary to employ a judiciously subtracted parton density function (PDF) for the $`b`$ quark in order to avoid double counting . The evaluation of $`b\overline{b}\varphi _1\varphi _2`$ with an unsubtracted $`b`$-quark PDF is expected to slightly overestimate the true cross section . For simplicity, we shall nevertheless adopt this effective approach in our analysis, keeping in mind that a QCD-correction factor below unity is to be applied.

In order to reduce the number of unknown supersymmetric input parameters, we adopt a scenario where the MSSM is embedded in a grand unified theory (GUT) involving supergravity (SUGRA) . The MSSM thus constrained is characterized by the following parameters at the GUT scale, which come in addition to $`\mathrm{tan}\beta `$ and $`m_{A^0}`$: the universal scalar mass $`m_0`$, the universal gaugino mass $`m_{1/2}`$, the trilinear Higgs-sfermion coupling $`A`$, the bilinear Higgs coupling $`B`$, and the Higgs-higgsino mass parameter $`\mu `$. Notice that $`m_{A^0}`$ is then not an independent parameter anymore, but it is fixed through the renormalization group equation. The number of parameters can be further reduced by making additional assumptions. Unification of the $`\tau `$-lepton and $`b`$-quark Yukawa couplings at the GUT scale leads to a correlation between $`m_t`$ and $`\mathrm{tan}\beta `$. Furthermore, if the electroweak symmetry is broken radiatively, then $`B`$ and $`\mu `$ are determined up to the sign of $`\mu `$. Finally, it turns out that the MSSM parameters are nearly independent of the value of $`A`$, as long as $`|A|\mathrm{\Gamma }<\mathrm{\hspace{0.17em}500}`$ GeV at the GUT scale.

This paper is organized as follows. In Sec. 2, we list analytic results for the tree-level cross sections of $`q\overline{q}\varphi _1\varphi _2`$, including the Yukawa-enhanced contributions for $`q=b`$, and the squark loop contributions to the $`gg\varphi _1\varphi _2`$ amplitudes in the MSSM. The relevant MSSM coupling constants and the squark loop form factors are relegated to Appendices A and B, respectively. In Sec. 3, we present quantitative predictions for the inclusive cross section of $`pp\varphi _1\varphi _2+X`$ at the LHC adopting a favorable SUGRA-inspired MSSM scenario. Sec. 4 contains our conclusions.

## 2 Analytic Results

In this section, we present the tree-level cross sections of the partonic subprocesses $`q\overline{q}\varphi _1\varphi _2`$, where $`\varphi _i=h^0,H^0,A^0`$, and the transition ($`T`$) matrix elements of $`gg\varphi _1\varphi _2`$ arising from squark triangle and box diagrams.

We work in the parton model of QCD with $`n_f=5`$ active quark flavors, which we take to be massless. However, we retain the $`b`$-quark Yukawa couplings at their finite values, in order not to suppress possibly sizeable contributions. We adopt the MSSM Feynman rules from Ref. . In Appendix A, we list the trilinear self-couplings of the neutral Higgs bosons and their couplings to gauge bosons and quarks. For each quark flavor $`q`$ there is a corresponding squark flavor $`\stackrel{~}{q}`$, which comes in two mass eigenstates $`i=1,2`$. The masses $`m_{\stackrel{~}{q}_i}`$ of the squarks and their trilinear couplings to the $`h^0`$ and $`H^0`$ bosons are listed in Eq. (A.5) of Ref. 111In Ref. , $`m_{\stackrel{~}{q}_i}`$ is called $`M_{\stackrel{~}{Q}a}`$. and Eq. (A.2) of Ref. , respectively. Their trilinear couplings to the $`A^0`$ boson and their quartic couplings to the $`h^0`$, $`H^0`$, and $`A^0`$ bosons may be found in Appendix A.

Considering the generic partonic subprocess $`ab\varphi _1\varphi _2`$, we denote the four-momenta of the incoming partons, $`a`$ and $`b`$, and the outgoing Higgs bosons, $`\varphi _1`$ and $`\varphi _2`$, by $`p_a`$, $`p_b`$, $`p_1`$, and $`p_2`$, respectively, and define the partonic Mandelstam variables as $`s=(p_a+p_b)^2`$, $`t=(p_ap_1)^2`$, and $`u=(p_bp_1)^2`$. The on-shell conditions read $`p_a^2=p_b^2=0`$ and $`p_i^2=h_i`$, where $`h_i`$ denotes the square of the $`\varphi _i`$-boson mass. Four-momentum conservation implies that $`s+t+u=h_1+h_2`$. Furthermore, we have $`sp_T^2=tuh_1h_2`$, where $`p_T`$ is the absolute value of transverse momentum common to $`\varphi _1`$ and $`\varphi _2`$ in the center-of-mass (c.m.) frame.

The tree-level diagrams for $`b\overline{b}\varphi _1\varphi _2`$, with $`\varphi _1\varphi _2=h^0h^0,h^0H^0,H^0H^0,A^0A^0`$ and $`\varphi _1\varphi _2=h^0A^0,H^0A^0`$, are depicted in Figs. 1(a) and (b), respectively. The cross sections for the first class of partonic subprocesses may be generically written as

$$\frac{d\sigma }{dt}(b\overline{b}\varphi _1\varphi _2)=\frac{1}{1+\delta _{\varphi _1\varphi _2}}\frac{G_F^2m_W^4}{3\pi s}\left(|S|^2+p_T^2T_{}^2\right),$$
(1)
where the prefactor accounts for identical-particle symmetrization if $`\varphi _1=\varphi _2`$, $`G_F`$ is Fermi’s constant, $`m_W`$ is the $`W`$-boson mass,

$`S`$ $`=`$ $`g_{\varphi _1\varphi _2h^0}g_{h^0bb}𝒫_{h^0}(s)+g_{\varphi _1\varphi _2H^0}g_{H^0bb}𝒫_{H^0}(s),`$
$`T_\pm `$ $`=`$ $`g_{\varphi _1bb}g_{\varphi _2bb}\left({\displaystyle \frac{1}{t}}\pm {\displaystyle \frac{1}{u}}\right).`$ (2)
Here,

$$𝒫_X(s)=\frac{1}{sm_X^2+im_X\mathrm{\Gamma }_X}$$
(3)
is the propagator function of particle $`X`$, with mass $`m_X`$ and total decay width $`\mathrm{\Gamma }_X`$. For the second class of partonic subprocesses, we have

$$\frac{d\sigma }{dt}(b\overline{b}\varphi _1\varphi _2)=\frac{G_F^2m_W^4}{3\pi s}\left[|P|^2+p_T^2(|V|^2+|AT_+|^2)\right],$$
(4)
where

$`P`$ $`=`$ $`g_{\varphi _1\varphi _2A^0}g_{A^0bb}𝒫_{A^0}(s),`$
$`V`$ $`=`$ $`2g_{\varphi _1\varphi _2Z}v_{Zbb}𝒫_Z(s),`$
$`A`$ $`=`$ $`2g_{\varphi _1\varphi _2Z}a_{Zbb}𝒫_Z(s).`$ (5)
Here, $`v_{Zbb}=\left(I_b2s_w^2Q_b\right)/(2c_w)`$ and $`a_{Zbb}=I_b/(2c_w)`$, with $`c_w^2=1s_w^2=m_W^2/m_Z^2`$, are the vector and axial-vector couplings of the $`b`$ quark, with weak isospin $`I_b=1/2`$ and electric charge $`Q_b=1/3`$, to the $`Z`$ boson. As for $`h^0A^0`$ and $`H^0A^0`$ production, there are also sizeable contributions from $`q\overline{q}`$ annihilation via a $`Z`$ boson for the quarks of the first and second generations, $`q=u,d,s,c`$. The corresponding Drell-Yan cross sections are obtained from Eq. (4) by putting $`P=T_+=0`$ and substituting $`bq`$. The resulting expression agrees with Eq. (36) of Ref. . The full tree-level cross sections are then obtained by complementing the $`b\overline{b}`$-initiated cross sections of Eq. (4) with the Drell-Yan cross sections for $`q=u,d,s,c`$.

The one-loop diagrams for $`gg\varphi _1\varphi _2`$, with $`\varphi _1\varphi _2=h^0h^0,h^0H^0,H^0H^0,A^0A^0`$ and $`\varphi _1\varphi _2=h^0A^0,H^0A^0`$, are depicted in Figs. 2(a) and (b), respectively. As for the quark loops, our analytical results fully agree with those listed in Ref. , and there is no need to repeat them here. For the partonic subprocesses of class two, the squark loop contributions are zero . This may be understood as follows. (i) The $`g_{g\stackrel{~}{q}_i\stackrel{~}{q}_j}`$ and $`g_{Z\stackrel{~}{q}_i\stackrel{~}{q}_j}`$ couplings are linear in the squark four-momenta, while the $`g_{gg\stackrel{~}{q}_i\stackrel{~}{q}_j}`$ couplings are momentum independent. Thus, the diagrams in the third line of Fig. 2(b) each vanish upon adding their counterparts with the loop-momentum flows reversed. (ii) The $`g_{g\stackrel{~}{q}_i\stackrel{~}{q}_j}`$, $`g_{gg\stackrel{~}{q}_i\stackrel{~}{q}_j}`$, $`g_{h^0\stackrel{~}{q}_i\stackrel{~}{q}_j}`$, and $`g_{H^0\stackrel{~}{q}_i\stackrel{~}{q}_j}`$ couplings are symmetric in $`i`$ and $`j`$, while the $`g_{A^0\stackrel{~}{q}_i\stackrel{~}{q}_j}`$ coupling is antisymmetric. Thus, the diagrams in the last line of Fig. 2(b) vanish upon summation over $`i`$ and $`j`$. For the partonic subprocesses of class one, the $`T`$-matrix elements corresponding to the squark triangle and box diagrams are found to be

$`\stackrel{~}{𝒯}_{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{G_Fm_W^2}{\sqrt{2}}}{\displaystyle \frac{\alpha _s(\mu _r)}{\pi }}\epsilon _\mu ^c(p_a)\epsilon _\nu ^c(p_b)A_1^{\mu \nu }\stackrel{~}{F}_{\mathrm{}},`$
$`\stackrel{~}{𝒯}_{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{G_Fm_W^2}{\sqrt{2}}}{\displaystyle \frac{\alpha _s(\mu _r)}{\pi }}\epsilon _\mu ^c(p_a)\epsilon _\nu ^c(p_b)\left(A_1^{\mu \nu }\stackrel{~}{F}_{\mathrm{}}+A_2^{\mu \nu }\stackrel{~}{G}_{\mathrm{}}\right),`$ (6)
respectively, where $`\alpha _s(\mu _r)`$ is the strong-coupling constant at renormalization scale $`\mu _r`$, $`\epsilon _\mu ^c(p_a)`$ is the polarization four-vector of gluon $`a`$ and similarly for gluon $`b`$, it is summed over the color index $`c=1,\mathrm{},8`$,

$`A_1^{\mu \nu }`$ $`=`$ $`g^{\mu \nu }{\displaystyle \frac{2}{s}}p_b^\mu p_a^\nu ,`$
$`A_2^{\mu \nu }`$ $`=`$ $`g^{\mu \nu }+{\displaystyle \frac{2}{p_T^2}}\left({\displaystyle \frac{h_1}{s}}p_b^\mu p_a^\nu +{\displaystyle \frac{uh_1}{s}}p_1^\mu p_a^\nu +{\displaystyle \frac{th_1}{s}}p_b^\mu p_1^\nu +p_1^\mu p_1^\nu \right),`$ (7)
and the form factors $`\stackrel{~}{F}_{\mathrm{}}`$, $`\stackrel{~}{F}_{\mathrm{}}`$, and $`\stackrel{~}{G}_{\mathrm{}}`$ are listed in Appendix B. Due to Bose symmetry, $`\stackrel{~}{𝒯}_{\mathrm{}}`$ and $`\stackrel{~}{𝒯}_{\mathrm{}}`$ are invariant under the simultaneous replacements $`\mu \nu `$ and $`p_ap_b`$. Consequently, $`\stackrel{~}{F}_{\mathrm{}}`$, $`\stackrel{~}{F}_{\mathrm{}}`$, and $`\stackrel{~}{G}_{\mathrm{}}`$ are symmetric in $`t`$ and $`u`$. Our analytic results for the squark loop contributions agree with those given in Eqs. (8)–(10) of Ref. , which are expressed in terms of helicity amplitudes.222There are two obvious typographical errors on the right-hand side of Eq. (10e) in Ref. : There should be an overall minus sign, and $`V_{H_{(i,j)}\stackrel{~}{q}_k\stackrel{~}{q}_k}`$ should be replaced by $`V_{H_{(i,j)}\stackrel{~}{q}_k\stackrel{~}{q}_l}`$.

The parton-level cross section of $`gg\varphi _1\varphi _2`$ including both quark and squark contributions is then given by

$`{\displaystyle \frac{d\sigma }{dt}}(gg\varphi _1\varphi _2)`$ $`=`$ $`{\displaystyle \frac{1}{1+\delta _{\varphi _1\varphi _2}}}{\displaystyle \frac{G_F^2\alpha _s^2(\mu _r)}{256(2\pi )^3}}\left[\right|{\displaystyle \underset{Q=t,b}{}}C_{\mathrm{}}^QF_{\mathrm{}}^Q+F_{\mathrm{}}{\displaystyle \frac{2m_W^2}{s}}(\stackrel{~}{F}_{\mathrm{}}+\stackrel{~}{F}_{\mathrm{}})|^2`$ (8)
$`+|G_{\mathrm{}}{\displaystyle \frac{2m_W^2}{s}}\stackrel{~}{G}_{\mathrm{}}|^2+|H_{\mathrm{}}|^2],`$
where the generalized couplings $`C_{\mathrm{}}^Q`$ and $`C_{\mathrm{}}^Q`$ and the form factors $`F_{\mathrm{}}^Q`$, $`F_{\mathrm{}}`$, and $`G_{\mathrm{}}`$ may be found in Eq. (16)–(18) and Appendix A of Ref. , respectively.

The kinematics of the inclusive reaction $`ABCD+X`$, where $`A`$ and $`B`$ are hadrons, which are taken to be massless, and $`C`$ and $`D`$ are massive particles, is described in Sec. II of Ref. . Its double-differential cross section $`d^2\sigma /dydp_T`$, where $`y`$ and $`p_T`$ are the rapidity and transverse momentum of particle $`C`$ in the c.m. frame of the hadronic collision, may be evaluated from Eq. (2.1) of Ref. .

## 3 Phenomenological Implications

We are now in a position to explore the phenomenological implications of our results. The SM input parameters for our numerical analysis are taken to be $`G_F=1.16639\times 10^5`$ GeV-2, $`m_W=80.419`$ GeV, $`m_Z=91.1882`$ GeV, $`m_t=174.3`$ GeV , and $`m_b=4.6`$ GeV . We adopt the lowest-order (LO) proton PDF set CTEQ5L . We evaluate $`\alpha _s(\mu _r)`$ from the LO formula  with $`n_f=5`$ quark flavors and asymptotic scale parameter $`\mathrm{\Lambda }_{\mathrm{QCD}}^{(5)}=146`$ MeV . We identify the renormalization and factorization scales with the $`\varphi _1\varphi _2`$ invariant mass $`\sqrt{s}`$, $`M=\mu _r=\sqrt{s}`$. We vary $`\mathrm{tan}\beta `$ and $`m_{A^0}`$ in the ranges $`3<\mathrm{tan}\beta <38m_t/m_b`$ and 90 GeV$`<m_{A^0}<1`$ TeV, respectively. As for the GUT parameters, we choose $`m_{1/2}=150`$ GeV, $`A=0`$, and $`\mu <0`$, and tune $`m_0`$ so as to be consistent with the desired value of $`m_{A^0}`$. All other MSSM parameters are then determined according to the SUGRA-inspired scenario as implemented in the program package SUSPECT . For the typical example of $`\mathrm{tan}\beta =3`$ and $`m_{A^0}=300`$ GeV, the residual masses and total decay widths of the $`\varphi _i`$ bosons are $`m_{h^0}=90`$ GeV, $`m_{H^0}=306`$ GeV, $`\mathrm{\Gamma }_{h^0}=3`$ MeV, $`\mathrm{\Gamma }_{H^0}=186`$ MeV, and $`\mathrm{\Gamma }_{A^0}=72`$ MeV, and the squark masses are $`m_{\stackrel{~}{u}_1}=m_{\stackrel{~}{c}_1}=412`$ GeV, $`m_{\stackrel{~}{u}_2}=m_{\stackrel{~}{c}_2}=422`$ GeV, $`m_{\stackrel{~}{d}_1}=m_{\stackrel{~}{s}_1}=413`$ GeV, $`m_{\stackrel{~}{d}_2}=m_{\stackrel{~}{s}_2}=428`$ GeV, $`m_{\stackrel{~}{t}_1}=317`$ GeV, $`m_{\stackrel{~}{t}_2}=443`$ GeV, $`m_{\stackrel{~}{b}_1}=384`$ GeV, and $`m_{\stackrel{~}{b}_2}=413`$ GeV, We do not impose the unification of the $`\tau `$-lepton and $`b`$-quark Yukawa couplings at the GUT scale, which would just constrain the allowed $`\mathrm{tan}\beta `$ range without any visible effect on the results for these values of $`\mathrm{tan}\beta `$. We exclude solutions which do not comply with the present experimental lower mass bounds of the sfermions, charginos, neutralinos, and Higgs bosons . In our analysis, an $`s`$-channel resonance only occurs in the process $`pph^0h^0+X`$ if $`m_{H^0}>2m_{h^0}`$.

We now study the fully integrated cross sections of $`pp\varphi _1\varphi _2+X`$ at the LHC, with c.m. energy $`\sqrt{S}=14`$ TeV. Figures 38 refer to the cases $`\varphi _1\varphi _2=h^0h^0,h^0H^0,H^0H^0,A^0A^0,`$ $`h^0A^0,H^0A^0`$, respectively. In part (a) of each figure, the $`m_{A^0}`$ dependence is studied for $`\mathrm{tan}\beta =3`$ and 30 while, in part (b), the $`\mathrm{tan}\beta `$ dependence is studied for $`m_{A^0}=300`$ GeV. We note that the SUGRA-inspired MSSM with our choice of input parameters does not permit $`\mathrm{tan}\beta `$ and $`m_{A^0}`$ to be simultaneously small, due to the experimental lower bound on the selectron mass . This explains why the curves for $`\mathrm{tan}\beta =3`$ in Figs. 38(a) only start at $`m_{A^0}240`$ GeV, while those for $`\mathrm{tan}\beta =30`$ already start at $`m_{A^0}90`$ GeV. On the other hand, $`\mathrm{tan}\beta `$ and $`m_{A^0}`$ cannot be simultaneously large either, due to the experimental lower bounds on the chargino and neutralino masses . For this reason, the curves for $`\mathrm{tan}\beta =30`$ in Figs. 38(a) already end at $`m_{A^0}560`$ GeV. Finally, the experimental $`m_{h^0}`$ lower bound  enforces $`\mathrm{tan}\beta \mathrm{\Gamma }>\mathrm{\hspace{0.17em}3}`$ if $`m_{A^0}=300`$ GeV, which is reflected in Figs. 38(b).

In Figs. 36, the $`b\overline{b}`$-annihilation contributions (dashed lines), which originate from Yukawa-enhanced amplitudes, and the total $`gg`$-fusion contributions (solid lines), corresponding to the coherent superposition of quark and squark loop amplitudes, are presented separately. For a comparison with future experimental data, they should be added. For comparison, also the $`gg`$-fusion contributions due to quark loops only (dotted lines) are shown. We first assess the relative importance of the $`b\overline{b}`$-annihilation and $`gg`$-fusion contributions. In the case of $`h^0h^0`$ production, $`b\overline{b}`$ annihilation is more important than $`gg`$ fusion for intermediate values of $`m_{A^0}`$, around 300 GeV, except at the edges of the allowed $`\mathrm{tan}\beta `$ range, while it is greatly suppressed for large values of $`m_{A^0}`$, independent of $`\mathrm{tan}\beta `$ (see Fig. 3). In the case of $`h^0H^0`$ production, $`b\overline{b}`$ annihilation dominates for intermediate to large values of $`m_{A^0}`$ and large values of $`\mathrm{tan}\beta `$, while it yields an insignificant contribution for small values of $`\mathrm{tan}\beta `$, independent of $`m_{A^0}`$ (see Fig. 4). As for $`H^0H^0`$ and $`A^0A^0`$ production, $`b\overline{b}`$ annihilation is suppressed compared to $`gg`$ fusion. For $`\mathrm{tan}\beta \mathrm{\Gamma }>\mathrm{\hspace{0.17em}8}`$, the suppression factor is modest, ranging between 2 and 3, but it dramatically increases as $`\mathrm{tan}\beta `$ becomes smaller (see Figs. 5 and 6). In order to avoid confusion, we should mention that the $`b\overline{b}`$-annihilation contribution for $`\mathrm{tan}\beta =3`$ is too small to be visible in Figs. 46(a). We now investigate the size of the supersymmetric corrections to $`gg`$ fusion, i.e., the effect of including the squark loops. We observe that these corrections can be of either sign and reach a magnitude of up to 90%. Specifically, for $`h^0h^0`$, $`h^0H^0`$, $`H^0H^0`$, and $`A^0A^0`$ production, they vary within the ranges $`10\%`$ to $`+3\%`$, $`16\%`$ to $`+32\%`$, $`32\%`$ to $`+24\%`$, and $`0\%`$ to $`+90\%`$, respectively, for the values of $`m_{A^0}`$ and $`\mathrm{tan}\beta `$ considered in Figs. 36(a). The fact that they are relatively modest in most cases is characteristic for our SUGRA-inspired MSSM scenario. This is partly due to the destructive interference of quark and squark loop amplitudes and to the suppression of the latter by heavy-squark propagators. By contrast, in Ref. , the supersymmetric corrections were found reach values in excess of 100. Since the authors of Ref.  did not specify all their input parameters, we could not reproduce their numerical results for the squark loop contributions.

In Figs. 7 and 8, the total $`q\overline{q}`$-annihilation contributions (dashed lines), corresponding to the coherent superposition of Drell-Yan and Yukawa-enhanced amplitudes, and the $`gg`$-fusion contributions (solid lines), which now only receive contributions from quark loops, are given separately. For comparison, also the pure Drell-Yan contributions (dotted lines) are shown. Again, we first compare the total $`q\overline{q}`$-annihilation contributions with the $`gg`$-fusion ones. In the case of $`h^0A^0`$ production, $`q\overline{q}`$ annihilation dominates for large values of $`\mathrm{tan}\beta `$, independent of $`m_{A^0}`$, while, at the lower end of the allowed $`\mathrm{tan}\beta `$ range, it is suppressed by a factor of 40 and more, depending on the value of $`m_{A^0}`$. On the other hand, in the case of $`H^0A^0`$ production, the $`q\overline{q}`$-annihilation contribution always overshoots the $`gg`$-fusion one by at least one order of magnitude. We then examine the effect of including the Yukawa-enhanced amplitudes in the evaluation of the $`q\overline{q}`$-annihilation cross section. In the case of $`h^0A^0`$ production, there is a dramatic enhancement for large values of $`\mathrm{tan}\beta `$, which may reach several orders of magnitude for large values of $`m_{A^0}`$. In the case of $`H^0A^0`$ production, there is also an enhancement for large values of $`\mathrm{tan}\beta `$, but it is much more moderate, less than a factor of three. It is interesting to observe that the Drell-Yan cross section of $`H^0A^0`$ production is fairly independent of $`\mathrm{tan}\beta `$ unless $`m_{A^0}`$ is close to its lower bound. This may be understood by observing that $`\mathrm{sin}(\alpha \beta )`$, which governs the $`H^0A^0Z`$ coupling $`g_{H^0A^0Z}`$, defined in Eq. (A.2), is then always very close to $`1`$. This is also apparent from Fig. 2 of Ref. .

At this point, we should estimate the theoretical uncertainties in our predictions. As a typical example, we consider the cross section of $`pph^0h^0+X`$ for $`\mathrm{tan}\beta =3`$ and $`m_{A^0}=300`$ GeV. In order to obtain a hint on the size of the as-yet unknown next-to-leading-order (NLO) corrections, we define the renormalization and factorization scales as $`M=\mu _r=\xi \sqrt{s}`$ and vary the scale parameter $`\xi `$ in the range $`1/2<\xi <2`$. The resulting variation in cross section amounts to $`\pm 8\%`$ in the case of $`b\overline{b}`$ annihilation and to $`\pm 11\%`$ in the case of $`gg`$ fusion. At NLO, one also needs to specify a renormalization scheme for the definition of the $`b`$-quark mass, which enters our analysis through the $`b`$-quark Yukawa coupling. Our LO analysis is appropriate for the on-mass-shell scheme, which uses pole masses as basic parameters. The modified minimal-subtraction ($`\overline{\mathrm{MS}}`$) scheme  provides a popular alternative. For example, a pole mass of $`m_b=4.6`$ GeV  corresponds to an $`\overline{\mathrm{MS}}`$ mass of $`\overline{m}_b^{(5)}(\mu _r)=2.7`$ GeV for the typical choice of renormalization scale $`\mu _r=\sqrt{s}=300`$ GeV. Recalling that the leading behaviour of Eq. (1) in $`m_b`$ is quadratic, switching to the $`\overline{\mathrm{MS}}`$ scheme would thus, at first sight, lead to a suppression of the cross section by a factor of approximately 1/3. However, we must keep in mind that this reduction should be largely cancelled, up to terms that are formally beyond NLO, by a respective shift in the NLO correction. Another source of uncertainty is related to the choice of PDF’s. In fact, there exist significant differences in the extraction of the $`b`$-quark PDF among different PDF sets, which are related to the threshold treatment of the $`gb\overline{b}`$ splitting, the choice of the $`b`$-quark mass, the dependence of the evolution on the latter, etc. In the case of a LO analysis with $`n_f=5`$ massless quark favours, which is considered here, the use of CTEQ5L  together with the scale choice $`M=\mu _r=\sqrt{s}`$ should be appropriate in the sense that these issues can largely be bypassed. If we employ the LO PDF set by Martin, Roberts, Stirling, and Thorne (MRST) , with $`\mathrm{\Lambda }_{\mathrm{QCD}}^{(5)}=132`$ MeV, then the $`b\overline{b}`$-annihilation and $`gg`$-fusion cross sections increase by $`10\%`$ and $`4\%`$, respectively, relative to their default values.

## 4 Conclusions

We analytically calculated the cross sections of the partonic subprocesses $`q\overline{q}\varphi _1\varphi _2`$ and $`gg\varphi _1\varphi _2`$, where $`\varphi _i=h^0,H^0,A^0`$, to LO in the MSSM. We included the Drell-Yan and Yukawa-enhanced contributions to $`q\overline{q}`$ annihilation (see Fig. 1) and the quark and squark loop contributions to $`gg`$ fusion (see Fig. 2). We listed our formulas for the $`q\overline{q}`$-annihilation cross sections and the squark loop amplitudes, for which we found rather compact expressions. As for the quark loop contributions, we found complete agreement with Ref. .

We then quantitatively investigated the inclusive cross sections of $`pp\varphi _1\varphi _2+X`$ at the LHC adopting a favorable SUGRA-inspired MSSM scenario, varying the input parameters $`m_{A^0}`$ and $`\mathrm{tan}\beta `$. The results are presented in Figs. 38. We found that the Yukawa-enhanced $`q\overline{q}`$-annihilation contribution, which had previously been neglected, can play a leading role, especially for $`h^0h^0`$ production if $`m_{A^0}`$ is of order 300 GeV and for $`h^0H^0`$, $`h^0A^0`$, and $`H^0A^0`$ production if $`\mathrm{tan}\beta `$ is large. The supersymmetric corrections to $`gg`$ fusion, which are present for $`h^0h^0`$, $`h^0H^0`$, $`H^0H^0`$, and $`A^0A^0`$ production, can be of either sign and reach a magnitude of up to 90%. Our numerical results for these corrections disagree with those presented in Ref. . For each process $`pp\varphi _1\varphi _2+X`$, the combined cross section, i.e., the sum of the full $`q\overline{q}`$-annihilation and $`gg`$-fusion contributions, varies by several orders of magnitude as the values of $`m_{A^0}`$ and $`\mathrm{tan}\beta `$ are changed within their allowed ranges, and its maximum value is typically between $`10^2`$ fb-1 and $`10^3`$ fb-1. If we assume the integrated luminosity per year to be at its design value of $`L=100`$ fb-1 for each of the two LHC experiments, ATLAS and CMS, then this translates into a maximum of 20.000 to 200.000 events per year for each of these signal processes.

A comprehensive discussion of the background processes competing with the $`pp\varphi _1\varphi _2+X`$ signals at the LHC lies beyond the scope of our study. However, we should briefly mention them and quote the relevant literature. Without specifying the decay channels of the $`\varphi _i`$ bosons, one expects the major backgrounds to arise from the pair production of neutral gauge bosons, the associate production of a neutral gauge boson and a neutral Higgs boson, and the continuum production of the respective $`\varphi _1\varphi _2`$ decay products. Published signal-to-background analyses  have concentrated on the $`\varphi _1\varphi _2b\overline{b}b\overline{b}`$ signals and their irreducible continuum backgrounds, due to the partonic subprocesses $`gg,q\overline{q}b\overline{b}b\overline{b}`$, which are dominantly of pure QCD origin. It has been demonstrated that, after optimizing the acceptance cuts, the LHC experiments might discover a signal, with experimental significance in excess of 5, if $`\mathrm{tan}\beta \mathrm{\Gamma }<\mathrm{\hspace{0.17em}3}`$ or $`\mathrm{tan}\beta \mathrm{\Gamma }>\mathrm{\hspace{0.17em}50}`$.

Acknowledgements

We thank S. Moretti for drawing our attention to Ref. . The work of A.A.B.B. was supported by the Deutsches Elektronen-Synchrotron DESY. The work of B.A.K. was supported in part by the Deutsche Forschungsgemeinschaft through Grant No. KN 365/1-1, by the Bundesministerium für Bildung und Forschung through Grant No. 05 HT9GUA 3, and by the European Commission through the Research Training Network Quantum Chromodynamics and the Deep Structure of Elementary Particles under Contract No. ERBFMRX-CT98-0194.

## Appendix A Relevant Higgs and squark couplings

In this appendix, we list the trilinear self-couplings of the $`h^0`$, $`H^0`$, and $`A^0`$ bosons as well as their couplings to the $`Z`$ boson and the $`t`$ and $`b`$ quarks. Furthermore, we collect the couplings of these Higgs bosons to the squarks $`\stackrel{~}{q}_i`$, with $`q=t,b`$ and $`i=1,2`$, which are not contained in Appendix A of Ref. . For convenience, we introduce the short-hand notations $`s_\alpha =\mathrm{sin}\alpha `$, $`c_\alpha =\mathrm{cos}\alpha `$, $`s_\beta =\mathrm{sin}\beta `$, $`c_\beta =\mathrm{cos}\beta `$, $`s_{2\beta }=\mathrm{sin}(2\beta )`$, $`c_{2\beta }=\mathrm{cos}(2\beta )`$, $`s_\pm =\mathrm{sin}(\alpha \pm \beta )`$, and $`c_\pm =\mathrm{cos}(\alpha \pm \beta )`$.

The trilinear self-couplings of the $`h^0`$, $`H^0`$, and $`A^0`$ bosons are given by

$`g_{h^0h^0h^0}`$ $`=`$ $`{\displaystyle \frac{3m_Z}{2c_w}}c_{2\alpha }s_+,g_{h^0h^0H^0}={\displaystyle \frac{m_Z}{2c_w}}(2s_{2\alpha }s_+c_{2\alpha }c_+),`$
$`g_{h^0H^0H^0}`$ $`=`$ $`{\displaystyle \frac{m_Z}{2c_w}}(2s_{2\alpha }c_++c_{2\alpha }s_+),g_{H^0H^0H^0}={\displaystyle \frac{3m_Z}{2c_w}}c_{2\alpha }c_+,`$
$`g_{h^0A^0A^0}`$ $`=`$ $`{\displaystyle \frac{m_Z}{2c_w}}c_{2\beta }s_+,g_{H^0A^0A^0}={\displaystyle \frac{m_Z}{2c_w}}c_{2\beta }c_+.`$ (A.1)
Their couplings to the $`Z`$ boson are given by

$$g_{h^0A^0Z}=\frac{c_{}}{2c_w},g_{H^0A^0Z}=\frac{s_{}}{2c_w}.$$
(A.2)
Their couplings to the $`t`$ and $`b`$ quarks are given by

$`g_{h^0tt}`$ $`=`$ $`{\displaystyle \frac{m_tc_\alpha }{2m_Ws_\beta }},g_{h^0bb}={\displaystyle \frac{m_bs_\alpha }{2m_Wc_\beta }},`$
$`g_{H^0tt}`$ $`=`$ $`{\displaystyle \frac{m_ts_\alpha }{2m_Ws_\beta }},g_{H^0bb}={\displaystyle \frac{m_bc_\alpha }{2m_Wc_\beta }},`$
$`g_{A^0tt}`$ $`=`$ $`{\displaystyle \frac{m_t\mathrm{cot}\beta }{2m_W}},g_{A^0bb}={\displaystyle \frac{m_b\mathrm{tan}\beta }{2m_W}}.`$ (A.3)
The missing couplings of these Higgs bosons to the squarks are given by

$`\left(\begin{array}{cc}g_{A^0\stackrel{~}{t}_1\stackrel{~}{t}_1}& g_{A^0\stackrel{~}{t}_1\stackrel{~}{t}_2}\\ g_{A^0\stackrel{~}{t}_2\stackrel{~}{t}_1}& g_{A^0\stackrel{~}{t}_2\stackrel{~}{t}_2}\end{array}\right)`$ $`=`$ $`{\displaystyle \frac{m_t(\mu +A_t\mathrm{cot}\beta )}{2m_W}}\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),`$ (A.8)
$`\left(\begin{array}{cc}g_{A^0\stackrel{~}{b}_1\stackrel{~}{b}_1}& g_{A^0\stackrel{~}{b}_1\stackrel{~}{b}_2}\\ g_{A^0\stackrel{~}{b}_2\stackrel{~}{b}_1}& g_{A^0\stackrel{~}{b}_2\stackrel{~}{b}_2}\end{array}\right)`$ $`=`$ $`{\displaystyle \frac{m_b(\mu +A_b\mathrm{tan}\beta )}{2m_W}}\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),`$ (A.13)
$`\left(\begin{array}{cc}g_{h^0h^0\stackrel{~}{t}_1\stackrel{~}{t}_1}& g_{h^0h^0\stackrel{~}{t}_1\stackrel{~}{t}_2}\\ g_{h^0h^0\stackrel{~}{t}_2\stackrel{~}{t}_1}& g_{h^0h^0\stackrel{~}{t}_2\stackrel{~}{t}_2}\end{array}\right)`$ $`=`$ $`^{\stackrel{~}{t}}\left(\begin{array}{cc}\frac{c_{2\alpha }\left(I_ts_w^2Q_t\right)}{2c_w^2}\frac{m_t^2c_\alpha ^2}{2m_W^2s_\beta ^2}& 0\\ 0& \frac{c_{2\alpha }s_w^2Q_t}{2c_w^2}\frac{m_t^2c_\alpha ^2}{2m_W^2s_\beta ^2}\end{array}\right)\left(^{\stackrel{~}{t}}\right)^T,`$ (A.18)
$`\left(\begin{array}{cc}g_{h^0h^0\stackrel{~}{b}_1\stackrel{~}{b}_1}& g_{h^0h^0\stackrel{~}{b}_1\stackrel{~}{b}_2}\\ g_{h^0h^0\stackrel{~}{b}_2\stackrel{~}{b}_1}& g_{h^0h^0\stackrel{~}{b}_2\stackrel{~}{b}_2}\end{array}\right)`$ $`=`$ $`^{\stackrel{~}{b}}\left(\begin{array}{cc}\frac{c_{2\alpha }(I_bs_w^2Q_b)}{2c_w^2}\frac{m_b^2s_\alpha ^2}{2m_W^2c_\beta ^2}& 0\\ 0& \frac{c_{2\alpha }s_w^2Q_b}{2c_w^2}\frac{m_b^2s_\alpha ^2}{2m_W^2c_\beta ^2}\end{array}\right)\left(^{\stackrel{~}{b}}\right)^T,`$ (A.23)
$`\left(\begin{array}{cc}g_{h^0H^0\stackrel{~}{t}_1\stackrel{~}{t}_1}& g_{h^0H^0\stackrel{~}{t}_1\stackrel{~}{t}_2}\\ g_{h^0H^0\stackrel{~}{t}_2\stackrel{~}{t}_1}& g_{h^0H^0\stackrel{~}{t}_2\stackrel{~}{t}_2}\end{array}\right)`$ $`=`$ $`^{\stackrel{~}{t}}\left(\begin{array}{cc}\frac{s_{2\alpha }\left(I_ts_w^2Q_t\right)}{2c_w^2}\frac{m_t^2s_{2\alpha }}{4m_W^2s_\beta ^2}& 0\\ 0& \frac{s_{2\alpha }s_w^2Q_t}{2c_w^2}\frac{m_t^2s_{2\alpha }}{4m_W^2s_\beta ^2}\end{array}\right)\left(^{\stackrel{~}{t}}\right)^T,`$ (A.28)
$`\left(\begin{array}{cc}g_{h^0H^0\stackrel{~}{b}_1\stackrel{~}{b}_1}& g_{h^0H^0\stackrel{~}{b}_1\stackrel{~}{b}_2}\\ g_{h^0H^0\stackrel{~}{b}_2\stackrel{~}{b}_1}& g_{h^0H^0\stackrel{~}{b}_2\stackrel{~}{b}_2}\end{array}\right)`$ $`=`$ $`^{\stackrel{~}{b}}\left(\begin{array}{cc}\frac{s_{2\alpha }(I_bs_w^2Q_b)}{2c_w^2}+\frac{m_b^2s_{2\alpha }}{4m_W^2c_\beta ^2}& 0\\ 0& \frac{s_{2\alpha }s_w^2Q_b}{2c_w^2}+\frac{m_b^2s_{2\alpha }}{4m_W^2c_\beta ^2}\end{array}\right)\left(^{\stackrel{~}{b}}\right)^T,`$ (A.33)
$`\left(\begin{array}{cc}g_{H^0H^0\stackrel{~}{t}_1\stackrel{~}{t}_1}& g_{H^0H^0\stackrel{~}{t}_1\stackrel{~}{t}_2}\\ g_{H^0H^0\stackrel{~}{t}_2\stackrel{~}{t}_1}& g_{H^0H^0\stackrel{~}{t}_2\stackrel{~}{t}_2}\end{array}\right)`$ $`=`$ $`^{\stackrel{~}{t}}\left(\begin{array}{cc}\frac{c_{2\alpha }\left(I_ts_w^2Q_t\right)}{2c_w^2}\frac{m_t^2s_\alpha ^2}{2m_W^2s_\beta ^2}& 0\\ 0& \frac{c_{2\alpha }s_w^2Q_t}{2c_w^2}\frac{m_t^2s_\alpha ^2}{2m_W^2s_\beta ^2}\end{array}\right)\left(^{\stackrel{~}{t}}\right)^T,`$ (A.38)
$`\left(\begin{array}{cc}g_{H^0H^0\stackrel{~}{b}_1\stackrel{~}{b}_1}& g_{H^0H^0\stackrel{~}{b}_1\stackrel{~}{b}_2}\\ g_{H^0H^0\stackrel{~}{b}_2\stackrel{~}{b}_1}& g_{H^0H^0\stackrel{~}{b}_2\stackrel{~}{b}_2}\end{array}\right)`$ $`=`$ $`^{\stackrel{~}{b}}\left(\begin{array}{cc}\frac{c_{2\alpha }(I_bs_w^2Q_b)}{2c_w^2}\frac{m_b^2c_\alpha ^2}{2m_W^2c_\beta ^2}& 0\\ 0& \frac{c_{2\alpha }s_w^2Q_b}{2c_w^2}\frac{m_b^2c_\alpha ^2}{2m_W^2c_\beta ^2}\end{array}\right)\left(^{\stackrel{~}{b}}\right)^T,`$ (A.43)
$`\left(\begin{array}{cc}g_{A^0A^0\stackrel{~}{t}_1\stackrel{~}{t}_1}& g_{A^0A^0\stackrel{~}{t}_1\stackrel{~}{t}_2}\\ g_{A^0A^0\stackrel{~}{t}_2\stackrel{~}{t}_1}& g_{A^0A^0\stackrel{~}{t}_2\stackrel{~}{t}_2}\end{array}\right)`$ $`=`$ $`^{\stackrel{~}{t}}\left(\begin{array}{cc}\frac{c_{2\beta }\left(I_ts_w^2Q_t\right)}{2c_w^2}\frac{m_t^2\mathrm{cot}^2\beta }{2m_W^2}& 0\\ 0& \frac{c_{2\beta }s_w^2Q_t}{2c_w^2}\frac{m_t^2\mathrm{cot}^2\beta }{2m_W^2}\end{array}\right)\left(^{\stackrel{~}{t}}\right)^T,`$ (A.48)
$`\left(\begin{array}{cc}g_{A^0A^0\stackrel{~}{b}_1\stackrel{~}{b}_1}& g_{A^0A^0\stackrel{~}{b}_1\stackrel{~}{b}_2}\\ g_{A^0A^0\stackrel{~}{b}_2\stackrel{~}{b}_1}& g_{A^0A^0\stackrel{~}{b}_2\stackrel{~}{b}_2}\end{array}\right)`$ $`=`$ $`^{\stackrel{~}{b}}\left(\begin{array}{cc}\frac{c_{2\beta }(I_bs_w^2Q_b)}{2c_w^2}\frac{m_b^2\mathrm{tan}^2\beta }{2m_W^2}& 0\\ 0& \frac{c_{2\beta }s_w^2Q_b}{2c_w^2}\frac{m_b^2\mathrm{tan}^2\beta }{2m_W^2}\end{array}\right)\left(^{\stackrel{~}{b}}\right)^T.`$ (A.53)
Here, $`^{\stackrel{~}{q}}`$ denotes the mixing matrix which rotates the left- and right-handed squark fields, $`\stackrel{~}{q}_L`$ and $`\stackrel{~}{q}_R`$, into the mass eigenstates $`\stackrel{~}{q}_i`$. Its definition may be found in Eq. (A.1) of Ref. . Relations similar to Eq. (A.53) are valid for the squarks of the first and second generations, which are also included in our analysis. However, in these cases, we neglect terms which are suppressed by the smallness of the corresponding light-quark masses.

## Appendix B Squark loop form factors

In this appendix, we express the squark triangle and box form factors, $`\stackrel{~}{F}_{\mathrm{}}`$, $`\stackrel{~}{F}_{\mathrm{}}`$, and $`\stackrel{~}{G}_{\mathrm{}}`$, in terms of the standard scalar three- and four-point functions, which we abbreviate as $`C_{ijk}^{ab}(c)=C_0(a,b,c,m_{\stackrel{~}{q}_i}^2,m_{\stackrel{~}{q}_j}^2,m_{\stackrel{~}{q}_k}^2)`$ and $`D_{ijkl}^{abcd}(e,f)=D_0(a,b,c,d,e,f,m_{\stackrel{~}{q}_i}^2,m_{\stackrel{~}{q}_j}^2,m_{\stackrel{~}{q}_k}^2,m_{\stackrel{~}{q}_l}^2)`$, respectively. The definitions of the latter may be found in Eq. (5) of Ref. .

We have

$`\stackrel{~}{F}_{\mathrm{}}`$ $`=`$ $`{\displaystyle \underset{\stackrel{~}{q}}{}}{\displaystyle \underset{i=1}{\overset{2}{}}}\left(g_{\varphi _1\varphi _2h^0}g_{h^0\stackrel{~}{q}_i\stackrel{~}{q}_i}𝒫_{h^0}(s)+g_{\varphi _1\varphi _2H^0}g_{H^0\stackrel{~}{q}_i\stackrel{~}{q}_i}𝒫_{H^0}(s)g_{\varphi _1\varphi _2\stackrel{~}{q}_i\stackrel{~}{q}_i}\right)F_1(s,m_{\stackrel{~}{q}_i}^2),`$
$`\stackrel{~}{F}_{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{2}{s}}{\displaystyle \underset{\stackrel{~}{q}}{}}{\displaystyle \underset{i,j=1}{\overset{2}{}}}g_{\varphi _1\stackrel{~}{q}_i\stackrel{~}{q}_j}g_{\varphi _2\stackrel{~}{q}_j\stackrel{~}{q}_i}F_2(s,t,h_1,h_2,m_{\stackrel{~}{q}_i}^2,m_{\stackrel{~}{q}_j}^2),`$
$`\stackrel{~}{G}_{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{2}{sp_T^2}}{\displaystyle \underset{\stackrel{~}{q}}{}}{\displaystyle \underset{i,j=1}{\overset{2}{}}}g_{\varphi _1\stackrel{~}{q}_i\stackrel{~}{q}_j}g_{\varphi _2\stackrel{~}{q}_j\stackrel{~}{q}_i}F_3(s,t,h_1,h_2,m_{\stackrel{~}{q}_i}^2,m_{\stackrel{~}{q}_j}^2).`$ (B.1)
Here, we have introduced the auxiliary functions

$`F_1(s,m_{\stackrel{~}{q}_i}^2)`$ $`=`$ $`2+4m_{\stackrel{~}{q}_i}^2C_{iii}^{00}(s),`$
$`F_2(s,t,h_1,h_2,m_{\stackrel{~}{q}_i}^2,m_{\stackrel{~}{q}_j}^2)`$ $`=`$ $`t_1C_{ijj}^{h_10}(t)t_2C_{ijj}^{h_20}(t)+2sm_{\stackrel{~}{q}_i}^2D_{ijii}^{h_1h_200}(s,t)+s\left({\displaystyle \frac{p_T^2}{2}}+m_{\stackrel{~}{q}_i}^2\right)`$
$`\times D_{ijji}^{h_10h_20}(t,u)+(tu),`$
$`F_3(s,t,h_1,h_2,m_{\stackrel{~}{q}_i}^2,m_{\stackrel{~}{q}_j}^2)`$ $`=`$ $`s\left(t+m_{\stackrel{~}{q}_i}^2\right)C_{iii}^{000}(s)+sm_{\stackrel{~}{q}_i}^2C_{jjj}^{000}(s)tt_1C_{ijj}^{h_10}(t)tt_2C_{ijj}^{h_20}(t)`$ (B.2)
$`+(t^2h_1h_2)C_{iji}^{h_1h_2}(s)+\left[st^22t_1t_2m_{\stackrel{~}{q}_i}^2+2sm_{\stackrel{~}{q}_i}^2\left(m_{\stackrel{~}{q}_i}^2m_{\stackrel{~}{q}_j}^2\right)\right]`$
$`\times D_{ijii}^{h_1h_200}(s,t)2stm_{\stackrel{~}{q}_i}^2D_{jijj}^{h_1h_200}(s,t)+sm_{\stackrel{~}{q}_i}^2\left(p_T^2+m_{\stackrel{~}{q}_i}^2m_{\stackrel{~}{q}_j}^2\right)`$
$`\times D_{ijji}^{h_10h_20}(t,u)+(tu),`$
where $`t_i=th_i`$.


format: text

text:
# Effect of an inhomogeneous external magnetic field on a quantum dot quantum computer

## I Introduction

It has long been proposed that quantum mechanical computers might provide much more computing power than machines based on classical physics. However, it was the discovery of fast quantum algorithms for tasks such as prime factorization and searching disordered databases together with effective error correction schemes that motivated a frenetic search for a suitable physical system on which reliable quantum hardware could be built.

Among the many proposed quantum hardware schemes the ion trap quantum computer, based on electronic states of laser cooled trapped ions, have not only demonstrated full control of elementary one and two-qubit (quantum bit, the basic unit of quantum information) logic gates, but also recently produced 4-qubit entanglement. Although the ion trap represents the current state of the art in quantum computing architecture, its scalability prospect is unclear, since moving a large number of ions between memory and computing sectors of the ion traps without destroying quantum coherence is a difficult problem. Another proposed scheme, the liquid state Nuclear Magnetic Resonance (NMR) quantum computer, based on nuclear spin states of molecules in solutions, recently demonstrated a 5-qubit order finding algorithm. However, the NMR quantum computer also has a scalability problem, since the signal to noise ratio decreases exponentially as the size of the molecules increases. Furthermore, there has been a theoretical proof that no entanglement is present in the current ensemble-averaged NMR experiments, although its implication for NMR quantum computation has not yet been fully explored.

With the scalability problem in mind it has been suggested that advanced semiconductor technology might provide the basis for building a large (i.e. many qubits) quantum computer (QC), just as it has been proved to be effective in building today’s classical digital computers. Two of the main proposals for such semiconductor QCs are based on qubits using the spin of a phosphorous donor nucleus in bulk silicon and the spin of a single electron trapped in a GaAs quantum dot. Even though our work here refers specifically to the latter, our results are valid for both solid state QC proposals since the two are similar in the sense that they use the same universal quantum gate to perform quantum computation: The exchange or swap gate which acts on two qubits coupled by an effective Heisenberg spin interaction and one-qubit gates that perform rotations on two orthogonal spin axes by applying external magnetic fields. Experimentally demonstrating the operation of any of these quantum gates is a formidable task and has not been carried out yet, and all semiconductor based QCs are at best promising proposals at this stage. Nevertheless, the hope is that after the first step of successfully demonstrating a single qubit in a semiconductor-based QC, scaling up from one qubit to many would be feasible based on the established infrastructure of semiconductor microelectronics industry. Moreover, precisely controlling the quantum dynamics of a few spins in a semiconductor nanostructure represents a major experimental challenge, and therefore any success will lead to new insight into mesoscopic spin properties of nanostructures, which by itself is a subject of intrinsic fundamental interest.

Quantum computation can be considered as a complicated unitary evolution of a N-qubit system, governed by a N-qubit Hamiltonian, where each qubit is a quantum mechanical two-level system. It has been proved that this usually extremely complicated process can be factored into a number (that is a polynomial in N) of simple one and two-qubit processes with sufficiently high accuracy, in analogy with a classical network of logic gates. To improve efficiency multiple simple operations might be lumped together, as was proposed in various parallel pulse schemes. In these operations for a quantum dot quantum computer (QDQC) an inhomogeneous magnetic field inevitably appears, as it is a powerful technique to rotate spins while controlling spin exchange. The inhomogeneity in the field here arises necessarily from the fact that the parallel scheme employs local magnetic field pulses on individual dots to achieve different spin rotations, and a magnetic field localized on one dot is essentially an inhomogeneous field for a different dot. Furthermore, inhomogeneous magnetic fields can be produced by unwanted magnetic impurities or external currents in the sample, which will be invariably present in any real QDQC architecture. Since in the proposed parallel schemes the exchange coupling J (which controls the swap gate) is treated as independent of the magnetic field inhomogeneity, it is important to explore whether this dependence is truly weak enough, or if it is necessary to adjust the inter-dot barrier continuously in order to maintain the desired J in the presence of a varying inhomogeneous field (which will be a formidable, if not impossible, task).

We showed in a previous paper that an inhomogeneous magnetic field leads to errors in the swap operation of two qubits because the Zeeman term in the effective Hamiltonian is not proportional to the total spin. For example, if we try to swap the state $`|`$ we would obtain $`|`$ plus an entangled component proportional to the magnetic field inhomogeneity. In other words, once a product state is entered, the output will always be entangled, although for practical purposes this error can be kept within the currently estimated error correction bounds ($`10^4`$) for usual field inhomogeneities. This result is based on the assumption that one knows J exactly in order to achieve the best control, otherwise there can be additional errors in swap because of the uncertainty in J.

Motivated by the above concerns and problems arising from inhomogeneous magnetic fields, we calculate in this paper the exchange coupling J of a two-electron double dot as a function of the inhomogeneity in the external magnetic field. We employ a single envelope function approach and calculate the orbital energies of the two electrons, thus relating the general two-electron problem to a special two-spin Hamiltonian, and analyze the effect of field inhomogeneities on the orbital degrees of freedom of the double dot structure. The current work is thus complementary to our recent publication dealing with the field inhomogeneity induced errors in the swap operation.

Quite recently it has been demonstrated that QDQC can be done without any applied magnetic field as long as a price is paid in increasing the number of exchange operations: Two ($`S=\frac{1}{2},S_z=+\frac{1}{2}`$) states of three spins are used as a qubit, and the nearest-neighbor exchange interaction can provide all the required one-qubit and two-qubit operations of a quantum computer. In this scheme applications of external magnetic fields for single qubit operations is not required, as exchange gates provide both 1- and 2-qubit operations. This proposal removes the problem of having strong applied inhomogeneous fields. However, stray fields due to impurities and external currents can still be present, and must be taken into account using appropriate error corrections. In this paper we focus on clarifying the issues that arise in the simpler sequential scheme that may contain an applied inhomogeneous magnetic field as well as pointing out potential errors in all exchange QC schemes. Our results are also directly relevant to the proposed parallel pulse schemes as discussed above.

## II Model

We consider two electrons in a horizontally coupled GaAs double quantum dot structure (Fig. 1). Since these electrons are at the bottom of the conduction band the single envelope function approach provides a sufficiently complete physical picture of their energy states. Suppose there is an inhomogeneous magnetic field applied so that the electron in the left well is subject to a field which is on average $`2B_{dif}`$ lower than the one in the right. For simplicity we assume a linearly changing field profile

$$𝐁(x,z)=(\overline{B}+B_{dif}\frac{x}{a_m})\widehat{}z+B_{dif}\frac{z}{a_m}\widehat{}x,$$
(1)
where $`\overline{B}=\frac{1}{2}(B_R+B_L)`$, $`B_{dif}=\frac{1}{2}(B_RB_L)`$, and $`2a_m`$ is the separation between the dots. The complete Hamiltonian is given by

$$H=H_{orb}+H_Z,$$
(2)
$$H_{orb}=h_1+h_2+C$$
(3)
$$h_i=\frac{1}{2m}\left(𝐩_i+\frac{e}{c}𝐀\right)^2+V(𝐱_i),$$
(4)
$$C=\frac{e^2}{\epsilon }\frac{1}{\left|𝐱_1𝐱_2\right|},$$
(5)
$$H_Z=g^{}\mu _B\underset{i=1}{\overset{2}{}}𝐁(𝐱_i)𝐒_i.$$
(6)
Here $`𝐱_i`$ is the position of the $`i`$th electron in the xy plane, $`m=0.067m_e`$ is the electron effective mass in GaAs, $`\epsilon =13.1`$ is the dielectric constant, $`g^{}=0.44`$ is the GaAs gyromagnetic ratio, and $`𝐀`$ is the vector potential for the inhomogeneous field:

$$𝐀=\left(\frac{\overline{B}}{2}y\right)\widehat{}x+\left[\frac{\overline{B}}{2}x+\frac{B_{dif}}{2a_m}(x^2z^2)\right]\widehat{}y.$$
(7)
Note that since the quantum dots are fabricated using a high mobility 2D electron gas (taken to be along the xy plane), the dot dynamics is assumed to be constrained to the $`z=0`$ plane because the z-width of the 2D electron gas is typically much smaller than the QD confinement in the x-y plane. The z-dependent components in $`𝐀`$ and $`𝐁`$ do not have any effect other than to guarantee the $`\times 𝐁=0`$ condition. Hence the effect of $`B_x=B_{dif}z/a_m`$ will be neglected in our calculation (the 2D electron gas has a typical thickness of the order of a few nm and neglecting $`B_x`$ is an excellent approximation for realistic QD systems). The confining potential $`V`$ in the x-y plane is modeled in our calculations by a linear combination of gaussians characterized by the parameters $`V_0`$, $`a`$, $`V_b`$, $`l_x`$, $`l_y`$, $`l_{bx}`$ and $`l_{by}`$ which can be adjusted to control the size of the dots and switch on and off the exchange coupling:

$$V=V_0\left[\mathrm{exp}\left(\frac{(xa)^2}{l_x^2}\right)+\mathrm{exp}\left(\frac{(x+a)^2}{l_x^2}\right)\right]\mathrm{exp}\left(\frac{y^2}{l_y^2}\right)+V_b\mathrm{exp}\left(\frac{x^2}{l_{bx}^2}\right)\mathrm{exp}\left(\frac{y^2}{l_{by}^2}\right).$$
(8)
In our calculation we choose $`V_0=50meV`$ as the potential well depth and $`V_b=30meV`$as the central barrier height. In addition, $`a=15nm`$, $`l_x=`$ $`l_y=30nm`$, $`l_{bx}=15nm`$, and $`l_{by}=80nm`$, which leads to an effective barrier height (difference in energy between the top of the central barrier and the bottom of one of the wells) of 9.65meV, and a slightly shifted potential minimum $`a_m=16.41nm`$. We note that our choice of Eq. (8) as the QD confinement potential model is actually quite a general approximation, and allows considerable flexibility in our QDQC Hilbert space structure calculations.

An important feature of the Zeeman Hamiltonian (6) is that it couples coordinate $`𝐱_i`$ with spin $`𝐒_i`$ in an inhomogeneous field. This interaction introduces a coupling between the singlet and the $`S_z=0`$ triplet state, and can be well described (in the Heitler-London approximation) by an effective spin Hamiltonian

$$H_Z\gamma _RS_{1z}+\gamma _LS_{2z},\gamma _R=g^{}\mu _BB_R,\gamma _L=g^{}\mu _BB_L.$$
(9)
But note that Eq. (9) is not invariant under particle permutation $`P_{12}`$, which is a natural consequence of the fact that it has a non-vanishing matrix element between singlet and triplet states (which are respectively odd and even under a $`P_{12}`$ operation acting only in the two-spin space). This lifting of the permutation symmetry is a direct consequence of the inhomogeneous external magnetic field in this problem.

Our main objective is to study the dependence of the exchange coupling J on the inhomogeneity of the magnetic field. The parameter J describes the splitting between the lowest energy singlet and triplet states of the Hamiltonian (2), which at low energy reduces to the effective two-spin Heisenberg-type Hamiltonian that enables the operation of the quantum gates,

$$H_{eff}=J𝐒_1𝐒_2+\gamma _RS_{1z}+\gamma _LS_{2z}.$$
(10)
We use four different approximations to calculate J from the orbital Hamiltonian given in Eq. (3): Heitler-London (HL), Molecular Orbital with S states (MO-S), with up to P states (MO-SP), and with up to D states (MO-SPD). It should be noted that the Zeeman splitting term given in Eq. (9) does not affect the calculation of J in any way and it is only added to Eq. (10) after J has been determined in order to obtain the effective 2-spin exchange Hamiltonian controlling the quantum gates.

## III Methods and results

To diagonalize the Hamiltonian given in Eq. (3), we use as our one particle basis the so called Fock-Darwin energy states centered at the bottom of the left and right wells ($`x=\pm a_m`$, $`y=0`$). The Fock-Darwin QD states, briefly described in Appendix A of this paper for the sake of completeness and notational clarity, are simply the quantized states of the 2D orbital motion of the electrons in the presence of an external parabolic confinement potential and magnetic field. Fock-Darwin states have a principal quantum number n which describes a shell structure: n=0 is referred to as S states, n=1 as P, n=2 as D, and so on. Each shell is composed of n+1 levels labeled by the magnetic quantum number m=-n, -n+2, …, n-2, n. Their energy is given by $`E_{nm}^{R/L}=\mathrm{}\mathrm{\Omega }^{R/L}(n+1)+\frac{1}{2}\mathrm{}\omega _c^{R/L}m`$ where $`\mathrm{\Omega }^{R/L}=\sqrt{\omega _0^2+(\omega _c^{R/L}/2)^2}`$, $`\omega _c^{R/L}=eB_{R/L}/mc`$ is the cyclotron frequency and $`\omega _0=\sqrt{V_0/ml_x^2}`$ approximates the harmonic confinement frequency near the bottom of the dot potential. Each Fock-Darwin state has a characteristic width $`l_{0R/L}=\sqrt{\mathrm{}/m\mathrm{\Omega }^{R/L}}`$, which assumes different values depending on whether the state is in the right well subject to a field $`B_R`$ or in the left well where the field is $`B_L`$, which are different in the inhomogeneous field case we are studying. Notice that the vector potential (7) is centered in between the dots which is located $`\pm a_m`$ to the side of each well. This means that we have to use basis functions which are gauge transformed Fock-Darwin states with an additional magnetic phase $`\mathrm{exp}\left(\pm i\frac{eB_{L/R}a_my}{2\mathrm{}c}\right)`$ (sign $`+`$ and $``$ refer to the left and right wells respectively).

In what follows we use this Fock-Darwin single electron basis to form two-electron states and solve our two-particle problem in two approximations: Heitler-London which can be solved analytically and Molecular Orbital which is solved numerically.

### A Heitler-London

Heitler-London (HL) is the simplest method to calculate the exchange coupling J in two coupled dots. It consists of the singlet and triplet wave functions formed from the lowest Fock-Darwin states (n=m=0):

$$\mathrm{\Psi }_{S/T}(1,2)=\frac{1}{\sqrt{2(1\pm |S|^2)}}\left(\mathrm{\Psi }_L(1)\mathrm{\Psi }_R(2)\pm \mathrm{\Psi }_R(1)\mathrm{\Psi }_L(2)\right)\chi _{s/t},$$
(11)
$$\mathrm{\Psi }_{R/L}(𝐫)=\frac{1}{\sqrt{\pi }l_{0R/L}}\mathrm{exp}\left(\frac{(xa_m)^2+y^2}{2l_{0R/L}^2}\right)\mathrm{exp}\left(i\frac{eB_{R/L}a_my}{2\mathrm{}c}\right),$$
(12)
where the spin functions are $`\chi _s=\frac{1}{\sqrt{2}}\left(\right||)`$ and $`\chi _t`$ is one of the three triplet states $`\frac{1}{\sqrt{2}}\left(\right|+|),|,|`$. The exchange energy is then calculated as the difference between the expectation values of the Hamiltonian (3) calculated with the singlet and triplet trial wave functions:

$$J=T\left|H_{orb}\right|TS\left|H_{orb}\right|S.$$
(13)
Using this formula in an inhomogeneous field, we can obtain an exact expression for J, which is quite complicated when both $`\overline{B}`$ and $`B_{dif}`$ are nonzero, and is given in its entirety in Appendix B. Here we show the expression for J in the special situation when $`\overline{B}=0`$, which catches the main qualitative features of the inhomogeneous case:

$`J(\overline{B}=0,B_{dif})`$ $`=`$ $`{\displaystyle \frac{2S^2}{1S^4}}\{{\displaystyle \frac{\mathrm{}^2a_m^2}{ml_0^4}}+{\displaystyle \frac{1}{4}}{\displaystyle \frac{e^2}{mc^2}}B_{dif}^2(3l_0^2a_m^2){\displaystyle \frac{2V_0l_xl_y}{\sqrt{(l_x^2+l_0^2)(l_y^2+l_0^2)}}}`$ (16)
$`\times \left[\mathrm{exp}\left({\displaystyle \frac{(aa_m)^2}{l_x^2+l_0^2}}\right)+\mathrm{exp}\left({\displaystyle \frac{(a+a_m)^2}{l_x^2+l_0^2}}\right)2\mathrm{exp}\left({\displaystyle \frac{a^2}{l_x^2+l_0^2}}\right)\right]`$
$`{\displaystyle \frac{2V_bl_{bx}l_{by}}{\sqrt{(l_{bx}^2+l_0^2)(l_{by}^2+l_0^2)}}}[1\mathrm{exp}({\displaystyle \frac{a_m^2}{l_{bx}^2+l_0^2}})]\sqrt{{\displaystyle \frac{\pi }{2}}}{\displaystyle \frac{e^2}{\epsilon l_0}}[1SI_0({\displaystyle \frac{a_m^2}{l_0^2}})]\},`$
where $`S=\mathrm{exp}\left(a_m^2/l_0^2\right)`$ is the overlap, $`l_0=l_{0R}=l_{0L}`$ the Fock-Darwin characteristic width (the same for the left and right wells in the $`\overline{B}=0`$ case as $`\left|B_L\right|=\left|B_R\right|`$), and $`I_0`$ the zeroth order modified Bessel function. Note that in general $`a_ma`$ since the location of the minimum in the wells is usually different from the parameter $`a=15nm`$. The main difference between Eq. (16) and the homogeneous field $`B_{dif}=0,\overline{B}0`$ case (see Appendix B) is that the latter contains a number of phase interference terms proportional to $`\mathrm{exp}\left[(e\overline{B}/2\mathrm{}c)^2a_m^2l_0^2\right]`$, responsible for the faster decay of J as a function of $`\overline{B}`$ as will be discussed below. In Fig. 2 we depict, primarily for the purpose of comparison with the inhomogeneous field results shown in Figs. 3 and 4, the calculated exchange coupling J as a function of the external magnetic field B in the homogeneous field case with $`B_{dif}=0`$ in all our approximations.

Two curves calculated for different inhomogeneities using the HL approximation are shown in Fig. 3. The HL method is qualitatively reasonable, however, quantitatively it differs significantly from the more sophisticated molecular orbital calculations described below in section III B because J depends sensitively on the wave function overlap, which can be significantly affected by including higher excited states (as is done in MO calculations, but neglected in the HL theory). Even though the maximum field gradient obtainable in a laboratory is around 1 T in a 30nm range, we show values of $`B_{dif}`$ that go all the way up to 5 T for completeness.

### B Molecular orbital calculation

The basic idea of the MO calculation is to mix in higher excited states in the set of basis states so as to improve the energetics. To go beyond the HL approximation, we can assemble more two-particle wave functions from the Fock-Darwin states, which we call molecular orbitals

$$\mathrm{\Psi }_i(1,2)=\frac{1}{\sqrt{2}}\left(\mathrm{\Psi }_{i(1)}(1)\mathrm{\Psi }_{i(2)}(2)\pm \mathrm{\Psi }_{i(2)}(1)\mathrm{\Psi }_{i(1)}(2)\right)\chi _{s/t}.$$
(17)
Electrons are fermions so the singlet has to be symmetric in space while the triplet is antisymmetric. Here we use this molecular orbital (MO) basis to solve a generalized eigenvalue problem $`𝐇𝚿=E𝐒𝚿`$ where $`H_{ij}=\mathrm{\Psi }_i\left|H\right|\mathrm{\Psi }_j`$, $`𝚿=(\mathrm{\Psi }_1,\mathrm{\Psi }_2,\mathrm{})`$ is the vector state and the overlap matrix $`S_{ij}=<\mathrm{\Psi }_i|\mathrm{\Psi }_j>`$ is needed for orthonormalization since the MO basis is in general not orthonormal. We have done such a calculation in a constant homogeneous field, using both S and P Fock-Darwin states. Here we perform calculations with the progressively larger basis sets of S, SP, and SPD states in the inhomogeneous field problem.

In the simplest MO calculation we use only the Fock-Darwin S states to form the two-electron orbitals. Such a calculation improves upon the HL method by incorporating the doubly occupied states $`\mathrm{\Psi }_R(1)\mathrm{\Psi }_R(2)`$ and $`\mathrm{\Psi }_L(1)\mathrm{\Psi }_L(2)`$ and we call this the MO-S calculation. We then have a $`3\times 3`$ singlet matrix to diagonalize, and the triplet energy is the same as the HL case (a $`1\times 1`$ matrix). Note that we cannot use parity symmetry here to break up the matrices in a block diagonal form as was done before since the field is inhomogeneous in the x direction. The resulting MO calculation is carried out numerically. The results of such calculations are shown in Figures 2 (homogeneous case), and 3 (inhomogeneous case) where there is a noticeable quantitative difference between the HL and MO-S exchange energies.

To improve the accuracy of our molecular orbital calculations, we also perform an MO-SP calculation, adding to our basis the $`n=1`$, $`m=\pm 1`$ states on each side so that we have 6 one particle levels in total. The singlet MO-SP matrix is $`21\times 21`$ ($`\left(\begin{array}{c}6\\ 2\end{array}\right)+6=21`$) and the triplet matrix is $`15\times 15`$. We are able to efficiently calculate the Coulomb matrix elements $`C_{ij}=\mathrm{\Psi }_i\left|C\right|\mathrm{\Psi }_j`$ by reducing the 4-dimensional two-particle integrals to a 1D integral using the identity

$$\frac{1}{\left|𝐱_1𝐱_2\right|}=\frac{1}{r_{12}}=\frac{1}{\sqrt{\pi }}_0^{\mathrm{}}\frac{ds}{\sqrt{s}}\mathrm{exp}(sr_{12}^2)$$
(18)
and integrating over $`𝐱_1`$ and $`𝐱_2`$ in the gaussian basis. We then further increase the size of the MO basis by including Fock-Darwin D states (n=2, m=-2,0,2). We now have 12 single-particle levels, leading to a $`78\times 78`$ and a $`66\times 66`$ two-electron Hamiltonian matrix in the singlet and triplet basis, respectively.

Before we analyze our results it is important to determine the range of validity of our approximations. A strong inhomogeneity across the region where our basis functions vary appreciably leads to the necessity of including more and more orbitals to calculate the energy levels accurately. The extra terms that a non-vanishing inhomogeneity ($`B_{dif}0`$) adds to the one particle Hamiltonian can be written explicitly as a one particle operator $`\mathrm{\Delta }H`$. We write the vector potential as $`𝐀=𝐀_R+\delta _R`$ where $`\delta _R`$ represents the deviation from $`𝐀_R=\frac{1}{2}yB_R\widehat{}x+\frac{1}{2}xB_R\widehat{}y`$ which describes the field in the right well. The one particle Hamiltonian becomes

$$h_i=\frac{1}{2m}(𝐩_i+\frac{e}{c}𝐀_R)^2+V+\underset{\mathrm{\Delta }H}{\underset{}{\frac{e}{mc}\delta _R.(𝐩_i+\frac{e}{c}𝐀_R)+\frac{1}{2m}(\frac{e}{c})^2\delta _R^2}}.$$
(19)
$`\mathrm{\Delta }H`$ is the part of the orbital Hamiltonian (3) that is most sensitive to inhomogeneities in the magnetic field, since it is proportional to $`B_{dif}^2`$ (the Coulomb matrix elements do not depend strongly on $`B_{dif}`$). In an inhomogeneous field, $`\mathrm{\Delta }H`$ strongly couples the S orbitals with the P and D orbitals, as can be seen from Table 1, where this orbital mixing is quantitatively represented. Each individual term in $`\mathrm{\Delta }H`$ contributes approximately the same to the coupling, so we cannot neglect any of the terms in the calculation.

Note that as $`B_{dif}`$ increases, the corrections to the MO energy levels grows very fast, reaching 10% of the exchange coupling J ($`J0.2meV`$ at $`\overline{B}=0T`$ and $`J0.1meV`$ at $`\overline{B}=4T`$) for $`B_{dif}1T`$ in the second and third columns of Table 1. This suggests that the HL, MO-S and MO-SP approximations are only reliable up to $`B_{dif}=1T`$ for both the $`\overline{B}=0T`$ and $`\overline{B}=4T`$ cases. Moreover, a good estimate for the reliability of our best approximation, MO-SPD can be set as $`B_{dif}=4T`$ for the $`\overline{B}=0T`$ case and $`B_{dif}=2T`$ for the $`\overline{B}=4T`$ case. In Fig. 3 we show the exchange coupling J calculated only up to $`B_{dif}=5T`$, which is an upper limit on the reliability of all our approximations. For completeness, Fig. 4 displays the same J calculated up to $`B_{dif}=25T`$ (In this context, we should emphasize that in realistic experimental situations $`B_{dif}`$ is likely to be much less than 1 T).

On the other hand, our approximations are increasingly better at low inhomogeneity, as can be seen in Table I. Fig. 3 shows that the exchange coupling J decreases only by 0.1% in the $`\overline{B}=0T`$ case when $`B_{dif}=0.5T`$, which is the maximum field gradient obtainable in an experiment ($`B_z/x300Gauss/nm`$). Also, in the $`\overline{B}=4T`$ case it increases by 0.6% when the same field gradient is present. Using these values we can estimate that the error correction upper bound ($`10^4`$) corresponds to an inhomogeneity of $`30Gauss/nm`$. In section IV below we will argue that quantum interference in the inhomogeneous case is responsible for such a small change in J, as contrasted with the homogeneous situation in Fig. 2. Moreover, Fig. 3b shows that the exchange coupling J slightly increases as we switch on the inhomogeneity in a finite field ($`\overline{B}=4T`$). All our approximations yield an increasing J, which can be understood in a simple one-particle picture. As we increase $`B_{dif}`$ the magnetic field on the right ($`B_R`$) also increases, squeezing the Fock-Darwin state located on the right quantum dot. However, the left magnetic field ($`B_L`$) decreases to less than 4 T, which leads to a broader electronic wave function on the left quantum dot. When $`B_{dif}=4T`$, the electronic state on the left dot reaches its maximum width, and starts to be squeezed for greater inhomogeneities. This suggests a maximum of the exchange coupling J located around $`B_{dif}=4T`$. This maximum is not evident in our approximations simply because they lose their accuracy for $`B_{dif}`$ below 4 T as discussed above, but an increasing trend in J is apparent from Fig. 3b. Moreover, by looking at the $`\overline{B}=0T`$ case (Fig. 3a) we can infer that J is probably not increasing for $`B_{dif}<2T`$, since, as the MO approximations become more accurate, J first increase (MO-S and MO-P) and then decreases (MO-SPD). No maximum in J is expected here since as $`B_{dif}`$ increases both the left and right electrons are squeezed.

## IV Discussion

Comparing the results presented in Figs. 2 and 3 it becomes evident that J decreases much faster (Fig. 2) in the homogeneous magnetic field, even becoming negative at around 5T, in both HL and MO approximation schemes. This behavior can be roughly understood as arising from the approximate proportionality of J to the wave function overlap $`S=\mathrm{\Psi }_R^{}\mathrm{\Psi }_L𝑑𝐫`$ ($`\mathrm{\Psi }_{R/L}`$ are defined in Eq. (12)) multiplied by a polynomial in $`\overline{B}`$ and $`B_{dif}`$ which explains features such as the singlet-triplet crossing in Fig. 2 (the overlap S is plotted in dimensionless units together with our results for J). It turns out that $`S`$ is very sensitive to the magnetic phase $`\mathrm{exp}\left(\pm i\frac{eB_{L/R}a_my}{2\mathrm{}c}\right)`$ which appears in the right and left basis functions discussed above. In a homogeneous field ($`B_R=B_L=B`$) those phases add up, leading to destructive interference which gives an extra exponential decay proportional to $`B^2`$ in the overlap: $`S_{Hom}=\mathrm{exp}\left((a_m/l_0)^2\right)\mathrm{exp}\left(\left(\frac{1}{2}l_0a_m\frac{eB}{\mathrm{}c}\right)^2\right)`$. On the other hand in the inhomogeneous case with $`\overline{B}=0`$ , $`B_R=B_L`$ the phases cancel out in the S integral leaving just a single exponential decay, $`S_{Inh}=\mathrm{exp}\left((a_m/l_0)^2\right)`$ which depends mainly on the interdot distance. This absence of quantum interference in the inhomogeneous field case explains why J hardly changes when there is a magnetic field difference over the two quantum dots (Fig. 3), which is in fact good news for the operation of the swap gate in an environment that contains some magnetic impurities for example (or a field inhomogeneity arising from some other unintended source).

Note that in the homogeneous case the exchange coupling constant in Fig. 2 decreases extremely strongly at a few Tesla which leads to a singlet-triplet crossing. This crossing arises from the competition between the long range Coulomb interaction (which favors the triplet state) and the kinetic energy (favoring the singlet state). There is no such singlet-triplet crossing in an inhomogeneous field (except when $`\overline{B}`$ is such that $`J<0`$, where the application of $`B_{dif}`$ leads to a positive J), because the kinetic energy term (which favors the singlet state) now dominates as quantum interference terms are suppressed. In the homogeneous field case on the other hand the potential energy term always wins at high enough B fields leading to the triplet state becoming lower in energy than the singlet state.

## V Conclusion

In conclusion, we study the singlet-triplet splitting (i.e. the exchange coupling) J of a horizontally coupled double quantum dot system with two electrons. In particular, we calculate J when the magnetic field is inhomogeneous with different field strengths ($`\mathrm{\Delta }B=2B_{dif}`$) in each dot, and show that realistic values of field inhomogeneity in double quantum dot structures should have quantitatively small (but qualitatively interesting) effect on J. We perform an analytical Heitler-London calculation, and a numerical diagonalization of the orbital Hamiltonian in a set of increasingly larger two-electron MO basis to conclude that the cancellation of the magnetic phase difference between the two electrons on the dots is responsible for the relatively weak dependence of J on $`B_{dif}`$ for realistic values of $`B_{dif}`$. The swap gate can thus be operated in an inhomogeneous environment as long as the exchange coupling J is properly chosen and quantum error correction is performed to account for the inhomogeneity-induced errors arising from the Zeeman term as discussed earlier. We conclude therefore that electron spin based QDQC operations should in principle be feasible in an inhomogeneous magnetic field environment (with the inhomogeneity arising either from external magnetic fields needed for single qubit manipulation, or for running the parallel pulse scheme, or from unintended magnetic impurities or currents invariably present in the system) as long as the field inhomogeneity is not too large.

This work is supported by ARDA and the US-ONR. RdS acknowledges useful discussions with Dr. J. Fabian.

## A Fock-Darwin states

We briefly describe the Fock-Darwin states, which are used as the single particle basis in our calculations. Consider an electron of mass $`m^{}`$ restricted to move in the x-y plane and confined in a parabolic well of frequency $`\omega _0`$. A magnetic field B is applied in the z direction, which can be described by the vector potential $`𝐀=\frac{1}{2}yB\widehat{}x+\frac{1}{2}xB\widehat{}y`$. The Hamiltonian then becomes

$$H=\frac{1}{2m^{}}\left(𝐩+\frac{e}{c}𝐀\right)^2+\frac{1}{2}m^{}\omega _0^2r^2,$$
(A1)
where $`r=\sqrt{x^2+y^2}`$ is the distance between the electron and the origin, e the electronic charge and c the speed of light. Since this Hamiltonian commutes with the z component of the angular momentum operator $`L_z`$, its eigenfunctions can be chosen with well defined angular momentum $`m\mathrm{}`$. Apart from the magnetic quantum number m, a radial quantum number n completely describes the eigenenergies of the system, given by $`E_{nm}=(n+1)\mathrm{}\mathrm{\Omega }+\frac{1}{2}m\mathrm{}\omega _c`$, with $`n=0,1,2,\mathrm{}`$ and $`m=n,n+2,\mathrm{},n2,n`$. Here $`\omega _c=\frac{eB}{m^{}c}`$ is the electron’s cyclotron frequency whereas $`\mathrm{\Omega }=\sqrt{\omega _o^2+(\omega _c/2)^2}`$ is the harmonic oscillator frequency renormalized by the magnetic field. For small magnetic fields $`(\omega _c\omega _0)`$ the quantum number n describes a shell structure, with $`n=0,1,2,\mathrm{}`$ being the $`S,P,D,\mathrm{}`$ shells respectively. In table 2 below we write the angular momentum and energy eigenfunctions for the S, P and D shells that are used in calculations throughout this paper. The negative m wave functions can be easily obtained from the positive ones using the relation $`\varphi _{nm}(r,\theta )=\varphi _{nm}^{}(r,\theta )`$. Note that a characteristic length $`l_0=\sqrt{\mathrm{}/m\mathrm{\Omega }}`$ sets the range of these wave functions.

## B Exchange coupling J in the Heitler-London approximation

Here we write down our analytical results in our calculation of the exchange coupling J in the Heitler-London approximation. $`\overline{B}`$ and $`B_{dif}`$ define the inhomogeneous external magnetic field in the system (Fig. 1). As described in section III A, J is calculated by subtracting the average energies of the Heitler-London states $`\mathrm{\Psi }_{S/T}`$ (Eq. (11)) leading to

$$J=\frac{2}{1\left|S\right|^4}\left\{\left|S\right|^2\left[\mathrm{\Psi }_R|h|\mathrm{\Psi }_R+\mathrm{\Psi }_L|h|\mathrm{\Psi }_L+D\right]\left[2Re(S\mathrm{\Psi }_R|h|\mathrm{\Psi }_L)+I\right]\right\},$$
(B1)
where h is the single particle Hamiltonian (3), and $`S=\mathrm{\Psi }_L|\mathrm{\Psi }_R`$ is the overlap. We also define the direct (D) and exchange (I) Coulomb integrals

$$D=|\mathrm{\Psi }_L(1)|^2\frac{e^2/\epsilon }{|𝐱_1𝐱_2|}|\mathrm{\Psi }_R(2)|^2𝑑𝐱_1𝑑𝐱_2,$$
(B2)
$$I=\mathrm{\Psi }_L^{}(1)\mathrm{\Psi }_R^{}(2)\frac{e^2/\epsilon }{|𝐱_1𝐱_2|}\mathrm{\Psi }_L(2)\mathrm{\Psi }_R(1)𝑑𝐱_1𝑑𝐱_2,$$
(B3)
where e is the electron charge and $`\epsilon `$ the dielectric constant of the material. We list below all the analytic expressions for those quantities as a function of the effective mass m, the characteristic lengths in the right and left wells $`l_{0R/L}=\sqrt{\mathrm{}/m\mathrm{\Omega }^{R/L}}`$, the magnetic fields $`B_{R/L}=\overline{B}\pm B_{dif}`$, and the confining potential parameters $`a`$, $`l_x`$, $`l_y`$, $`l_{bx}`$, $`l_{by}`$ (Eq. (8)). Note that $`a_m`$ is the minimum location of each well, which is in general different than the parameter $`a`$. Throughout the expressions we use the zeroth order modified Bessel function $`I_0(x)`$, $`\gamma =\frac{1}{l_{0R}^2}+\frac{1}{l_{0L}^2}`$, $`\lambda =\frac{1}{l_{0R}^2}\frac{1}{l_{0L}^2}`$ , $`b=\frac{e\overline{B}}{\mathrm{}c},`$ and the overlap $`S`$, which are listed below.

$$S=\frac{2l_{0R}l_{0L}}{l_{0L}^2+l_{0R}^2}\mathrm{exp}\left(\frac{2a_m^2}{l_{0L}^2+l_{0R}^2}\frac{b^2a_m^2}{2\gamma }\right),$$
(B4)
$$D=\sqrt{\pi }\frac{e^2}{\epsilon }\frac{1}{\sqrt{l_{0L}^2+l_{0R}^2}}\mathrm{exp}\left(\frac{2a_m^2}{l_{0L}^2+l_{0R}^2}\right)I_0\left(\frac{2a_m^2}{l_{0L}^2+l_{0R}^2}\right),$$
(B5)
$$I=\sqrt{\pi }\frac{e^2}{\epsilon }\frac{S}{\sqrt{l_{0L}^2+l_{0R}^2}}\mathrm{exp}\left(\frac{2a_m^2}{l_{0L}^2+l_{0R}^2}\right)I_0\left(\frac{b^2a_m^2}{2\gamma }\right),$$
(B6)
$`\mathrm{\Psi }_L|h|\mathrm{\Psi }_L`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2ml_{0L}^2}}+{\displaystyle \frac{1}{8}}{\displaystyle \frac{e^2}{mc^2}}\left[\overline{B}^2(B_RB_L)B_L+{\displaystyle \frac{3}{4}}{\displaystyle \frac{l_{0L}^2}{a_m^2}}B_{dif}^2\right]l_{0L}^2`$ (B9)
$`{\displaystyle \frac{V_0l_xl_y}{\sqrt{\left(l_{0L}^2+l_x^2\right)\left(l_{0L}^2+l_y^2\right)}}}\left[\mathrm{exp}\left({\displaystyle \frac{(a_m+a)^2}{l_{0L}^2+l_x^2}}\right)+\mathrm{exp}\left({\displaystyle \frac{(a_ma)^2}{l_{0L}^2+l_x^2}}\right)\right]`$
$`+{\displaystyle \frac{V_bl_{bx}l_{by}}{\sqrt{\left(l_{0L}^2+l_{bx}^2\right)\left(l_{0L}^2+l_{by}^2\right)}}}\mathrm{exp}\left({\displaystyle \frac{a_m^2}{l_{0L}^2+l_{bx}^2}}\right).`$
For the average of the single particle Hamiltonian h in the right side, just substitute L for R and vice-versa in the expression above.

$`\mathrm{\Psi }_R|h|\mathrm{\Psi }_L`$ $`=`$ $`{\displaystyle \frac{1}{2\gamma }}S\{{\displaystyle \frac{\mathrm{}^2}{m}}{\displaystyle \frac{1}{l_{0L}^2l_{0R}^2}}({\displaystyle \frac{4a_m^2}{l_{0L}^2+l_{0R}^2}}2)+{\displaystyle \frac{e^2a_m^2}{mc^2}}[{\displaystyle \frac{\gamma }{4}}(\overline{B}^2+B_{dif}^2){\displaystyle \frac{\overline{B}^2}{l_{0L}^2+l_{0R}^2}}{\displaystyle \frac{\lambda }{2}}\overline{B}B_{dif}]`$ (B15)
$`+{\displaystyle \frac{2}{\gamma }}{\displaystyle \frac{e^2a_m^2}{mc^2}}\left[{\displaystyle \frac{\overline{B}^2}{l_{0L}^2l_{0R}^2}}+{\displaystyle \frac{1}{4}}\left(\lambda \overline{B}+{\displaystyle \frac{B_{dif}}{a_m^2}}+{\displaystyle \frac{\lambda ^2}{\gamma }}B_{dif}\right)\left(\lambda \overline{B}\gamma B_{dif}\right)\right]`$
$`+{\displaystyle \frac{1}{4}}{\displaystyle \frac{e^2\overline{B}^2}{mc^2}}\left[2{\displaystyle \frac{b^2a_m^2}{\gamma }}+{\displaystyle \frac{3}{\gamma a_m^2}}\left({\displaystyle \frac{B_{dif}}{\overline{B}}}\right)^2+{\displaystyle \frac{\lambda }{\gamma }}\left(1+{\displaystyle \frac{\lambda }{\gamma }}{\displaystyle \frac{B_{dif}}{\overline{B}}}\right)\left(\lambda a_m^2+{\displaystyle \frac{B_{dif}}{\overline{B}}}\left({\displaystyle \frac{\lambda ^2}{\gamma }}a_m^2+6\right)\right)\right]`$
$`{\displaystyle \frac{V_0}{\sqrt{\frac{1}{2}\left(1+\frac{l_{0L}^2}{l_{0R}^2}\right)+\frac{l_{0L}^2}{l_y^2}}\sqrt{\frac{1}{2}\left(1+\frac{l_{0R}^2}{l_{0L}^2}\right)+\frac{l_{0R}^2}{l_x^2}}}}\mathrm{exp}\left({\displaystyle \frac{b^2l_y^2a_m^2}{2\left(2+\gamma l_y^2\right)}}{\displaystyle \frac{2l_x^2a_m^2}{l_{0L}^2l_{0R}^2(2+\gamma l_x^2)}}\right)`$
$`\times \left[\mathrm{exp}\left({\displaystyle \frac{\left(a_ma\right)^2/l_{0R}^2+\left(a_m+a\right)^2/l_{0L}^2}{2+\gamma l_x^2}}\right)+\mathrm{exp}\left({\displaystyle \frac{\left(a_ma\right)^2/l_{0L}^2+\left(a_m+a\right)^2/l_{0R}^2}{2+\gamma l_x^2}}\right)\right]`$
$`+{\displaystyle \frac{V_b}{\sqrt{\frac{1}{2}\left(1+\frac{l_{0L}^2}{l_{0R}^2}\right)+\frac{l_{0L}^2}{l_{by}^2}}\sqrt{\frac{1}{2}\left(1+\frac{l_{0R}^2}{l_{0L}^2}\right)+\frac{l_{0R}^2}{l_{bx}^2}}}}\mathrm{exp}({\displaystyle \frac{b^2l_{by}^2a_m^2}{2\left(2+\gamma l_{by}^2\right)}}{\displaystyle \frac{2\frac{l_{bx}^2}{l_{0L}^2l_{0R}^2}+\gamma }{2+\gamma l_{bx}^2}}a_m^2)\}.`$
A particularly interesting expression is the homogeneous field exchange coupling $`J(\overline{B},B_{dif}=0)`$, which we list below (in addition to Eq. (16)) which is the inhomogeneous case $`J(\overline{B}=0,B_{dif})`$). Since for this case $`B_R=B_L`$, we simplify our notation by using $`l_0=l_{0R/L}`$ for the characteristic length. Also the overlap becomes $`S(\overline{B},B_{dif}=0)=\mathrm{exp}\left((a_m/l_0)^2\right)\mathrm{exp}\left(\left(\frac{1}{2}l_0ba_m\right)^2\right)`$.

$`J(\overline{B},B_{dif}=0)`$ $`=`$ $`{\displaystyle \frac{2S^2}{1S^4}}\{{\displaystyle \frac{\mathrm{}^2a_m^2}{ml_0^4}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{e^2\overline{B}^2}{mc^2}}a_m^2(1{\displaystyle \frac{1}{8}}b^2l_0^4){\displaystyle \frac{2V_0l_xl_y}{\sqrt{(l_0^2+l_x^2)(l_0^2+l_y^2)}}}`$ (B19)
$`\times \left[\mathrm{exp}\left({\displaystyle \frac{(aa_m)^2}{l_0^2+l_x^2}}\right)+\mathrm{exp}\left({\displaystyle \frac{(a+a_m)^2}{l_0^2+l_x^2}}\right)2\mathrm{exp}\left({\displaystyle \frac{a^2}{l_0^2+l_x^2}}+{\displaystyle \frac{l_0^2(bl_0a_m)^2}{4(l_0^2+l_y^2)}}\right)\right]`$
$`{\displaystyle \frac{2V_bl_{bx}l_{by}}{\sqrt{(l_0^2+l_{bx}^2)(l_0^2+l_{by}^2)}}}\left[\mathrm{exp}\left({\displaystyle \frac{l_0^2(bl_0a_m)^2}{4(l_0^2+l_{by}^2)}}\right)\mathrm{exp}\left({\displaystyle \frac{a_m^2}{l_0^2+l_{bx}^2}}\right)\right]`$
$`\sqrt{{\displaystyle \frac{\pi }{2}}}{\displaystyle \frac{e^2}{\epsilon l_0}}e^{(\frac{1}{2}bl_0a_m)^2}[I_0\left(\left({\displaystyle \frac{1}{2}}bl_0a_m\right)^2\right)SI_0\left(\left({\displaystyle \frac{a_m}{l_0}}\right)^2\right)]\}.`$