format: text

text:
# On the Problem of Programming Quantum Computers

## I Introduction

The logical operation of a conventional digital computer does not depend on the details of the hardware implementation although the speed of operation and the cost of the machine obviously do. Conventional computers are in one particular state at any time. Furthermore from the point of view of programming the computer, the internal machinery of the basic units comprising the computer is irrelevant. This is very important because it implies that on a conceptual level, algorithms designed to run on a conventional computer will produce answers that do not depend on the hardware that is used.

A quantum computer (QC) is very different in this respect. A QC exploits the fact that a quantum system can be in a superposition of states and that interference of these states allows exponentially many computations to be done in parallel . The presence of this superposition is a manifestation of the internal quantum dynamics of the elementary units (i.e. the qubits). In other words, the quantum dynamics is essential to the operation of a physically realizable QC.

The operation of an ideal QC does not depend on the intrinsic dynamics of the physical qubits: One imagines that the qubits are ideal two-state quantum systems that perform their task instantaneously and perfectly. From a theoretical point of view this situation is very similar to that of computers built from conventional digital circuits. However, in practice there is a fundamental difference. The fact that the logical operation of conventional digital circuits does not depend on their hardware implementation (e.g. semiconductors, relays, vacuum tubes, etc.) is directly related to the presence of dissipative processes that drive the circuits into regions of stable operation. Dissipation suppresses the effects of the internal, non-ideal (chaotic) dynamics of these circuits.

The quantum dynamics of small physical devices is usually very sensitive to small perturbations and this holds for qubits as well. Unfortunately, in contrast to the case of conventional circuits, dissipation usually has a devastating effect on the coherent quantum dynamical motion of the qubits, i.e. on the very essence of QC’s. Therefore the specific physical realization of a QC is intimately related to the stability of its operation.

In this paper we study the relation between the physical realization of QC’s and their logical operation and explore physical mechanisms to enlarge the region(s) of stable operation. We demonstrate that programming a physical implementation of a QC is non-trivial, even if the QC consists of only two or four qubits. In most theoretical work on QC’s and quantum algorithms (QA’s)  one considers theoretically ideal (but physically unrealizable) QC’s and therefore this problem of programming QC’s (which we will call the Quantum Programming Problem (QPP)) is not an issue. As far as we know no experimental data has been published that specifically addresses this, for potential applications, very important and intrinsic problem of programming QC’s. The aim of this paper is to investigate various aspects of the QPP by simulating QC hardware.

How does a QPP reveal itself? Consider two logically independent operations ($`O_1`$ and $`O_2`$) of the machine. On a conventional computer or ideal QC, the order in which we execute these two mutually independent instructions does not matter: $`O_1O_2=O_2O_1`$. However, it turns out that on a physically realizable QC sometimes the order does matter, even if there are no dependencies in these two program steps. In some cases $`O_1O_2O_2O_1`$ and the QC may give the wrong answers. The QPP is due to the fact that we are dealing with interacting quantum mechanical objects (as communication between qubits is essential for computation), technical difficulties to manipulate a qubit without disturbing others and the fundamental physical fact that the state of a qubit cannot be frozen during the time that other qubits are being addressed. Also note the qualifier sometimes. There seems to be no general rule to decide beforehand which operation and at what stage of the QA the QPP leads to incorrect results. At present the only way to find out seems to be to actually carry out the calculations and check the results.

It does not require a lot of imagination to see that the QPP implies that it may be very difficult to develop a non-trivial quantum program for a physical QC. Moreover there is no guarantee that a QA that works well on one QC will perform well on other QC’s. Marginal changes in the qubit hardware may affect the interchangeability of logically independent operations. There are several factors that contribute to the QPP:

* Differences between the theoretically perfect and physically realizable one- and two-qubit operations, e.g. the one-qubit operations affecting other qubits and inaccuracies on the time-interval used to perform operations.
* Physical qubits cannot be kept still while others are being addressed.
* The effect of coupling of the qubits to other degrees of freedom (dissipation, decoherence).

In this paper we address these issues through case studies. In Section II we describe the physical model that will be the starting point of our investigations. Our choice of physical models is largely inspired by NMR-QC experiments , only because other candidate technologies  for building QC’s are not yet developed to the point that they can execute computationally non-trivial QA’s.

As an example of such a QA we will take Grover’s database search algorithm (see Section III) and implement it on various physical models for 2- and 4-qubit QC’s (Sections IV to VI). Our approach for analysing the QPP is to run Grover’s QA by simulating the time-evolution of the physical model representing the QC. Thereby we strictly stick to the rules of quantum mechanics, i.e. we solve the time-dependent Schrödinger equation that describes the evolution of the physical apparatus representing the QC. The main vehicle for performing these simulations is a Quantum Computer Emulator (QCE). A detailed description of this software tool is given elsewhere . Our work is fundamentally different from those of others who also address questions related to error propagation in QA’s  in that we execute QA’s on physical models of QC hardware. The influence of non-resonant effects on the quantum computations using Ising spin quantum computers has been studied by Berman et al. . This work is similar in spirit to the one of the present paper as it explores the consequences of the difference between the ideal and physically realizable QC’s. In this paper we focus on the QPP, not on methods to suppress non-resonant effects. For simplicity most of our calculations (Sections IV and VI) are done for systems at zero temperature, in the absence of interactions with the environment. Simulations of QC’s operating at a non-zero temperature, in contact with a heat bath, are discussed in Section V.

## II Physical Model

The simplest qubit is a two-state quantum system, e.g. the spin of electrons or the polarization of photons. The basic operations in a meaningful computation are the manipulation of each qubit (e.g. by applying external fields) and the exchange of information between the qubits. In physical terms, the latter implies that there should be some interaction between the qubits. A non-trivial QC contains at least two qubits. It is known that the most simple spin-1/2 system, i.e. the Ising model, can be used for quantum computing . In the presence of an external magnetic field, the Hamiltonian of the two-spin Ising model reads

$`H=J_zS_1^zS_2^zh_z(g_1S_1^z+g_2S_2^z),`$ (1)
where $`J_z`$, $`g_1`$ and $`g_2`$ are material-specific constants, $`h_z`$ represents the applied magnetic field, and $`S_j^\alpha `$ denotes the $`\alpha `$-th component ($`\alpha =x,y,z`$) of the spin-1/2 operators describing the nuclear spins. In this paper we use units such that $`\mathrm{}=1`$.

According to the rules of quantum mechanics, the state of the QC at time $`t`$ is completely described by the wave function $`|\mathrm{\Phi }(t)`$. Executing a program on a QC is equivalent to solving the time-dependent Schrödinger equation (TDSE)

$`i{\displaystyle \frac{}{t}}|\mathrm{\Phi }(t)=H(t)|\mathrm{\Phi }(t).`$ (2)
For the 2-qubit QC the most general linear combination reads

$`|\mathrm{\Phi }(t)=a(,;t)|+a(,;t)|+a(,;t)|+a(,;t)|,`$ (3)
where the $`a(b_1,b_2;t)`$ are complex numbers. The probability that, upon measurement, the QC is in one of the four basis states $`|,\mathrm{},|`$ is given by $`|a(,;t)|^2,\mathrm{},|a(,;t)|^2`$ respectively.

The final results of a QC calculation can be read off by performing an experiment that measures the expectation value(s) of the spin(s). The value of a qubit is related to the expectation value of the $`z`$-component of the spin operator:

$`Q_j1/2\mathrm{\Phi }(t)|S_j^z|\mathrm{\Phi }(t);j=1,2,\mathrm{},N,`$ (4)
where $`N`$ denotes the total number of qubits of the QC. In this paper we will denote the state of a qubit by

$`|0|;|1|,`$ (5)
and if necessary we will add a subscript, e.g. $`|0_1`$, to label the qubit.

The Ising interaction between the spins is sufficient to implement control-NOT (CNOT) gates. Single qubit operations can be performed by applying additional external fields in the $`x`$ and $`y`$ direction (see below). It has been shown that any operation involving an arbitrary number of qubits can be written as a sequence of these elementary operations, in other words the Ising model is a universal QC.

## III Grover’s database search algorithm

Finding the needle in a haystack of $`N`$ elements requires $`𝒪(N)`$ queries on a conventional computer. Grover has shown that a QC can find the needle using only $`𝒪(\sqrt{N})`$ attempts. Assuming a uniform probability distribution for the needle, for $`N=4`$ the average number of queries required by a conventional algorithm is 9/4. With Grover’s QA the correct answer can be found in a single query. In general the reduction from $`𝒪(N)`$ to $`𝒪(\sqrt{N})`$ is due to the intrinsic massive parallel operation of the QC.

Experimentally this QA has been implemented on a 2-qubit NMR-QC for the case of a database containing four items. This implementation uses elementary rotations about 90 degrees (clock-wise) around the $`x`$ and $`y`$-axis (e.g. $`X_1|00=(|00+i|10)/\sqrt{2}`$ etc.) and an interaction-controlled phase shift. It is convenient to express these procedures in matrix notation, for example

$`X_1\left(\begin{array}{c}|00\hfill \\ |10\hfill \\ |01\hfill \\ |11\hfill \end{array}\right){\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{cccc}\hfill 1& \hfill i& \hfill 0& \hfill 0\\ \hfill i& \hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1& \hfill i\\ \hfill 0& \hfill 0& \hfill i& \hfill 1\end{array}\right)\left(\begin{array}{c}|00\hfill \\ |10\hfill \\ |01\hfill \\ |11\hfill \end{array}\right),`$ (18)
where $`|b_1b_2|b_1|b_2`$ and $`b_i=0,1`$.

$`Y_2\left(\begin{array}{c}|00\hfill \\ |10\hfill \\ |01\hfill \\ |11\hfill \end{array}\right){\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{cccc}\hfill 1& \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0& \hfill 1\\ \hfill 1& \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0& \hfill 1\end{array}\right)\left(\begin{array}{c}|00\hfill \\ |10\hfill \\ |01\hfill \\ |11\hfill \end{array}\right).`$ (31)
The two qubits “communicate” with each other through the interaction-controlled phase shift

$`I(a)\left(\begin{array}{c}|00\hfill \\ |10\hfill \\ |01\hfill \\ |11\hfill \end{array}\right){\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{cccc}\hfill e^{ia/4}& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill e^{+ia/4}& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill e^{+ia/4}& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill e^{ia/4}\end{array}\right)\left(\begin{array}{c}|00\hfill \\ |10\hfill \\ |01\hfill \\ |11\hfill \end{array}\right).`$ (44)
This unitary transformation can be used to implement a CNOT gate. Grover’s algorithm can be expressed entirely by a sequence of the basic operations (18),(31), and (44).

To see how Grover’s QA works it is instructive to consider an example of a database with $`N=4`$ positions, labeled 0,…,3. Let us assume that the item we are searching for is located at position 2. First we put the QC in its initial state $`|00`$. Then we transform this state into the uniform superposition state

$`|u{\displaystyle \frac{1}{2}}(|00+|10+|01+|11),`$ (45)
by letting the sequence $`\overline{X}_2\overline{X}_2\overline{Y}_2\overline{X}_1\overline{X}_1\overline{Y}_1`$ act on $`|00`$:

$`|u=\overline{X}_2\overline{X}_2\overline{Y}_2\overline{X}_1\overline{X}_1\overline{Y}_1|00,`$ (46)
where $`\overline{X}_j`$ ($`\overline{Y}_j`$) denotes the inverse of $`X_j`$ ($`Y_j`$).

Next we apply to $`|u`$ the sequence of elementary operations

$`F_2\overline{Y}_1X_1\overline{Y}_1Y_2\overline{X}_2\overline{Y}_2I(\pi ),`$ (47)
to encode the content of the database as

$`|\mathrm{\Psi }_2={\displaystyle \frac{1}{2}}(|00+|10|01+|11),`$ (48)
where we adopt the notation in which the basis states are labeled by the binary representation of integers with the order of the bits reversed. In (48) the position of the item in the database (i.e. 2 in this example) is encoded by modifying $`|u`$ such that the amplitude of the corresponding basis state changes sign.

The key ingredient of Grover’s algorithm is an operation that determines which of the basis state contributes to (48) with the minus-one amplitude. In matrix notation

$`G={\displaystyle \frac{1}{2}}\left(\begin{array}{cccc}\hfill 1& \hfill 1& \hfill 1& \hfill 1\\ \hfill 1& \hfill 1& \hfill 1& \hfill 1\\ \hfill 1& \hfill 1& \hfill 1& \hfill 1\\ \hfill 1& \hfill 1& \hfill 1& \hfill 1\end{array}\right),`$ (53)
and in terms of elementary operations

$`G=X_1X_1\overline{Y}_1X_2X_2\overline{Y}_2Y_1\overline{X}_1\overline{Y}_1Y_2\overline{X}_2\overline{Y}_2I(\pi )X_1X_1\overline{Y}_1X_2X_2\overline{Y}_2.`$ (54)
The sequence (54) is by no means unique: Various alternative expressions can be written down by using the algebraic properties of the $`X`$’s and $`Y`$’s. This feature can be exploited to eliminate redundant elementary operations. Starting from the uniform superposition state $`|u`$, one choice for the optimized sequences that implement the four different states of the database and Grover’s search algorithm is

$`U_0`$ $`=X_1\overline{Y}_1X_2\overline{Y}_2I(\pi )X_1\overline{Y}_1X_2\overline{Y}_2I(\pi ),`$ (55a)
$`U_1`$ $`=X_1\overline{Y}_1X_2\overline{Y}_2I(\pi )X_1\overline{Y}_1\overline{X}_2\overline{Y}_2I(\pi ),`$ (55b)
$`U_2`$ $`=X_1\overline{Y}_1X_2\overline{Y}_2I(\pi )\overline{X}_1\overline{Y}_1X_2\overline{Y}_2I(\pi ),`$ (55c)
$`U_3`$ $`=X_1\overline{Y}_1X_2\overline{Y}_2I(\pi )\overline{X}_1\overline{Y}_1\overline{X}_2\overline{Y}_2I(\pi ),`$ (55d)
where the $`U_n`$ correspond to the case where the needle is in position $`n`$.

The sequences of unitary transformations, e.g. $`X_1,Y_2`$,etc., are easily emulated on an ordinary computer. We use our software package, Quantum Computer Emulator (QCE) , to perform these calculations. An example will be shown in Appendix A. On an ideal QC, sequences (55) return the correct answer, i.e. the position of the searched-for item. This is easily verified on the QCE by selecting the elementary operations (called micro instructions on the QCE) that implement an ideal QC.

## IV Two-qubit QC’s

The energy-level structure of the nuclear spins of molecules such as deuterated cytosin and carbon-13 labeled chloroform can be described by model (1) and hence they can be used as physical realizations of 2-qubit QC’s.

### A Resonant pulses

NMR-QC experiments on carbon-13 labeled chloroform use resonant pulses to manipulate the quantum state of the nuclear spins of the 1H and 13C atoms. In the presence of a static magnetic field along the $`+z`$ direction this NMR-QC system is described by (1) with $`h_zg_1/2\pi 500\text{MHz}`$, $`h_zg_2/2\pi 125\text{MHz}`$, and $`J_z/2\pi 215\text{Hz}`$. A detailed account of simulations for this case have been published elsewhere. Simulations for physical model (1) confirm that sequences (55) yield the correct answers for the database search problem. However, we also demonstrated that the outcome of these calculations may be very sensitive to the order in which logically independent operations are carried out.. Some of these results are reproduced in Table I (first six rows, for details see ).

The results marked with a tilde are obtained by using a logically identical but physically different uniform superposition, i.e.

$`|u^{}=\overline{X}_1\overline{X}_1\overline{Y}_1\overline{X}_2\overline{X}_2\overline{Y}_2|00.`$ (56)
On an ideal QC, $`|u^{}=|u`$ but on a physical realizable machine this is unlikely to be the case. In an experiment it is simply impossible to freeze spin 2 (1) during the time that resonant pulses are being applied to spin 1 (2). Unless the length of these pulses is chosen judiciously, the wave function will acquire an additional phase. The corresponding unitary transformation does not necessarily commute with the operations that follow, potentially leading to an incorrect final result (as shown in Table I), as we pointed out in a previous paper.

For the two-spin system (56) one may optimize the pulse durations such that the effect of these phase errors yields qualitatively correct answers. A basic step is to make the pulse lengths commensurate with all relevant time scales. The main idea of the method suggested in and demonstrated in is to choose the parameters of the electromagnetic pulse such that the non-resonant spin still rotates but returns to its initial position at the end of the pulse . The results of the calculations are shown in Fig. 4 and summarized in Table I (primed symbols). In the numerical calculations we used $`J=10^6`$, $`h_zg_1=1`$, $`h_zg_2=0.25`$ , ($`h_yg_1=0.025`$, $`h_yg_2=0.0625`$, $`t_{\text{pulse}}=4\pi `$), for the pulses on the first spin and ($`h_yg_1=0.05`$, $`h_yg_2=0.0125`$, $`t_{\text{pulse}}=8\pi `$) for the pulses on the second spin, all values in dimensionless units. It is clear that optimization has the desired effect on the sequences that operate on $`|u^{}`$ (symbols with a tilde).

In spite of this optimization, Table I shows that there are significant quantitative differences between the theoretically exact results (rows (1,2)) and those obtained by simulating a physical model of a QC (e.g. rows 7 to 10). Even if the pulse length is taken to be commensurate with the relevant time scales of the QC, changing the state of the qubits by way of resonant pulses yields quantum states that are different from those obtained by means of the unitary transformations used in the analysis of the ideal QC. This is because spin 2 also interacts with the field applied to spin 1 and vice versa. With each program step, the non-ideal unitary operation may or may not result in the proliferation of errors. These errors are systematic (there is no “random” error source in our calculations) and directly linked to the structure of the QA. This is a clear case of a QPP, although we managed to let the different QA’s produce the correct answer. Note that the QPP cannot be solved by means of error correction: The operations on the extra qubits required for error correction will suffer from exactly the same QPP.

### B Stability of Grover’s quantum algorithm

So far we studied the stability of quantum algorithms by perturbing the input to the database encoding part of the algorithm. In this subsection we will study the QPP of the database query part of Grover’s search algorithm (the operation $`G`$, see (54)) itself. First we will assume that the input provided to $`G`$ is exact (i.e. of the form (46) for example) and we will compare the output of logically identical but physically different implementations of $`G`$. As examples we take the original sequence

$`G=X_1X_1\overline{Y}_1X_2X_2\overline{Y}_2Y_1\overline{X}_1\overline{Y}_1Y_2\overline{X}_2\overline{Y}_2I(\pi )X_1X_1\overline{Y}_1X_2X_2\overline{Y}_2,`$ (57)
and two, logically identical, sequences

$`\widehat{G}=X_1X_1\overline{Y}_1X_2X_2\overline{Y}_2Y_1\overline{X}_1Y_2\overline{X}_2\overline{Y}_1\overline{Y}_2I(\pi )X_1X_1\overline{Y}_1X_2X_2\overline{Y}_2,`$ (58)
$`\stackrel{~}{G}=X_1X_1\overline{Y}_1X_2X_2\overline{Y}_2Y_2\overline{X}_2\overline{Y}_2Y_1\overline{X}_1\overline{Y}_1I(\pi )X_2X_2\overline{Y}_2X_1X_1\overline{Y}_1.`$ (59)
Note that on purpose we did not “optimize” these sequences by using e.g. $`\overline{X}_1X_1=1`$. On an ideal QC we have

$`G|\mathrm{\Psi }_2=G^3|\mathrm{\Psi }_2=\widehat{G}|\mathrm{\Psi }_2=\widehat{G}^3|\mathrm{\Psi }_2=\stackrel{~}{G}|\mathrm{\Psi }_2=\stackrel{~}{G}^3|\mathrm{\Psi }_2=|01,`$ (60)
providing another test of the stability of the query operation $`G`$ on a physical QC. Table II contains the numerical results obtained by running the sequences (57), (58), and (59) on the QCE using the exact state $`|\mathrm{\Psi }_2`$ (see (54)) as input. In the case of the NMR-like QC optimized resonant pulses were used. The ideal QC performs as expected but the physical implementation (symbols with a prime) does not. In fact even one application of $`\stackrel{~}{G}`$ apparently returns an answer that is close to being useless ($`Q_1^{}0.5`$). As the three sequences (57), (58), and (59) are logically identical this is a clear case of a QPP.

The occurrence of a QPP seems to be a generic feature of QA’s running on QC’s. Therefore it is of interest to try to quantify the QPP. We now describe a simple procedure for this purpose, using $`G`$ and the case where the item is located in position 0 (the exact input state being $`|\mathrm{\Psi }_0`$) as an example. In general there are two sources of errors in this calculation: The input $`|\mathrm{\Psi }^{}`$ to $`G`$ and $`G`$ itself, the latter depending on the particular hardware implementation of the QC. As before we will use the resonant pulse technique in our numerical experiments.

We write $`|\mathrm{\Psi }^{}`$ as

$`|\mathrm{\Psi }^{}=\alpha _0|\mathrm{\Psi }_0+\alpha _1|\mathrm{\Psi }_1+\alpha _2|\mathrm{\Psi }_2+\alpha _3|\mathrm{\Psi }_3,`$ (61)
where the amplitude $`\alpha _0`$ can always be taken real ($`1\alpha _01`$) and is chosen at random. The other three complex coeffients are chosen randomly too, subject to the constraint $`|\alpha _1|^2+|\alpha _2|^2+|\alpha _3|^2=1|\alpha _0|^2`$. The real variable $`x\mathrm{\Psi }^{}|\mathrm{\Psi }_0=\alpha _0`$ is a measure for how much the input state deviates from the exact reference input $`|\mathrm{\Psi }_0`$. On an ideal QC, $`G|\mathrm{\Psi }_0=|00`$. Thus we can use the state $`|00`$ as reference to determine how much the output state $`|\mathrm{\Theta }G|\mathrm{\Psi }^{}`$ deviates from the exact answer. We quantify this deviation by the variable $`y=|00|G\mathrm{\Psi }^{}|`$.

The result of a numerical experiment using 20000 random input states $`|\mathrm{\Psi }_0`$ is shown in Fig. 1. Plots for the three other cases are nearly identical. We classify input-output pairs as “good” or “bad” as follows. First we choose a confidence level $`0c1`$ ($`c=0.7`$ for the data shown in Fig. 1). A particular input-output pair is considered to be good if $`x^2c`$ and $`y^2c`$. In Fig. 1 the good (bad) pairs are shown by black (gray) markers.

At a fairly low confidence level of $`c=0.7`$, the region of stable operation of the $`G`$ operation is rather small. This corroborates our earlier finding that successive applications of $`G`$, e.g. $`G^3`$, rapidly drive the system into a region of instability. In general, quantum systems are very sensitive to noise and become more sensitive as the number of operations on qubits is increased. In the absence of dissipation, it is easy for the system to leave the relatively small manifold of good input states, a characteristic feature of almost chaotic dynamics.

### C Hard nonselective pulses

Another physical implementation of a 2-qubit QC employs the nuclear spins of two 1H spin-1/2 nuclei in a solution of cytosine in D2O. This system can also be described by Hamiltonian (1). In the NMR experiments hard nonselective pulses are used to address the qubits. In this section we study the stability of QC operation for this physical realization of a QC.

As usual it is expedient to transform to another frame of reference that rotates with a constant frequency. This is accomplished by substituting in the TDSE

$`|\mathrm{\Phi }(t)=e^{ith_z(g_1+g_2)(S_1^z+S_2^z)/2}|\mathrm{\Psi }(t).`$ (62)
The time evolution of $`|\mathrm{\Psi }(t)`$ is then governed by the Hamiltonian

$`H=J_zS_1^zS_2^z{\displaystyle \frac{\mathrm{\Omega }}{2}}S_1^z+{\displaystyle \frac{\mathrm{\Omega }}{2}}S_2^z,`$ (63)
where $`\mathrm{\Omega }=h_z(g_1g_2)`$. Guided by experiment in our numerical work we will set $`\mathrm{\Omega }/2=1`$ and $`J_z/\pi \mathrm{\Omega }=0.01887`$ (in dimensionless units).

We now consider the time evolution of the two spins when we apply a static magnetic field $`h_x`$ along the $`x`$-axis. The Hamiltonian in the laboratory frame is given by

$`H=J_zS_1^zS_2^zh_z(g_1S_1^z+g_2S_2^z)h_x(g_1S_1^x+g_2S_2^x),`$ (64)
and the corresponding expression in the rotating frame of reference reads

$`H=`$ $`J_zS_1^zS_2^z{\displaystyle \frac{\mathrm{\Omega }}{2}}S_1^z+{\displaystyle \frac{\mathrm{\Omega }}{2}}S_2^z`$ (65)
$`h_x(g_1S_1^x+g_2S_2^x)\mathrm{cos}\mathrm{\Omega }th_x(g_1S_1^y+g_2S_2^y)\mathrm{sin}\mathrm{\Omega }t.`$ (66)
If the duration of the pulse is short, i.e. $`\mathrm{\Omega }t1`$, it is a good approximation to drop the time-dependence in (65) and we obtain

$`\stackrel{~}{H}={\displaystyle \frac{\mathrm{\Omega }}{2}}S_1^z+{\displaystyle \frac{\mathrm{\Omega }}{2}}S_2^zh_x(g_1S_1^x+g_2S_2^x),`$ (67)
where we used the fact that since $`|J_z||\mathrm{\Omega }/2|`$, for short pulses the effect of the spin-spin interaction on the time evolution is small and may be neglected. It is instructive to compute the time evolution of the spins, initially in state $`|0`$ ($`|0=|`$ by convention), under these circumstances. A straightforward calculation yields

$`Q_j=0|e^{it\stackrel{~}{H}}S_j^ze^{it\stackrel{~}{H}}|0={\displaystyle \frac{1}{2}}{\displaystyle \frac{[2h_xg_j(\mathrm{\Omega }/2+\lambda )]^2}{[(\mathrm{\Omega }/2+\lambda _j)^2+h_x^2g_j^2]^2}}\mathrm{sin}^2{\displaystyle \frac{\lambda _jt}{2}},`$ (68)
where $`\lambda _j=(\mathrm{\Omega }^2/4+h_x^2g_j^2)^{1/2}`$. Expression (68) shows that for hard pulses ($`|h_xg_j|\mathrm{\Omega }/2`$) and $`\lambda _jt=\pi `$, the effect of the pulse is to change qubit $`j`$ from $`|0`$ to approximately $`|1`$. Note however that a sequence of such pulses can never turn $`|0`$ into $`|1`$ exactly and that both spins are affected by the pulse.

In the rotating frame of reference the two spins rotate around the $`z`$-axis in the opposite sense. This fact can be exploited to perform operations that leave one spin intact while flipping the other spin . For instance, to rotate spin 1 around 90 degrees (clockwise) about the $`x`$-axis and without changing the state of spin 2, an operation we call $`X_1`$, we can use the sequence of pulses

$`X_1^{}=X_{1,2}^{}\overline{Y}_{1,2}^{}Z_{1,2}^{}Y_{1,2}^{},`$ (69)
where $`X_{1,2}^{}`$ denotes a 45-degree pulse around the $`x`$-axis acting simultaneously on both spins (see (67)), $`Y_{1,2}`$ a 90-degree pulse around the $`y`$-axis, and $`Z_{1,2}^{}`$ represents a pulse during which the spins evolve in time according to (67) and make a rotation of 45 degrees around the $`z`$-axis, (clockwise in the case of spin 1, anti-clockwise in the case of spin 2). There are many sequences of two-spin pulses that yield $`X_1`$ (or, as a matter of fact, any other rotation). For example the inverse of $`X_1^{}`$, $`\overline{X}_1^{}`$, can be obtained directly from (69), i.e.

$`\overline{X}_1^{}=\overline{Y}_{1,2}^{}\overline{Z}_{1,2}^{}Y_{1,2}^{}\overline{X}_{1,2}^{},`$ (70)
or can also be written as

$`\overline{X}_1^{}=\overline{X}_{1,2}^{}\overline{Y}_{1,2}^{}\overline{Z}_{1,2}^{}Y_{1,2}^{}.`$ (71)
For our numerical calculations we have chosen to work with the sequences

$`X_1^{}`$ $`=X_{1,2}^{}\overline{Y}_{1,2}^{}Z_{1,2}^{}Y_{1,2}^{},`$ (72a)
$`\overline{X}_1^{}`$ $`=\overline{X}_{1,2}^{}\overline{Y}_{1,2}^{}\overline{Z}_{1,2}^{}Y_{1,2}^{},`$ (72b)
$`X_2^{}`$ $`=X_{1,2}^{}\overline{Y}_{1,2}^{}\overline{Z}_{1,2}^{}Y_{1,2}^{},`$ (72c)
$`\overline{X}_2^{}`$ $`=\overline{X}_{1,2}^{}\overline{Y}_{1,2}^{}Z_{1,2}^{}Y_{1,2}^{},`$ (72d)
$`Y_1^{}`$ $`=Y_{1,2}^{}\overline{X}_{1,2}^{}Z_{1,2}^{}X_{1,2}^{},`$ (72e)
$`\overline{Y}_1^{}`$ $`=\overline{Y}_{1,2}^{}\overline{X}_{1,2}^{}\overline{Z}_{1,2}^{}X_{1,2}^{},`$ (72f)
$`Y_2^{}`$ $`=Y_{1,2}^{}\overline{X}_{1,2}^{}\overline{Z}_{1,2}^{}X_{1,2}^{},`$ (72g)
$`\overline{Y}_2^{}`$ $`=\overline{Y}_{1,2}^{}\overline{X}_{1,2}^{}Z_{1,2}^{}X_{1,2}^{},`$ (72h)
to perform the single-qubit operations that are used in the implementation of Grover’s search algorithm. As in the case of resonant pulses the natural time-evolution of the system provides the tool to carry out the conditional phase-shift operation $`I(\pi )`$.

The numerical values for the qubits 1 and 2 in the final state of the machine, obtained by running the QCE for this particular hardware realization of a QC and QA are are given in Table III.

From a comparison with the data of Table I we conclude that this hardware realization of a 2-qubit QC seems to be much less sensitive to QPP errors. As a matter of fact, in our numerical experiments (data not shown but included in the QCE download) we found none. The main reason for the good performance of this QC can be traced back to the fact that the QC contains two spins with exactly the same precession frequency, precessing in opposite direction. It is clear that it will become rather cumbersome if not impossible to design a $`N`$-qubit QC (if $`N2`$) using this approach: We would have to solve a complicated optimization problem that for each QA we would like to run on the QC.

### D Reducing systematic errors

For QC hardware to be of practical use, a basic requirement is that logically identical but physically different implementations of the same QA yield the same answers. The examples discussed above show that optimization of the pulses is essential to achieve this goal. For 2-qubit QC’s, the optimization problem is relatively easy to solve because only two frequencies are involved. Alternatively one may choose to use the same duration for all pulses and optimize the strengths of the pulses. However, in the case of $`N2`$ qubits, choosing pulse durations commensurate with $`N`$ different frequencies will considerably slow down the operation of the whole machine. In addition, the longer the pulses, the more important the effect of imperfections of the pulses becomes. Although there are ingenious techniques to optimize pulses in this respect, it may well be that in order to successfully run a QA on $`N`$-qubit QC hardware one first has to solve a rather complicated optimization problem and then simulate the QC by running the QA on a conventional machine.

## V Can dissipation reduce the QPP ?

Elsewhere we reported on the effect of dissipation on the quantum dynamics of nanoscale magnets. In these systems dissipation usually causes decoherence. Decoherence limits the time for performing quantum computations. Above we have demonstrated that programming non-ideal QC’s is difficult due to intrinsic instabilities of the physical device. Therefore we wonder if dissipation can provide the stabilizing processes, perhaps at the cost of reducing the time of coherent quantum operation.

Let us consider the procedure to generate the uniform superposition state from the state with all spins up. Above we made use of sequences such as $`XXY`$ but it is easy to see that from a programming-point-of-view $`Y`$ is just as good: In a classical picture both sequences change the direction of a spin from the $`z`$ to the $`x`$-axis. Thus, on an ideal QC we may omit $`XX`$ if we want. However, in an experiment (in which there always is some amount of dissipation), $`XX`$ can help to stabilize the direction of the spin once it points in the $`x`$-direction. Although in the pure quantum case $`XX`$ only causes the spin to precess around the $`x`$-axis, the damping of this motion by dissipation will let the spin relax to the $`x`$-direction.

We have studied the effects of dissipation on quantum computations using the master equation

$`{\displaystyle \frac{\rho (t)}{t}}=i[H,\rho (t)]\lambda \left([C,R\rho (t)]+[C,R\rho (t)]^{}\right),`$ (73)
that describes the time-evolution of the reduced (to the subspace of the qubits) density matrix $`\rho (t)`$. The operator $`R`$ is defined as

$`k|R|m=\zeta (E_kE_m)N_\beta (E_kE_m)k|C|m,`$ (74)
where $`H|k=E_k|k`$, $`H`$ being the Hamiltonian describing the QC. We take for the spectral density of the bosonic thermal bath $`\zeta (E)=I_0E^2\text{sign}(E)`$, i.e. an super-Ohmic reservoir, and $`N_\beta (E)=(e^{\beta E}1)^1`$ denotes the boson occupation number. The operator $`C`$ specifies the coupling between the spins and the bosons, e.g. $`C=_j(S_j^x+S_j^z)/2,`$. We set the inverse temperature $`\beta =100`$. The parameter $`\lambda `$ controls the flow of heat to and from the QC.

The procedure for simulating a QC in the presence of dissipation is very similar to the one used in the pure quantum case: The only difference is that we solve the master equation (73) instead of TDSE (2). For $`\lambda =0`$ the two completely different numerical methods used to solve these two equations (see  and  respectively) yield identical results.

We have studied the effect of dissipation by performing calculations for the same cases as described in Section IV A. Our simple-minded model of dissipation is found to cause only decoherence. We don’t find qualitative differences between the time-evolutions corresponding to the two different initial conditions. The presence of dissipation does not seem to affect the QPP. For $`\lambda =10^5`$ the value of the qubits in the final states is approximately 1/2, i.e. all useful information is lost (see Fig. 2). Apparently our choice does not lead to a relaxation during the $`XX`$ pulses. Although the dissipation processes incorporated through the use of master equation (73) lead to the correct thermal equilibrium state, there is no unique prescription to choose the bath Hamiltonian and the coupling between the bath and the spin system, i.e. the form of $`C`$. We are currently studying more realistic mechanisms for dissipation and will report on the results elsewhere.

## VI Four-qubit QC’s

In general practical applications of QC’s will require complicated QA’s. Therefore it is of interest to study how the instabilities due to imperfections in the operations affect the performance in such complicated QA’s. In this section we use a copy (or swap) procedure in a four-qubit QC as one of the simplest examples.

As discussed above, an implementation of Grover’s search algorithm on a QC consist of two separate parts: 1) The preparation of the uniform superposition state $`|u`$ and encoding of the database information and 2) the search (query) of the database. The algorithms described above use the same two qubits for these two different tasks. First they are used to store the information contained in the database and second they are used to carry out Grover’s algorithm to query the very same two qubits. Although this is sufficient to demonstrate the realization of 2-qubit QC hardware, from the point of view of computation this demonstration is of little use. Indeed, instead of organizing the computation in two steps (bringing the two qubits into a state that reflects the content of the database and then using the same two qubits to perform the query) we could have written the information from the database directly into the qubits in a form which reveals the position of the item in a trivial manner.

A computationally non-trivial demonstration of a QC running Grover’s database search algorithm requires the database and quantum processor to have their own qubits. This obviously means using at least four qubits, two to hold the database information (in superposition state) and two to process the query. Therefore we will need an operation to copy the state of the database into the quantum processor. A schematic picture of how the calculation is organized is shown in Fig. 3. Let us now study how instabilities in the physical operations affect the outcome of this more complicated computational procedure.

The copy operation transfers the amplitudes of an arbitrary linear combination of spin-up and spin-down of e.g. spin 1 to e.g. spin 2, initially in a state of spin up. More specifically

$`C_{1,2}(a|0_1+b|1_1)|0_2=|0_1(a|0_2+b|1_2).`$ (75)
In principle this can be done by a network of CNOT gates. One possible realization of a CNOT gate is the pulse sequence

$`C_{1,2}=Y_1X_1Y_1X_1\overline{Y}_1X_2\overline{Y}_2I(\pi )Y_2,`$ (76)
up to irrelevant phase factors, but there are many others. A fundamentally different sequence that performs the same task is

$`C_{1,2}=Y_2\overline{X}_2\overline{Y}_2I^{}(\pi ),`$ (77)
where

$`I^{}(2a)\left(\begin{array}{c}|00\hfill \\ |10\hfill \\ |01\hfill \\ |11\hfill \end{array}\right)\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \mathrm{cos}a& i\mathrm{sin}a& 0\\ 0& i\mathrm{sin}a& \mathrm{cos}a& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}|00\hfill \\ |10\hfill \\ |01\hfill \\ |11\hfill \end{array}\right).`$ (90)
A simple calculation shows that the two-qubit operation $`I^{}(2a)`$ corresponds to the time-evolution of a two-spin XY model:

$`I^{}(2\tau J_{xy})=e^{i\tau J_{xy}(S_1^xS_2^x+S_1^yS_2^y)}.`$ (91)
The other pulses in sequence (77) serve to remove unwanted phase factors. Note that (75) and (91) destroy the state of the database, a manifestation of the “no-cloning theorem”. As Fig. 3 shows this is not really a problem as the state of the database after the copy operation took place is the same state that was used to initialize the state of the database.

As before our aim is to address the fidelity of the results obtained by running a quantum algorithm on a hardware realization of a QC. Our QCE can simulate various candidate technologies and architectures (ways of interconnecting different units). However, to analyse the QPP it is at present sufficient to consider marginal extensions of two-qubit QC’s. With this in mind we have chosen to “implement” a 4-qubit QC as two identical Ising models

$`H=J_zS_1^zS_2^zh_z(g_1S_1^z+g_2S_2^z)J_zS_3^zS_4^zh_z(g_1S_3^z+g_2S_4^z).`$ (92)
In our numerical work we take for the values of the model parameters those corresponding to the chloroform molecule and we use selective resonant pulses to address the individual spins (see above). Thereby we assume, again for simplicity of analysis, that we can address pairs (1,2) and (3,4) separately (i.e. a pulse tuned to spin 1 does not affect spin 3 etc.). In the case where we use CNOT’s to perform the copy-qubit operation we add to (92) the Ising interactions

$`H_{1,3}`$ $`=J_zS_1^zS_3^z,`$ (93a)
$`H_{2,4}`$ $`=J_zS_2^zS_4^z,`$ (93b)
or, in the case the XY-model is used to copy the qubits from the database into the query processor, we add to (92)

$`H_{1,3}^{}`$ $`=J_{xy}(S_1^xS_3^x+S_1^yS_3^y),`$ (94a)
$`H_{2,4}^{}`$ $`=J_{xy}(S_2^xS_4^x+S_2^yS_4^y).`$ (94b)
It is obvious that these operations cannot be realized in terms of NMR technology. However this extension from a 2-qubit to a 4-qubit QC is sufficient to study in detail the stability of QA’s.

### A Resonant pulses

In our numerical simulations on the QCE we take the same model parameters as those of the chloroform molecule (see above) for both the database and query machine. We use resonant pulses to address the single qubits. Thereby we assume that qubits 1 and 2 can be shielded from qubits 3 and 4 during the application of these pulses.

The results of running Grover’s search algorithm on the QCE are summarized in Table IV. In these calculations we used the pulses optimized in the sense discussed above (i.e. when used to run on a 2-qubit QC they yield the correct answers, independent of the initialization sequence). Nevertheless we observe that the implementation that uses the CNOT gates performs rather poorly, at the border of giving the wrong answers. As in the 2-qubit case, the fact that the duration of the pulses is no longer commensurate with the Larmor periods of the spins is one source of errors. Another mechanism that causes errors to accumulate is discussed below.

### B Rotating field

Finally, in an attempt to further reduce the uncertainty on the outcome of the quantum computations, we have repeated the simulations using resonant pulses of rotating fields. More specifically, to apply a pulse to rotate for example spin 1, we add to the Hamiltonian a term of the form

$`H_1(t)=h_1S_j^y\mathrm{sin}f_1t+h_1S_j^x\mathrm{cos}f_1t,`$ (95)
where the frequency $`f_1`$ is tuned to the resonance frequency of spin 1 and the intensity of the pulse $`h_1`$ is determined by the desired angle of rotation. A simple calculation shows that such a pulse transforms the single qubit in exactly the same manner as the corresponding, ideal transformation. It is evident that this kind of optimization should improve the stability of the QA with respect to logically-allowed interchanges of operations. Further improvements can be made by reordering some of the elementary operations (recall that the sequence implementing e.g. a CNOT is not unique). The results of the simulations are presented in Table V. Compared to the results of Table IV there is one major change: The CNOT based implementation performs much better.

In this implementation there are two sources of errors: The effect of the pulses on non-resonant spins and the interaction-controlled phase shift. In the case of the chloroform molecule there is a large difference in time scale (roughly a factor $`10^6`$) between the precession frequencies and the spin-spin interaction. Therefore, during the execution of an interaction-controlled phase shift, the spins make a large number of rotations about the $`z`$-axis. We have seen that it really matters for the stability of the calculation whether the angles of rotation are a multiple of $`2\pi `$ or not. However, in view of the difference in time scale, this implies that the duration of an interaction-controlled phase shift must be specified with a sufficiently high accuracy in order to recover the correct final results.

## VII Summary

We have studied the difficulties that are encountered in programming quantum computer hardware. Taking Grover’s database quantum search algorithm as an example, we ran various, logically identical versions of the algorithm on different, physically realizable quantum computer hardware. We demonstrated that the choice of the physical processes used for quantum computation has direct consequences in terms of programming the machine.

We have shown that in the absence of dissipation non-ideal, physically realizable quantum computers operate in a regime of extreme sensitivity. This high sensitivity reflects itself in problems of programming quantum computers such that they perform correct calculations. These problems cannot be solved by error correction but may be reduced, and in some cases almost eliminated, at the cost of extensive, machine and program-specific optimization. A possible route to reduce this sensitivity may be to introduce some form of dissipation in a well-controlled manner. Dissipation can extend the region of stable operation of the quantum computer but also limits the time interval for quantum computer operation due to decoherence. We are currently exploring ideas along this line.

## Acknowledgement

Support from the Dutch NCF, FOM and NWO, and from the Grant-in-Aid for Research from the Japanese Ministry of Education, Science and Culture is gratefully acknowledged.

## A Quantum Computer Emulator (QCE)

In this appendix, we show an example of QCE. In Fig. 4 we show a snapshot of a window of QCE with a set of operations used to produce the results of Table I for the optimized pulses.

The left panel in Fig. 4 lists the fundamental operations, where MI stands for micro instruction. For example MI ‘-X1’ means the operation $`\overline{X}_1`$, etc. MI’s are defined by a table with interaction parameters, duration time, etc. This table appears when we double-click an MI. More details can be found in . The set of MI’s for the case of resonant optimized pulses is stored under the name ‘nmr-insta’. A QP represents a subroutine which consist of sequence of MI’s. For example QP ‘f2’ means the sequence of operations $`F_2`$, see Eq. (47). The MI ‘initialize’ simply prepares the state $`|00`$, MI ‘inv-means’ performs the inversion about the mean (the essential step in Grover’s algorithm, see Eq. (53)), MI ‘tau’ corresponds to the operation $`I(\pi )`$ (see Eq. (44)), QP ‘prepare’ denotes the preparation of $`|u`$ while QP ‘prapera’ denotes the preparation of $`|u^{}`$, etc. The MI’s ‘tau2’, ‘tau\_insta’, and ‘w1’ are not used in this example.

The eight small windows show the programs corresponding to $`U_0,\mathrm{}U_3`$ using $`|u`$ (upper level) and those using $`|u^{}`$ (lower level), which are named QP ‘g0’ to QP ‘g3prap1̃’. The result of running each program is shown in the grid below the QP. The expectation value of $`(1S_i^\alpha )/2`$, $`i=1`$ or 2, and $`\alpha =x,y`$ or z, can be read off from a particular cell by moving the mouse over the cell (not shown). In Fig. 4 the approximate value can be determined from the darkness of the cell: Light gray corresponds to zero, dark gray to one. For example, the value of the cell $`(z,1)`$ in QP ‘g0’ is 0.027 and that of $`(z,2)`$ is 0.152, as given in Table I ($`Q_1^{}`$ and $`Q_2^{}`$ respectively). For other cases, e.g. the non-optimized pulses, we use another set of MI’s (‘nmr’) with parameters given in .


text:
# Opportunities for future supernova studies of cosmic acceleration

\[

## Abstract

We investigate the potential of a future supernova dataset, as might be obtained by the proposed SNAP satellite, to discriminate among different “dark energy” theories that describe an accelerating Universe. We find that

many such models can be distinguished with a fit to the effective pressure-to-density ratio, $`w`$, of this energy. More models can be distinguished when the effective slope, $`dw/dz`$, of a changing $`w`$ is also fit, but only if our

knowledge of the current mass density, $`\mathrm{\Omega }_m`$, is improved. We investigate the use of “fitting functions” to interpret luminosity distance data from supernova searches, and argue in favor of a particular preferred method, which we use in our analysis.

\]

Increasing evidence that the Universe is accelerating  confronts us with deep unresolved questions. As yet no compelling understanding of the acceleration has been achieved but many models have been proposed, typically introducing a form of “dark energy” to account for cosmic acceleration. Further resolving the observational evidence for acceleration is certain to have a great impact on our understanding of the “dark energy” and of the underlying fundamental questions.

Type Ia supernovae can be used as standard candles to infer the luminosity distance ($`d_L`$) as a function of redshift , and such data provide a key element in the case for cosmic acceleration. By analyzing a simulated dataset, as might be obtained by the proposed SNAP satellite, we can test the ability of such experiments to distinguish among currently attractive models. Our work will be reported more extensively in . We note that where our calculations overlap we are in complete numerical agreement with the recent paper by Maor et al. . However our perspective and emphasis is considerably different.

We consider a future supernova dataset that has been converted to a table of effective magnitudes and redshifts of objects with a single fiducial absolute magnitude $`M`$. We will consider both statistical and systematic uncertainties in the magnitudes. Typically the redshifts are known to sufficiently high precision that their uncertainties can be ignored.

The magnitudes and redshifts are related by the luminosity distance according to

$$m(z)=+5\mathrm{log}𝒟_L,$$
(1)
where

$$𝒟_L=H_0d_L=H_0(1+z)_0^z\frac{c}{H\left(z^{}\right)}𝑑z^{}$$
(2)
and $`M5\mathrm{log}H_0+25`$. When nearby supernovae are used to determine $``$ the uncertainties in the Hubble constant $`H_0`$ drop out of Eqn. 1. It is through $`𝒟_L(z)`$ that specific models of the evolution of the universe enter the picture. If we had just a few well-specified models of dark energy to compare, then their predictions for $`𝒟_L`$ could be fit to the supernova dataset and their relative likelihoods determined. Even if a model had a free parameter, that parameter could be determined to some accuracy. However, it is not understood what fundamental physical principles specify the actual dark energy in the Universe, so we are left with a potentially infinite set of families of theoretical models to compare, each with a potentially large number of free parameters. No finite dataset will ever select among all of these theoretical models (see, for example, Maor et al.), who show nine models, from this infinite set, that cannot be distinguished by a proposed supernova dataset).

In the face of this theoretical uncertainty it is useful to select a fitting function whose parameters can characterize the $`m(z)`$ dataset independent of any specific theory. By choosing a function that also fits most extant theories’ $`𝒟_L`$ predictions well, with the fewest possible parameters, we can then use these parameters as indices that label the theories so they can be compared with the dataset. Although it is also desirable for the fitting parameters to represent physically meaningful concepts of the current meta-theories, this is not necessary. This approach has already been undertaken by several authors in the context of supernova datasets. Each of these papers uses different types of functions to fit $`𝒟_L(z)`$. Here we describe a reason to choose one of these functions.

Original reports of supernova data used (two-parameter) models based on adding a cosmological constant term to Einstein’s equations. This approach was generalized in  to an effective, but constant, equation of state factor. Huterer and Turner  modeled the luminosity distance according to the simple power law $`d_L(z)=_{i=1}^Nc_iz^i`$. Saini et al.  used the form111 As we completed this work a preprint came out using a fitting function with even more parameters. We will discuss this function in

$$𝒟_L(z)=2(1+z)\left[\frac{z\alpha \sqrt{1+z}+\alpha }{\beta z+\gamma \sqrt{1+z}+2\alpha \gamma }\right].$$
(3)
Our preferred fitting function is motivated by cosmic acceleration driven by an extra energy component (“quintessence” or “dark energy”) with the pressure $`p`$ and density $`\rho `$ related by an equation of state $`p=w\rho `$. The cosmological constant case is reproduced when $`w=1`$. For more general constant values of $`w`$ this equation of state implies $`\rho (1+z)^{3(w+1)}`$ To allow for a $`z`$-dependence of $`w`$, it has been expanded in a power series in $`z`$ (as in ) or,equivalently, in $`(1+z)`$:

$$w(z)=\underset{i=0}{\overset{N}{}}w_iz^i=\underset{i=0}{\overset{N}{}}\stackrel{~}{w}_i\left(1+z\right)^i.$$
(4)
We assume a flat universe (theoretically preferred by inflation and currently favored by CMB data). Thus the remaining parameter required to fix the cosmology is $`\mathrm{\Omega }_m`$, the density of ordinary (pressureless) matter today in units of the critical density. For these models we get

$`𝒟_L(z)`$ $`=`$ $`c(1+z){\displaystyle \underset{0}{\overset{z}{}}}(1+z^{})^{3/2}[\mathrm{\Omega }_\mathrm{m}+\mathrm{\Omega }_Q(1+z^{})^{3\stackrel{~}{w}_0}`$ (6)
$`\times \mathrm{exp}\left\{3{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{\stackrel{~}{w}_i}{i}}[(1+z^{})^i1]\right\}]^{1/2}dz^{}`$
where for a flat universe $`\mathrm{\Omega }_Q=1\mathrm{\Omega }_\mathrm{m}`$. Equation 5 takes on a simple form in terms of $`\stackrel{~}{w_i}`$, but throughout the rest of this letter we use the more intuitive quantity $`w_i`$.

In  we survey essentially all the published quintessence models, and in every case the $`w`$-expansion gave a better fit than either Eqn. 3 or the power law expansion of $`d_L`$, when taken to the same order of fit variables (see Fig. 1). The quadratic order $`d_L`$-expansion (top dashed line) results in a worse fit than the linear order $`w`$-expansion (top solid line), and the quadratic order $`w`$-expansion (lower solid line) fits as well as the cubic order $`d_L`$-expansion (lower dashed line). (This result also holds when we weight the data points according to the dataset specification described below, although all the fits improve.)

The $`w`$-expansion parameters also have an intuitive meaning in their own right. Given the uncertainty regarding causes of cosmic acceleration, it might be reasonable to simply think of the $`w_i`$’s as the fundamental parameters for now. They tell us something relatively straightforward about the nature of the dark energy. The same cannot be said of the $`c_i`$ of  or the $`\alpha ,\beta `$ and $`\gamma `$ of . We are fortunate that the twin criteria of intuitive interpretation of parameters and good fits to published models point to the same choice of fitting functions. It is worth noting however that we have not proven that the $`w`$-expansion provides the best  possible fits to current models. It is possible that an even better choice could be discovered in the future.

Maor et al.  also favor the same $`w`$-expansion formalism. However we do not agree completely with their discussion of the alternative methods . Despite claims to the contrary in  all the approaches use “fitting functions”. If for example a new understanding of dark energy emerged in which the $`c_i`$ of  were fundamental parameters, then the $`w`$-expansion would provide a bad fit to those models and would be a bad choice for the fitting functions (whereas $`d_L(z)=_{i=1}^Nc_iz^i`$ would be an “exact expression”).

So we have argued that the $`w`$-expansion offers a good means of interpreting supernova data. Now we apply this expansion to simulated future supernova datasets.

For concreteness, we assume a future dataset similar to one proposed for the SNAP satellite project, with SNe observed uniformly within four different redshift ranges with the following different sampling rates: In the first range from $`z=00.2`$ we assume that their are $`50`$ observations, in the second and largest redshift range from $`z=0.21.2`$ there are $`1800`$ SNe and in the two high redshift bins, $`z=1.21.4`$ and $`z=1.41.7`$, there are $`50`$ SNe and $`15`$ SNe observations respectively. The statistical error in magnitude is assumed to be $`\sigma _{\mathrm{mag}}=0.15`$, including both measurement error and any residual intrinsic dispersion after calibration. It is worth noting that considerations such as those raised in  could change this redshift-distribution strategy to further optimize the impact of the data. We construct Monte Carlo simulated datasets, first for a flat, cosmological constant model with $`\mathrm{\Omega }_\mathrm{m}=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. Figure 2 shows an example of simulated data. The points show binned data with statistical error bars, where we combine $`56`$ SNe in each bin.

Figure 2 shows the magnitude difference, $`\mathrm{\Delta }m`$, between several dark energy models (solid) and the Monte Carlo fiducial model, along with the simulated binned data for that model. To guide the eye we also show a grid of constant $`w`$ models (dotted). Note how for the supergravity model  (labeled “SUGRA”), adding the linear term, $`w_1`$, improves the fit considerably (as shown by the dashed line, mostly covered by the solid line), compared to the best constant fit with $`w_0=0.74`$. For the two exponential model  (“2EXP”) the constant fit seems already sufficient. The line labeled “TOY” corresponds to a toy model with $`w_0=0.6`$ and $`w_1=0.8`$, which is not well represented by a constant fit. All the models have $`\mathrm{\Omega }_m=0.3`$.

Figure 2 illustrates that these models can be easily

differentiated from one another by the simulated data, and that for the toy model the contribution

from $`w_1`$ looks large enough to be discriminated by the data. The targeted SNAP systematic uncertainties can be estimated by considering for example a linear drift up or down in the data (as a function of $`z`$) which spans one error bar height over $`0<z<1.5`$ (Note that a constant systematic uncertainty would not affect the results.) Further investigation is required to determine if these targets are realistic. When the targeted systematic uncertainties are taken into account the “2EXP” and $`w=1`$ models can no longer be discriminated. In  we report more thoroughly on the models which can and cannot have $`w_1`$ differentiated from zero by such datasets. In most published models $`w_1`$ is below the level observable by supernovae alone.

We shall see below that the effect of introducing $`w_1`$ in the fit for this dataset is to substantially increase the uncertainty in $`w_0`$ (as emphasized in ). This is because with the additional fitting parameter there are many more ways of fitting the simulated data to good accuracy. This fact has nothing to do with the ability of the (in this case simulated) dataset to discriminate among the models shown. The ability to discriminate among the models is a fundamental property of the dataset, whereas the choice of fitting function is a matter of choosing suitable analysis techniques. As with any experimental data, for a given dataset each of the shown models is a member of a degenerate family of models that cannot be differentiated by the data. Opening up the space of fitting functions allows one to probe this degeneracy. Furthermore, we show below how realistic improvements in our independent knowledge of $`\mathrm{\Omega }_m`$ can break this degeneracy to an interesting degree.

Figure 3 shows joint likelihood contours in the $`\mathrm{\Omega }_mw_0`$ plane for the simulated data, when these are the only parameters used to fit the data. We show also the shift due to a linearly drifting systematic error of $`\pm 0.02\mathrm{mag}`$ per 1.5 units in redshift. Clearly the data will discriminate among the three models shown.

We now analyze the same simulated data using a fitting function that includes $`w_1`$. For this analysis (results shown in Table 1 and Figure 4), we consider different possible states of prior information about the value of $`\mathrm{\Omega }_m`$, ranging from relatively poorly known ($`\sigma _{\mathrm{\Omega }_m}0.15`$), roughly the case today, to well determined ($`\sigma _{\mathrm{\Omega }_m}0.04`$), a possible goal within the next decade. We show in  that the systematic errors lead to roughly the same contribution as the statistical errors for the case of a perfectly constrained $`\sigma _{\mathrm{\Omega }_m}=0.0`$. We conclude from Fig. 4 and Table 1 that expanding the fitting function to include $`w_1`$ will become useful when the value of $`\mathrm{\Omega }_m`$ is better constrained than it is today. This result adds to the case for producing an independent determination of $`\mathrm{\Omega }_m`$.

It may be possible to reduce the statistical uncertainties even further, by increasing the sample size and/or by improving the intrinsic luminosity calibration for each supernova. Figure 4 also shows the smaller likelihood region (dotted ellipse) which one gets when taking $`\sigma _{\mathrm{mag}}=0.09`$ and doubling the number of supernovae in each bin. This means (see Table 1) that the $`1\sigma `$ statistical error on $`w_1`$ improves from $`\sigma _{w_1}=0.6`$ for the current constraints on $`\mathrm{\Omega }_m`$ (solid line in Fig.4), to $`\sigma _{w_1}=0.2`$ for future $`\mathrm{\Omega }_m`$ constraints (dashed line), to very tight constraints of $`\sigma _{w_1}=0.12`$ for the “improved” scenario, with twice as many more-tightly-calibrated supernovae and $`\sigma _{\mathrm{mag}}=0.09`$ (dotted line).

Figure 4 is similar to Fig. 2 from . (To compare actual numbers for predicted uncertainties, note that Maor et al. quote the full-width range of 95%-confidence contours, i.e. four times a 1$`\sigma `$ uncertainty. Also, our different contour shapes are due to our choice of a Gaussian probability distribution for $`\mathrm{\Omega }_m`$ vs. their tophat distribution.) Maor et al. emphasize the uncertainties that emerge when the fitting function is expanded to include $`w_0`$ and $`w_1`$ given our current knowledge of $`\mathrm{\Omega }_m`$. We disagree with this emphasis and focus instead on the great impact a SNAP-type experiment will have. Even if our determination of $`\mathrm{\Omega }_m`$ were not to improve over the next decade, SNAP-class datasets would still be powerful discriminators among models. For example, such a dataset could easily discriminate between $`w=1`$ (a cosmological constant) as the source of cosmic acceleration and several interesting quintessence models (such as SUGRA or others discussed in ). A result either way would have profound

implications. Furthermore, significant improvements on the determination of $`\mathrm{\Omega }_m`$ are expected over the next decade and we have shown how sufficient improvements would allow the extra fitting parameter $`w_1`$ to become an “observable” and help to further differentiate theories .

We have investigated different analysis methods for supernova datasets, and used our preferred method to explore the prospects of further supernova searches. It is clear that a dataset of the sort proposed with SNAP  presents an exciting opportunity to constrain theories of cosmic acceleration.

ACKNOWLEDGMENTS: We thank S. Perlmutter for extensive discussions and in particular for emphasizing the importance of the improving knowledge of $`\mathrm{\Omega }_m`$; D. Huterer, A. Lewin, P. Steinhardt and M. Turner for helpful conversations; and

P. Astier for comparing our results with his independent calculations. We were supported by DOE grant DE-FG03-91ER40674 and U.C. Davis.


text:
# Magnetic field effects on quantum ring excitons

## I  Introduction

The fabrication of nanometer-size semiconductor structures by different techniques (including lithography, etching, direct chemical synthesis, and self-assembly, to name a few) has allowed a veritable explosion of activity in this area.  It is now well-known that carrier confinement into dimensions of a few tens of nanometers provides strong blue shift of the photoluminescence features from that in the original bulk material, a clear consequence of quantum confinement in these quantum dots. Currently, researchers are investigating a variety of magnetocapacitance and optical properties of dots,  including the role of inelastic light scattering and phonon confinement,  as well as Pauli blocking and other few- and many-particle effects in these systems.

In recent work, however, a new geometry of semiconductor quantum rings has been introduced in experiments of magnetocapacitance and infrared excitation for few-electrons.  Although metallic rings have been the subject of considerable attention for a number of years,  this geometry had not been achieved in semiconductors for sizes such that the electrons would propagate coherently (non-diffusively) throughout the ring. However, the self-assembled quantum rings now achieved are so small (with characteristic inner/outer radius of 20/100 nm and 2-3 nm in height), that they allow the study of a non–simply-connected geometry where carriers are coherent all throughout. It is clear that not only the single-particle states are interesting in this geometry (specially their behavior under magnetic flux), but also the role of interactions between particles (be it electrons or holes). Lorke et al. have shown that multiple-electron states in this geometry experience phase coherent effects in the presence of magnetic fields,  much as predicted by theory.  The question of the observability of similar coherent effects for electrons and holes around the loop, in the presence of a magnetic field, is fascinating, and some theoretical works have already begun to explore this regime.  Although beautiful experiments of optical emission in charge-tunable quantum rings have been recently presented,  they study the role of multiply-charged exciton complexes with no applied magnetic field.

The multiply-connected geometry of the semiconductor quantum rings adds an interesting dimension to the strong Coulomb effects in magnetic fields which have been explored in quantum confined systems. Excitons in magnetic fields have been investigated in structures such as quantum wells,  quantum wires,  and quantum dots,  as example of multi-carrier systems. For the ring geometry, one question that arises naturally is whether there is sensitivity of the exciton properties to an applied flux. This ‘Aharanov-Bohm effect’ (ABE) for an exciton is an interesting concept, since one would naturally associate the existence of the ABE with a net charge (as the coupling constant to the vector potential), and the net charge of this object is clearly zero. However, one could argue that the composite nature of the excitons would perhaps allow for a non-vanishing effect. In fact, for the case of particles constrained to move along a one-dimensional ring, rigorous derivations predict a non-zero ABE for the exciton, which will show in its various energy states and the associated dipole oscillator strength, for small enough rings.

In this work, we present calculations of the excitation spectrum and oscillator strength of excitons in rings pierced by magnetic fields. We explore the role of different confinement potentials and calculate binding energies, exciton sizes and their dependence on magnetic fields, as well as oscillator strengths which would be measurable in photoluminescence experiments, for example. Similar to the case of quantum dots, we find strong orbital effects from the magnetic field, which provides for effectively stronger confinement and accompanying diamagnetic level shifts, as well as splitting of some levels. These changes are found to be monotonic with field. Part of our motivation in this study is to explore the question of how robust is the ABE predicted in 1D rings, when one considers finite width and confinement potentials. The models we use are designed to mimic the situation in real semiconductor quantum rings achieved to date. Much to our chagrin, we find no trace of the predicted ABE for realistic values of the rings and fields in the problem. Although this negative result would suggest that it is difficult that experiments would measure this effect, it is still open to see to what extent highly sensitive experiments might be able to yield a positive result.

The remainder of this paper is organized as follows. In Sec. II, we present the model for the system and solution method. As a first approximation, the quantum ring structure is modeled by a parabolic confinement potential, in which the wavefunctions expressed in terms of the center of mass and relative coordinates are used as a basis set. This confinement potential has been experimentally confirmed by the recent experiments of Lorke et al. . In order to explore the role of different potentials, we also model the ring system with a hard-wall confinement, using wavefunctions expressed in terms of Bessel functions for electron and hole as a basis set. In Sec. III, we discuss the main effects of the magnetic field effects on the exciton characteristics, including the binding energy, electron-hole separation, and the linear optical susceptibility. Finally, we summarize our results in Sec. IV. The Appendix contains an outline of the derivation of the Coulomb matrix element with these basis functions. The analytical expressions presented there greatly simplify our calculations.

## II Theoretical Model

Our model is a two-dimensional exciton in a quantum ring and in a static magnetic field, simulating recent experimental quantum ring structures. The presence of magnetic fields oriented along the $`z`$ axis, perpendicular to the plane of the ring, induces the electron and hole to perform classical orbits along the circumference, which of course yield quantization of the angular momentum in that direction.  The ring structures are well approximated by using parabolic potentials, giving soft confinement barriers, appropriate to samples produced by self-assembly.  For narrow rings (with steep confinement), however, the parabolic confinement and associated wavefunctions fail in a real system, as the increased confinement may push the levels into the anharmonic part of the potential, and even produce deconfinement of one carrier (typically the electron), as found in some calculations in quantum dots.  We also consider the case of hard-wall confinement and analyze the different results.

As the quantum rings and excitonic states under consideration are much larger than the unit cell of the material, the effective-mass approximation is a suitable approach, and is given by $`H=H_e+H_h+H_{eh}`$, where the subscripts $`e`$ and $`h`$ represent electron and hole, and the last term is the electron-hole interaction. The expression for the Hamiltonian of the electron in magnetic fields (in a parabolic-band approximation) is given by

$$H_e=\frac{1}{2m_e}\left[𝐩_e+\frac{q_e}{c}𝐀_e\right]^2+V_e,$$
(1)
with a similar expression for the Hamiltonian of the hole. Here, $`V_e`$ (or $`V_h`$) is the ring confinement potential for electron (hole), and naturally $`q_e=|e|`$ and $`q_h=+|e|`$. For parabolic confinement potential across the width of the ring, we use

$$V_i=\frac{1}{2}m_i\omega _i^2(r_ir_o)^2,$$
(2)
where the mean radius of the ring is $`r_o`$, and the characteristic confinement energy is $`\mathrm{}\omega _i`$, giving a characteristic ring width $`2\sqrt{\mathrm{}/m_i\omega _i}`$ for each particle. Here, $`i=e,h`$ represents the different particles. We choose the fully symmetric gauge vector potentials $`𝐀_e=\frac{1}{2}𝐁\times (𝐫_e𝐫_h)`$ and $`𝐀_h=\frac{1}{2}𝐁\times (𝐫_h𝐫_e)`$, for electron and hole, respectively, following Ref. \[\]. The Coulomb interaction term between carriers is given by $`H_{eh}=e^2/ϵr_{eh}`$, screened by the background dielectric constant $`ϵ`$.

### A Parabolic confinement potential

For the parabolic confinement potential is convenient to separate the problem into center of mass and relative coordinates, described as usual by $`𝐫=𝐫_e𝐫_h`$, and $`𝐑=(m{}_{e}{}^{}𝐫_{e}^{}+m_h𝐫_h)/M`$, where the total and reduced masses are $`M=m{}_{e}{}^{}+m_h`$, and $`\mu =m{}_{e}{}^{}m{}_{h}{}^{}/M`$. The total Hamiltonian can then be re-expressed as $`H=H_{CM}+H_{rel}+H_{mix}`$, with individual terms

$`H_{CM}`$ $`=`$ $`{\displaystyle \frac{1}{2M}}P^2+{\displaystyle \frac{1}{2}}M\omega _o^2(Rr_o)^2,`$ (3)
$`H_{rel}^o`$ $`=`$ $`{\displaystyle \frac{1}{2\mu }}p^2+{\displaystyle \frac{1}{2}}\mu \omega ^2r^2,`$ (4)
and $`H_{rel}=H_{rel}^o+H_{rel}^{}`$, where

$`H_{rel}^{}`$ $`=`$ $`\omega _c\gamma L_z{\displaystyle \frac{e^2}{ϵr}},`$ (5)
and (6)
$`H_{mix}`$ $`=`$ $`{\displaystyle \frac{e}{Mc}}(𝐁\times 𝐫)𝐏`$ (9)
$`m_e\omega _o^2r_oR\left(1+{\displaystyle \frac{2m_h}{MR^2}}𝐑𝐫+{\displaystyle \frac{m_h^2}{M^2}}{\displaystyle \frac{r^2}{R^2}}\right)^{1/2}`$
$`m_h\omega _o^2r_oR\left(1{\displaystyle \frac{2m_e}{MR^2}}𝐑𝐫+{\displaystyle \frac{m_e^2}{M^2}}{\displaystyle \frac{r^2}{R^2}}\right)^{1/2},`$
where $`\gamma =(m_hm_e)/M`$ depends on the mass asymmetry of the carriers, and we have set $`\omega _e=\omega _h=\omega _o`$. We have also denoted the relative angular momentum in the $`z`$-direction as $`L_z=(𝐫\times 𝐩)_z`$, and the effective confinement frequency as $`\omega ^2=\omega _o^2+\omega _c^2`$, with $`\omega _c=eB/(2\mu c)`$, resulting from the combined confinement of the potential and the magnetic field.

The main purpose in the change of variables above is to use the solutions of $`H_{CM}`$ and $`H_{rel}^o`$ as a basis for the solution of the full Hamiltonian. The center of mass basis is essentially a harmonic oscillator, with wavefunction $`\psi _{N,l_{CM}}`$ centered about $`r_o`$,

$`\psi _{N,l_{CM}}(R)`$ $`=`$ $`\alpha \sqrt{{\displaystyle \frac{2N!}{(N+|l_{CM}|)!}}}{\displaystyle \frac{1}{\sqrt{2\pi }}}e^{il_{CM}\theta }e^{\alpha ^2(Rr_o)^2/2}`$ (10)
$`\times `$ $`\left[\alpha (Rr_o)\right]^{|l_{CM}|}L_N^{|l_{CM}|}\left(\alpha ^2(Rr_o)^2\right),`$ (11)
with $`\alpha =\sqrt{M\omega _o/\mathrm{}}`$, and eigenvalues $`E_{\mathrm{CM}}^o=\left(2N+1+|l_{CM}|\right)\mathrm{}\omega _o`$. Here, $`N`$ and $`l_{CM}`$ are quantum numbers of the radial and angular momentum part, respectively, for the center of mass coordinates, and $`L_N^{|l_{CM}|}`$ is the associated Laguerre polynomial.  Similarly, the eigenvalues and eigenfunctions for the non-interacting relative Hamiltonian are given by wavefunction $`\varphi _{n,l}`$ and energy $`E_{rel}^o`$, where

$`\varphi _{n,l}(r)`$ $`=`$ $`\beta \sqrt{{\displaystyle \frac{2n!}{(n+|l|)!}}}{\displaystyle \frac{1}{\sqrt{2\pi }}}e^{il\phi }e^{\beta ^2r^2/2}\left(\beta r\right)^{|l|}L_n^{|l|}\left(\beta ^2r^2\right),`$ (12)
$`E_{\mathrm{rel}}^o`$ $`=`$ $`\left(2n+1+|l|\right)\mathrm{}\omega ,`$ (13)
with $`\beta =\sqrt{\mu \omega /\mathrm{}}`$. Here, $`n`$ and $`l`$ are quantum numbers of the radial and angular momentum part, respectively, for the relative coordinates. With this harmonic basis set, matrix elements for the Coulomb interaction, magnetic field dependent and mixing terms can be calculated analytically.  These expressions are collected in the Appendix.

The total Hamiltonian given above is then diagonalized numerically, leading to the eigenvalues and eigenfunctions. All the physical properties of the exciton in the ring can in principle be extracted from these eigenvalues and eigenfunctions. Here, we present the binding energy, electron-hole separation, and the linear optical susceptibility of the quantum ring. These quantities are readily accessible in optical experiments of photoluminescence (PL) and photoluminescence excitation (PLE). Denoting the wavefunctions of the exciton as $`|\mathrm{\Psi }=_{Nl_{CM}nl}a_{Nl_{CM}nl}|N,l_{CM},n,l,`$ with coefficients $`a_{Nl_{CM}nl}`$ obtained from the diagonalization, one can write for example an expression for the electron-hole separation,

$`r_s^2`$ $`=`$ $`\mathrm{\Psi }|r^2|\mathrm{\Psi }=\delta _{N^{},N}\delta _{l_{CM}^{},l_{CM}}\delta _{l^{},l}{\displaystyle \underset{Nl_{CM}nlN^{}l_{CM}^{}n^{}l^{}}{}}a_{N^{}l_{CM}^{}n^{}l^{}}^{}a_{Nl_{CM}nl}`$ (14)
$`\times `$ $`\sqrt{{\displaystyle \frac{n!n^{}!}{(n+|l|)!(n^{}+|l^{}|)!}}}{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \underset{j=0}{\overset{n^{}}{}}}{\displaystyle \frac{(n+|l|)!(n^{}+|l^{}|)!}{(k+|l|)!(j+|l^{}|)!}}`$ (15)
$`\times `$ $`(1)^{k+j}{\displaystyle \frac{1}{k!(nk)!j!(n^{}j)!}}{\displaystyle \frac{1}{\beta ^2}}\left(|l|+k+j+1\right)!.`$ (16)
Similarly, the linear optical susceptibility is given by $`\chi (\omega )=_j\left|0|P|1_j\right|^2\left(\mathrm{}\omega E_ji\mathrm{}\mathrm{\Gamma }\right)^1`$, where $`0|P|1_j`$ is the dipole matrix element between one electron-hole pair $`j`$ state and the vacuum state. These are proportional to the bulk interband matrix elements, $`p_{cv}`$, and can be written in terms of the envelope function as,

$`\left|0|P|1\right|^2`$ $`=`$ $`\left|p_{cv}\right|^2\left|{\displaystyle \mathrm{\Psi }(𝐫_e,𝐫_h=𝐫_e)𝑑𝐫_e}\right|^2`$ (17)
$`=`$ $`\left|p_{cv}\right|^2\left\{\delta _{l,0}{\displaystyle \underset{n}{}}a_n\sqrt{{\displaystyle \frac{\mu \omega }{\pi \mathrm{}}}}\right\}^2\left\{\delta _{l_{CM},0}\sqrt{2\pi }\sqrt{{\displaystyle \frac{2N!}{(N+|l_{CM}|)!}}}(1)^N{\displaystyle \frac{1}{\alpha }}\right\}^2.`$ (18)
### B A hard-wall confinement potential

Given that the two-dimensional free exciton size (effective Bohr radius) in InAs is $`a_B^{2D}16`$ nm (6 nm for GaAs), the quantum rings with widths larger than $`2a_B^{2D}`$ would tend to yield highly symmetric (nearly circular) ground states of the exciton, with the confinement potential being a small perturbation. For narrower quantum rings, however, the symmetry of the exciton in the ring would be strongly affected, and become increasingly one-dimensional. This would be favorable for the appearance of the ABE, as predicted by theory.  To allow for this different case, and so as to test the possible bias of the numerical calculations in the parabolic potential, we have also implemented solutions of the problem in a hard-wall confinement potential basis.

In that situation, the basis set for the exciton problem is given by products of the radial and angular parts of electron and hole, $`\mathrm{\Psi }=\mathrm{\Psi }_e(r_e,\varphi _e)\mathrm{\Psi }_h(r_h,\varphi _h)`$, where the individual wavefunctions are given by (in the absence of magnetic field)

$$\mathrm{\Psi }(r_i,\varphi _i)=\psi _i(r_i)\frac{1}{\sqrt{2\pi }}e^{il_i\varphi _i},$$
(20)
where $`i=e,h`$, and $`l_i`$ is an integer. The wavefunctions of the radial part must satisfy the hard-wall boundary conditions and vanish at both the inner ($`a`$) and outer radius ($`a+2L`$) of the ring structure. As such, they are given by linear combinations of Bessel functions, $`\psi _i(r_i)=AJ_{l_i}(k_ir_i)+BN_{l_i}(k_ir_i)`$, for $`ar_ia+2L`$. Here, $`A`$ and $`B`$ are normalized constants, and $`J_{l_i}`$ and $`N_{l_i}`$ are the $`l_i`$th-order Bessel function of the first and second kind, respectively, with $`k_i=\sqrt{2m_iE_i/\mathrm{}^2}`$. The eigenvalue conditions are obtained from the secular equation $`J_{l_e}(k_ea)N_{l_e}(k_e(a+2L))=N_{l_e}(k_ea)J_{l_e}(k_e(a+2L))`$, with a similar expression for the hole states. These expressions yield the basis for the electron-hole pair problem without Coulomb interaction nor magnetic field, with eigenvalues $`E_o=\mathrm{}^2k_{e}^{}{}_{}{}^{2}/2m_e+\mathrm{}^2k_{h}^{}{}_{}{}^{2}/2m_h`$. One can conveniently write the Coulomb interaction potential matrix elements using this noninteracting pair basis via Fourier transform integrals, as done in Ref. \[\] (see Appendix).

Similarly, can obtain the matrix elements of the Hamiltonian which depend on magnetic fields,

$$H_B=\frac{e}{2m_ec}𝐩_e𝐁\times 𝐫_e+\frac{e^2B^2}{8m_ec^2}r_{e}^{}{}_{}{}^{2}+\frac{e}{2m_hc}𝐩_h𝐁\times 𝐫_h+\frac{e^2B^2}{8m_hc^2}r_{h}^{}{}_{}{}^{2},$$
(21)
by using straightforward finite domain integrals of the basis set given above. The energies and eigenfunctions for the exciton are calculated by numerical diagonalization of the total Hamiltonian with magnetic fields obtained from the summation of all terms above. The wavefunctions are then represented as $`|\mathrm{\Psi }=_{n_en_hl_el_h}a_{n_en_hl_el_h}|n_e,n_h,l_e,l_h`$, where $`a_{n_en_hl_el_h}`$ are the coefficients calculated from the diagonalization. In turn, the mean electron-hole separation $`r_s`$ and the linear optical susceptibility can be calculated.

## III Results

We present here characteristic results of our calculations. As mentioned before, they are scalable for different materials, in terms of the Bohr radius of the exciton and its relation to the size (specially the width) of the ring. The parameters employed here describe GaAs, yielding an effective 2D Bohr radius of 6 nm. Figure 1 compares the exciton binding energies obtained for parabolic confinement (triangles) with those for a hard-wall confinement (diamonds), as function of the quantum ring width. Notice that $`E_b=E_{eh}^0E_{grnd}^{exciton}`$, where the first term is only the confinement ground state of the electron and hole, ignoring the Coulomb interaction. For relatively narrow rings, the exciton binding energy for the parabolic confinement model is much larger than for the hard-wall confinement, as expected. This difference is a reflection of the relative strength of the kinetic energy to Coulomb attraction increasing for the hard-wall case over the parabolic potential. In other words, even though we use nominally the same width ($`=2\sqrt{\mathrm{}/\mu \omega }`$ for the parabolic potential and $`2L`$ for the hard wall), the parabolic potential solutions are effectively less confined due to the finite amplitude ‘leaking’ out of the ring. It is also interesting to emphasize that for smaller ring widths, the resulting wavefunctions are increasingly elongated along the ring, and this is more the case for the hard-wall confinement. In contrast, the exciton binding energy appears larger for the hard-wall confinement, for widths larger than $`a_B^{2D}`$, due to the known poor-convergence of the parabolic basis used (e.g., see a detailed discussion of this problem in Ref. ). In the range of widths shown, both approaches have converged numerically to within a few percent everywhere, and at least an order of magnitude better for the lower two-thirds of the range).

The inset in Fig. 1 compares the ring results with those of a quantum dot (with parabolic confinement ) with equal confinement area (solid line). For wider ring systems all energies are basically equal, as the confinement potential is a weak perturbation to the Coulomb interaction between carriers, as one would expect, be it ring or dot. On the other hand, the exciton binding energies in the narrower rings are larger than in the dot case with the same area, a reflection of the anisotropic confinement in the ring: For the narrow rings, the circular symmetry of the 2D free exciton (either free or in the dot) is strongly affected, and the exciton elongates along the ring, as described above. We should mention that the curvature of the ring has not much effect on the ground state or binding energies for the dimensions considered here, similar to the experimental values.

To indicate the role of the Coulomb potential on the exciton characteristics, Fig. 2 shows the ground state energy of the electron-hole pair in the parabolic-confinement ring with (triangles) and without (diamond) Coulomb interaction. For a smaller ring width, the Coulomb contribution clearly increases, but not as fast as the confinement energy itself. The inset shows the electron-hole separation vs. the ring width. For small width, the confinement energy is clearly dominant in determining the electron-hole separation, rather than the Coulomb interaction term. For widths larger than $`4`$ nm, however, the electron-hole separation depends mostly on the Coulomb interaction term. Notice, however, that the rapid vanishing of $`r_s`$ for thin rings is somewhat of a biasing artifact produced by the basis functions we use in the parabolic ring. The necessary truncation of the basis appears to favor a circularly-symmetric exciton, clearly not the case in very thin rings.

Figure 3 shows the exciton binding energy and electron-hole separation versus the external magnetic field, for several ring widths, for a ring with middle radius $`r_o=20`$ nm. One can see that for larger values of the confinement energy (i.e., smaller widths), the effect of the magnetic fields is weaker, yielding the slowly changing curves at the top. However, for the larger widths, the dependence of the exciton binding energy on magnetic fields increases, resulting in the strong enhancement of the binding energies and decreasing exciton sizes with field. Notice that for larger values of the field, the exciton binding energy changes little as function of the ring width, showing that the confinement provided by the magnetic field is dominant. This is to be expected, given that the magnetic length, $`l_c=\sqrt{\mathrm{}/\mu \omega _c}`$ overtakes the exciton radius at about 18 T.

For the radius of the ring in this figure, $`r_o=20`$ nm, one expects ABE oscillations with a periodicity given by multiples of $`B\pi r_0^2/\varphi _0`$ (where $`\varphi _0=hc/e=4.14\times 10^7`$ gauss$``$cm<sup>2</sup> is the flux quantum), or a period $`\mathrm{\Delta }B3.3`$ T. We find no appreciable evidence of ABE oscillations in either the binding energy or the exciton effective size, $`r_s`$. This result suggests that it is likely that ABE exciton oscillations will not be seen in measurements of the ground state properties of the exciton.

Since there is a prediction that the ABE oscillations are to be found much more easily in the case of excited states,  we have also looked for them in the linear optical susceptibility of a quantum ring. Figure 4 shows a typical result for different values of the magnetic field. This curve represents all the possible transitions of this excitonic state which would be measurable via photoluminescence excitation measurements (PLE; while the first peak gives the PL response). The higher peaks, starting from the one at lowest energy, correspond to electron-hole excitations involving the heavy-hole exciton ground state and various center-of-mass replicas (i.e., increasing excitations of the center of mass degree of freedom, without altering the ground state of the relative coordinate). On the other hand, the smaller amplitude peaks (at shift $`310`$, 330 meV for $`B=0`$), correspond to internal excited states of the exciton (its relative coordinate). These peaks are strongly upshifted with magnetic field, as the diamagnetic effect for each charge carrier pushes all relative energies upwards as well, and clearly these excited states shift even faster. This behavior is qualitatively similar to the excitons in a quantum dot.  We should also point out that in addition to the overall upward shift due to the diamagnetic effect, the $`\chi `$ traces show no discernible superimposed ABE oscillations with magnetic field in any of the excited states. It would appear that the finite width of the system suppresses the ABE predicted for the 1D ring.

We should also mention that higher excited states are likely to exhibit the ABE effect, as per earlier work.  However, the parabolic or hard-wall confinement models used in this calculation lose validity, since non-parabolic corrections to the effective mass hamiltonian, as well as finite confinement potential effects would become more important. Consequently, a quantitative estimate of the anticipated ABE effects for high-lying states is less reliable, and more subject to specific values of parameters. A weaker confinement might also enhance ABE, although experiments in that regime would be harder to identify and characterize uniquely.

## IV Conclusions

We have shown that magnetic field has strong effects on excitons in a quantum ring, for both parabolic and hard-wall confinement potentials. Using direct matrix diagonalization techniques, we have shown that at least for rings currently realizable, the excitons behave to a great extent as those in quantum dots of similar dimensions: There are strong diamagnetic shifts and restructuring of the overall excitation manifold, large shifts of internal excitations, and reduction of the effective exciton size. On the other hand, the predicted ABE oscillations in the various physical characteristics (including binding energy and oscillator strength of transitions) of 1D excitons, are not found in this more realistic calculation. Although we anticipated that the predicted ABE effects would be much weaker (due to the finite transverse size of the rings and ring radii larger than the exciton size), we have not been able to detect any oscillation ‘remnant’, in any of the features we analyzed.

This negative result is due to either of two reasons, we surmise: One, the result of exponentially small ABE amplitudes (given the somewhat larger ring radii).  More likely, perhaps, this is the result of destructive interference of relatively many transverse eigenstates (mixed by the Coulomb interaction), each with its own different phase and amplitude. Notice that this is quite different for multiple-electron states in the ring, as predicted by theory, and recently seen in experiment. The difference in result from the case of only electrons in the ring to that of an exciton, indicates the predicted fragility of the effect, since in this system the net charge (and then coupling to the magnetic vector potential) is zero. In fact, the delicate nature of the ABE suggests that smaller and narrower rings are needed in experiment, which may make for more one-dimensional–like excitons. Following our discussion above, one could also expect ABE oscillations to be more important for higher-excited states, even if more challenging in experiments. Perhaps one could also think of a technique that explores differences, and therefore is able to couple only to a modulation of the hole population, for example, as a sensitive way to access these coherent ABE oscillations for low-lying states.

###### Acknowledgements.

We would like to thank T. Shahbazyan, M. Cobb, A. Lorke, and C. Trallero for helpful discussions. This work has been supported in part by the US Department of Energy grant no. DE–FG02–91ER45334.

##

The matrix elements of $`H_{rel}^{}`$, including the Coulomb interaction term, can be calculated by using the parabolic ring basis,

$`N^{}l_{CM}^{}n^{}l^{}|`$ $`H_{rel}^{}`$ $`|Nl_{CM}nl=\delta _{N^{},N}\delta _{l_{CM}^{},l_{CM}}\delta _{l^{},l}\{\sqrt{{\displaystyle \frac{n!n^{}!}{(n+|l|)!(n^{}+|l^{}|)!}}}`$ (22)
$`\times `$ $`{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \underset{j=0}{\overset{n^{}}{}}}{\displaystyle \frac{(n+|l|)!(n^{}+|l^{}|)!}{(k+|l|)!(j+|l^{}|)!}}(1)^{k+j}{\displaystyle \frac{1}{k!(nk)!j!(n^{}j)!}}`$ (23)
$`\times `$ $`{\displaystyle \frac{e^2}{ϵ}}\beta \mathrm{\Gamma }(|l|+k+j+{\displaystyle \frac{1}{2}})+{\displaystyle \frac{eB\mathrm{}l}{2\mu ^{}c}}\delta _{n^{},n}\}.`$ (24)
The matrix elements of the mixing terms between the center of mass and relative coordinates, in the limit $`rR`$, are of the form

$`N^{}l_{CM}^{}n^{}l^{}|H_{mix}|Nl_{CM}nl=`$ (25)
$`\times `$ $`\sqrt{{\displaystyle \frac{n!n^{}!}{(n+|l|)!(n^{}+|l^{}|)!}}}\mathrm{\hspace{0.17em}2}\alpha ^2\sqrt{{\displaystyle \frac{N!N^{}!}{(N+|l_{CM}|)!(N^{}+|l_{CM}^{}|)!}}}`$ (26)
$`\times `$ $`{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \underset{j=0}{\overset{n^{}}{}}}{\displaystyle \frac{(n+|l|)!(n^{}+|l^{}|)!}{(k+|l|)!(j+|l^{}|)!}}(1)^{k+j}{\displaystyle \frac{1}{k!(nk)!j!(n^{}j)!}}`$ (27)
$`\times `$ $`{\displaystyle \underset{K=0}{\overset{N}{}}}{\displaystyle \underset{J=0}{\overset{N^{}}{}}}{\displaystyle \frac{(N+|l_{CM}|)!(N^{}+|l_{CM}^{}|)!}{(K+|l_{CM}|)!(J+|l_{CM}^{}|)!}}(1)^{K+J}{\displaystyle \frac{1}{K!(NK)!J!(N^{}J)!}}`$ (28)
$`\times `$ $`[{\displaystyle \frac{eB\mathrm{}}{Mc\beta }}\mathrm{\Gamma }\left({\displaystyle \frac{2j+2k+|l|+|l^{}|+3}{2}}\right)\{[{\displaystyle \frac{1}{4\alpha }}\mathrm{\Gamma }\left({\displaystyle \frac{2J+2K+|l_{CM}|+|l_{CM}^{}|+3}{2}}\right)`$ (29)
$`+`$ $`{\displaystyle \frac{1}{4\alpha }}\left(2K+|l_{CM}|\right)\mathrm{\Gamma }\left({\displaystyle \frac{2J+2K+|l_{CM}|+|l_{CM}^{}|+1}{2}}\right)`$ (30)
$``$ $`{\displaystyle \frac{r_o}{4}}\mathrm{\Gamma }\left({\displaystyle \frac{2J+2K+|l_{CM}|+|l_{CM}^{}|+2}{2}}\right)+{\displaystyle \frac{r_o}{4}}\left(2K+|l_{CM}|\right)`$ (31)
$`\times `$ $`\mathrm{\Gamma }\left({\displaystyle \frac{2J+2K+|l_{CM}|+|l_{CM}^{}|}{2}}\right)\left]\right[\delta _{l,l^{}+1}\delta _{l_{CM},l_{CM}^{}1}\delta _{l,l^{}1}\delta _{l_{CM},l_{CM}^{}+1}]`$ (32)
$``$ $`[\delta _{l,l^{}+1}\delta _{l_{CM},l_{CM}^{}1}+\delta _{l,l^{}1}\delta _{l_{CM},l_{CM}^{}+1}]{\displaystyle \frac{l_{CM}}{4\alpha }}\mathrm{\Gamma }\left({\displaystyle \frac{2J+2K+|l_{CM}|+|l_{CM}^{}|+1}{2}}\right)\}`$ (33)
$``$ $`{\displaystyle \frac{\mu \omega _o^2r_o}{4\alpha \beta ^2}}\delta _{l,l^{}}\delta _{l_{CM},l_{CM}^{}}(|l|+k+j)!\{\mathrm{\Gamma }(|l_{CM}|+K+J+{\displaystyle \frac{1}{2}})`$ (34)
$``$ $`e^{\alpha ^2r_o^2}{\displaystyle \underset{jj=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\left(|l_{CM}|+K+J\frac{1}{2}\right)!}{\left(|l_{CM}|+K+J+\frac{1}{2}+jj\right)!}}\left(\alpha ^2r_o^2\right)^{|l_{CM}|+K+J+\frac{1}{2}+jj}\}].`$ (35)
Coulomb interaction in hard wall confinement,  is given by

$$H_{eh}(𝐪)=\frac{e^2}{ϵ}(2\pi )^2\frac{1}{𝐫_{eh}}e^{i𝐪𝐫_{eh}}𝑑𝐫_{eh}=2\pi \frac{e^2}{ϵ}\frac{1}{q},$$
(37)
and

$$\mathrm{\Psi }_e(𝐫_e)=\frac{1}{(2\pi )^2}\varphi _e(𝐪)e^{i𝐪𝐫_e}𝑑𝐪.$$
(38)
The Coulomb interaction matrix elements,

$$n_e^{}l_e^{}n_h^{}l_h^{}|H_{eh}|n_el_en_hl_h=\mathrm{\Psi }_e^{}(𝐫_e)\mathrm{\Psi }_h^{}(𝐫_h)H_{eh}(𝐫_{eh})\mathrm{\Psi }_e(𝐫_e)\mathrm{\Psi }_h(𝐫_h)𝑑𝐫_e𝑑𝐫_h,$$
(39)
are then rewritten as

$`n_e^{}l_e^{}n_h^{}l_h^{}|H_{eh}|n_el_en_hl_h`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^6}}{\displaystyle \varphi _e^{}(𝐪_e)\varphi _h^{}(𝐪_h)H_{eh}(𝐪)}`$ (40)
$`\times `$ $`\varphi _e(𝐪_e𝐪)\varphi _h(𝐪_h+𝐪)d𝐪_ed𝐪_hd𝐪,`$ (41)
where the Fourier transform integrals used for electron and hole wavefunctions, respectively, are

$`\phi _e(𝐪)`$ $``$ $`{\displaystyle \varphi _e^{}(𝐪_e)\varphi _e(𝐪_e𝐪)𝑑𝐪_e}`$ (42)
$`=`$ $`(2\pi ){\displaystyle _a^{a+2L}}\psi _e^{}(r_e)\psi _e(r_e)r_e𝑑r_e{\displaystyle _0^{2\pi }}e^{im\varphi _e}e^{iqr_ecos(\varphi _q\varphi _e)}𝑑\theta _e`$ (43)
$`=`$ $`(2\pi )^2(i)^me^{im\varphi _q}{\displaystyle _a^{a+2L}}\psi _e^{}(r_e)\psi _e(r_e)J_m(r_eq)r_e𝑑r_e,`$ (44)
where $`m=l_el_e^{}`$, and the definition of the Bessel function

$$J_m(qr)=\frac{1}{2\pi }_\pi ^\pi e^{im\varphi +izsin\varphi }𝑑\varphi $$
(45)
has been used. A similar expression for the hole wavefunctions gives the interaction matrix elements

$`n_e^{}l_e^{}n_h^{}l_h^{}|`$ $`H_{eh}`$ $`|n_el_en_hl_h=\delta _{l_e+l_h,l_e^{}+l_h^{}}{\displaystyle \frac{e^2}{ϵ}}{\displaystyle _a^{a+2L}}\psi _e^{}(r_e)\psi _e(r_e)r_e𝑑r_e`$ (46)
$`\times `$ $`{\displaystyle _a^{a+2L}}\psi _h^{}(r_h)\psi _h(r_h)r_h𝑑r_h{\displaystyle _0^{\mathrm{}}}J_{|m|}(r_eq)J_{|m|}(r_hq)𝑑q.`$ (47)
In order to evaluate the integrals in Eq. (47), we use

$`{\displaystyle _0^{\mathrm{}}}J_{|m|}(r_eq)J_{|m|}(r_hq)𝑑q=`$ $`{\displaystyle \frac{r_<^{|m|}}{r_>^{|m|+1}}}{\displaystyle \frac{\mathrm{\Gamma }(|m|+\frac{1}{2})}{\mathrm{\Gamma }(|m|+1)\mathrm{\Gamma }(\frac{1}{2})}}`$ (48)
$`\times `$ $`F(|m|+{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}};|m|+1;\left({\displaystyle \frac{r_<}{r_>}}\right)^2),`$ (49)
where $`r_>`$ ($`r_<`$) is the larger (smaller) of $`r_e`$ and $`r_h`$, and $`F`$ is a hypergeometric function.  Inserting (49) into (47), the interaction matrix elements can be written as

$`n_e^{}l_e^{}n_h^{}l_h^{}|`$ $`H_{eh}`$ $`|n_el_en_hl_h=\delta _{l_e+l_h,l_e^{}+l_h^{}}{\displaystyle \frac{e^2}{ϵ}}{\displaystyle \frac{\mathrm{\Gamma }(|m|+\frac{1}{2})}{\mathrm{\Gamma }(|m|+1)\mathrm{\Gamma }(\frac{1}{2})}}`$ (50)
$`\times `$ $`\{{\displaystyle _a^{a+2L}}\psi _e^{}(r_e)\psi _e(r_e){\displaystyle \frac{1}{r_e^{|m|}}}dr_e{\displaystyle _a^{r_e}}\psi _h^{}(r_h)\psi _h(r_h)r_h^{|m|+1}`$ (52)
$`\times F(|m|+{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}};|m|+1;\left({\displaystyle \frac{r_h}{r_e}}\right)^2)dr_h`$
$`+`$ $`{\displaystyle _a^{a+2L}}\psi _h^{}(r_h)\psi _h(r_h){\displaystyle \frac{1}{r_h^{|m|}}}𝑑r_h{\displaystyle _a^{r_h}}\psi _e^{}(r_e)\psi _e(r_e)r_e^{|m|+1}`$ (54)
$`\times F(|m|+{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}};|m|+1;\left({\displaystyle \frac{r_e}{r_h}}\right)^2)dr_e\}.`$
This greatly simplified expression for the interaction matrix elements is easily evaluated numerically. The total Hamiltonian matrix is diagonalized numerically, providing all the eigenvalues and eigenfunctions. To improve numerical convergence, we use the transformation of hypergeometric functions given by

$`F(|m|+{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}};|m|+1;z)={\displaystyle \frac{\mathrm{\Gamma }(|m|+1)}{\mathrm{\Gamma }(|m|+\frac{1}{2})\mathrm{\Gamma }(\frac{1}{2})}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(|m|+\frac{1}{2})_n(\frac{1}{2})_n}{n!^2}}`$ (55)
$`\times \left\{2\psi (n+1)\psi (|m|+n+{\displaystyle \frac{1}{2}})\psi (n+{\displaystyle \frac{1}{2}})ln(1z)\right\}(1z)^n,`$ (56)
where $`\psi `$ is the digamma function.

The matrix elements dependent on magnetic fields are given by

$`n_e^{}l_e^{}n_h^{}l_h^{}|H_B|n_el_en_hl_h=\delta _{l_e^{},l_e}\delta _{l_h^{},l_h}\{({\displaystyle \frac{l_h}{m_h}}{\displaystyle \frac{l_e}{m_e}}){\displaystyle \frac{\mathrm{}eB}{2c}}`$ (58)
$`\times {\displaystyle _a^{a+2L}}\psi _e^{}(r_e)\psi _e(r_e)r_edr_e{\displaystyle _a^{a+2L}}\psi _h^{}(r_h)\psi _h(r_h)r_hdr_h`$
$`+`$ $`{\displaystyle \frac{e^2B^2}{8m_ec^2}}{\displaystyle _a^{a+2L}}\psi _h^{}(r_h)\psi _h(r_h)r_h𝑑r_h{\displaystyle _a^{a+2L}}\psi _e^{}(r_e)\psi _e(r_e)r_{e}^{}{}_{}{}^{3}𝑑r_e`$ (59)
$`+`$ $`{\displaystyle \frac{e^2B^2}{8m_hc^2}}{\displaystyle _a^{a+2L}}\psi _e^{}(r_e)\psi _e(r_e)r_edr_e{\displaystyle _a^{a+2L}}\psi _h^{}(r_h)\psi _h(r_h)r_{h}^{}{}_{}{}^{3}dr_h\}.`$ (60)
The size of the exciton is here given by

$`r_{s}^{}{}_{}{}^{2}`$ $`=`$ $`\mathrm{\Psi }|r^2|\mathrm{\Psi }={\displaystyle \underset{n_en_hl_el_hn_e^{}n_h^{}l_e^{}l_h^{}}{}}a_{n_e^{}n_h^{}l_e^{}l_h^{}}^{}a_{n_en_hl_el_h}`$ (61)
$`\times `$ $`\{\delta _{l_e^{},l_e}\delta _{l_h^{},l_h}({\displaystyle _a^{a+2L}}\psi _h^{}(r_h)\psi _h(r_h)r_hdr_h{\displaystyle _a^{a+2L}}\psi _e^{}(r_e)\psi _e(r_e)r_{e}^{}{}_{}{}^{3}dr_e`$ (62)
$`+`$ $`{\displaystyle _a^{a+2L}}\psi _e^{}(r_e)\psi _e(r_e)r_edr_e{\displaystyle _a^{a+2L}}\psi _h^{}(r_h)\psi _h(r_h)r_{h}^{}{}_{}{}^{3}dr_h)`$ (63)
$``$ $`\delta _{l_e^{}+l_h^{},l_e+l_h}\delta _{l_e^{}+1,l_e}{\displaystyle _a^{a+2L}}\psi _e^{}(r_e)\psi _e(r_e)r_{e}^{}{}_{}{}^{2}𝑑r_e`$ (65)
$`\times {\displaystyle _a^{a+2L}}\psi _h^{}(r_h)\psi _h(r_h)r_{h}^{}{}_{}{}^{2}dr_h\}.`$