text:
# Positive rational solutions to 𝑥^𝑦=𝑦^{𝑚𝑥}: a number-theoretic excursion

## 1 Prologue

Late in the last millennium, the second author ran a seminar course for undergraduates which was intended to introduce them to problem-solving and question-asking in the context of mathematical research. He led them through the classic “difficult” equation

$$x^y=y^x,x,y>0,$$
(1)
whose solution is much easier than one would think at first glance. Solutions were sought, first over $``$, then over $``$ and finally over $``$. In the context of this seminar, it was natural to consider the variant equation

$$x^y=y^{2x},x,y>0,$$
(2)
which does not seem to have appeared in the literature. It turns out that there are solutions to (2) which do not fit the well-known parametric pattern of (1); c.f. (5) below. For example,

$$x=\left(\frac{4}{5}\right)^{128},y=\left(\frac{4}{5}\right)^{125}$$
(3)
is a solution to (2). This preposterous fact is trivial to verify: simply substitute (3) into (2), take logs and transpose:

$$\frac{2x}{y}=2\left(\frac{4}{5}\right)^3=\frac{128}{125}=\frac{\mathrm{log}x}{\mathrm{log}y}.$$

As we shall see, the ultimate explanation for this identity is that $`24^3=5^3+3`$. Upon discovering (3), the second author realized he was in over his head and contacted the first author. This paper is the result.

## 2 $`x^y=y^x`$

First, we review the familiar, but beautiful solution to (1), reserving historical references to the last paragraph of the section. We acknowledge the solutions $`x=y`$ and now let $`y=tx`$, $`t1`$, so that

$$x^{tx}=\left(tx\right)^xx^t=txx=t^{\frac{1}{t1}}.$$
(4)
The positive real solutions to (1) are thus

$$(x,y)=(u,u)\text{and}(x,y)=(t^{\frac{1}{t1}},t^{\frac{t}{t1}}).$$
(5)
(We might equally well have set $`y=x^r`$ in (1), and drawn essentially the same conclusion.) Since (1) implies that $`x^{1/x}=y^{1/y}`$ and since $`f\left(u\right)=u^{1/u}`$ increases on $`(1,e)`$ and decreases on $`(e,\mathrm{})`$, for each $`x(1,e)`$, there is exactly one $`y(e,\mathrm{})`$ so that (1) holds. In particular, the only integer solution to (1) is $`2^4=4^2`$.

Euler already noted that if $`t=1+\frac{1}{n}`$ for integral $`n`$, then (4) gives a rational solution to (1); namely,

$$x_n=\left(1+\frac{1}{n}\right)^n,y_n=\left(1+\frac{1}{n}\right)^{n+1}.$$
(6)
(Observe that, as $`n\mathrm{}`$, we have $`t1`$ and $`\left(x_n\right)`$ and $`\left(y_n\right)`$ increase and decrease monotonically to $`e`$, as is familiar from calculus.)

To show that these are the only rational solutions, we need an elementary lemma, whose proof, relying upon the Fundamental Theorem of Arithmetic, we omit:

###### Lemma 1.

Suppose $`a,b,m,n`$ are integers, with $`\mathrm{gcd}(a,b)=\mathrm{gcd}(m,n)=1`$ and $`b,n0`$. Then $`\left(\frac{m}{n}\right)^{a/b}`$ is rational if and only if $`m`$ and $`n`$ are $`|b|`$-th powers of integers.

Let us now proceed to find all rational solutions to (1). By symmetry, we may assume that $`y>x`$, so $`t>1`$. If $`x`$ and $`y`$ are rational, then so is $`t=y/x`$. Write $`t`$ in lowest terms as

$$t=\frac{p}{q}:=\frac{q+d}{q},$$

($`d,q>0`$), so that $`\frac{1}{t1}=\frac{q}{d}`$ and $`\frac{t}{t1}=\frac{q+d}{d}`$. With this substitution, we have

$$x=\left(\frac{q+d}{q}\right)^{q/d},y=\left(\frac{q+d}{q}\right)^{q/d+1}.$$

Since $`\mathrm{gcd}(d,q)=\mathrm{gcd}(q,q+d)=\mathrm{gcd}(d,q+d)=1`$, the integers $`q`$ and $`q+d`$ must both be $`\left|d\right|=d`$-th powers by Lemma 1. This causes no problem when $`d=1`$, of course, and, setting $`q=n`$, we recover (6). Suppose $`d>1`$, and write $`q=a^d`$, $`q+d=b^d`$ for positive integers $`a<b`$, so that $`b^da^d=d`$. Observe that

$$b^da^d\left(a+1\right)^da^d1+da1+d>d$$

for positive integers $`a`$ and $`b`$, and so there are no solutions with $`d>1`$. Finally, we remark that if $`n=r`$ is allowed to be negative in (6), then $`(x_r,y_r)=(y_{r1},x_{r1})`$:

$$\left(1\frac{1}{r}\right)^r=\left(1+\frac{1}{r1}\right)^r,\left(1\frac{1}{r}\right)^{r+1}=\left(1+\frac{1}{r1}\right)^{r1}.$$

Thus, by removing the restriction that $`n`$ be positive, we can eliminate the constraint $`y>x`$ in (6).

According to Dickson \[3, p.687\], the first reference to (1) was in a letter from D. Bernoulli to C. Goldbach, dated June 29, 1728. Bernoulli asserts, without proof, (see \[2, p.263\]) that this equation has only one solution in positive integers, and infinitely many rational solutions. The first person to write about (1) in detail was Euler (see \[4, pp.340–341\]. Euler made the substitution $`y=tx`$ and solved the equation over $`_+`$ and $`_+`$, and presented the rational solutions (6), without claiming that they were the only ones. Dickson mentions other writers who went over the same ground, and then, mysteriously, writes “\*A. Flechsenhaar and R. Schimmack discussed the rational solutions”. The asterisk means that the paper was “not available for report” (\[3, p.xxii\]). These papers appeared in 1911 and 1912 in the journal *Unterrichtsblätter für Mathematik und Naturwissenschaften*, and the authors found them without difficulty on the shelves of the magnificent Mathematics Library of the University of Illinois at Urbana-Champaign, then and now only 200 kilometers south of Dickson’s office. (For graduates of the US school system, a kilometer is a Canadian mile, attractively priced at $`0.621`$ American miles). In any event, Flechsenhaar  appears to earn credit as the first author to solve (1) over positive rationals. In 1967, Hurwitz  gave the first readily accessible analysis of the rational solutions; subsequent work on this equation, including generalizations to algebraic solutions, can be found in , , , ,  and .

## 3 $`x^y=y^{mx}`$

Let us now consider the generalization of Euler’s equation to

$$x^y=y^{mx}$$
(7)
where $`m>1`$ is a fixed positive integer. We will again restrict our attention to positive solutions $`(x,y)`$. If $`x=1`$ or $`y=1`$, then necessarily $`(x,y)=(1,1)`$. Supposing that $`x,y1`$ and taking logarithms in (7), we find that

$$\frac{m\mathrm{log}y}{y}=\frac{\mathrm{log}x}{x},$$

so $`xy`$. Write $`y=x^r`$, $`x1,r1`$. Then $`x^{r1}=mr`$. Therefore, the positive real solutions to (7) are given by $`x=y=1`$ and

$$x=\left(mr\right)^{\frac{1}{r1}}\text{ and }y=\left(mr\right)^{\frac{r}{r1}}.$$

We now restrict our attention to positive rational solutions. Since $`y/x=mr`$, it follows that $`r>0`$, and, further, since $`x`$ and $`y`$ are rational, we have $`r`$. Let us write $`r=\frac{a}{b}`$ where $`a,b`$ with $`\mathrm{gcd}(a,b)=1`$ and set $`k=\left|ab\right|_+`$. To have $`x`$, we require

$$x=\left(\frac{ma}{b}\right)^{\frac{b}{ab}}.$$

Suppose that $`\mathrm{gcd}(b,m)=d`$, and write $`b=db^{}`$ and $`m=dm^{}`$. Then $`\mathrm{gcd}(am^{},b^{})=1`$ and so we need $`am^{}`$ and $`b^{}`$ to be $`k`$-th powers of integers. If

$$am^{}=u^k\text{ and }b^{}=v^k,$$

it follows that

$$\left|u^kmv^k\right|=\left|am^{}mb^{}\right|=\left|\frac{am}{d}\frac{mb}{d}\right|=\frac{m}{d}k=m^{}k.$$
(8)
For a given positive integer $`m`$, we will classify the set $`S\left(m\right)`$ of positive rational solutions to (7) as follows: write

$$S\left(m\right)=\underset{k=0}{\overset{\mathrm{}}{}}S_k\left(m\right),$$

where $`S_k\left(m\right)`$ represents the set of solutions $`(x,y)`$ corresponding to equation (8), for $`m^{}`$ ranging over all positive integral divisors of $`m`$. Here, $`S_0\left(m\right)`$ denotes the solutions with $`x=y`$ (i.e. just the set $`(x,y)=(1,1)`$ for $`m>1`$).

The remainder of this paper will be devoted to analyzing $`S\left(m\right)`$. As a consequence, we will show how this set may be completely characterized for any given $`m`$. In case $`m=2`$ or $`3`$ (where all features of interest for general $`m`$ may in fact be observed), we have the following results:

###### Theorem 2.

If $`x`$ and $`y`$ are positive rational numbers for which $`x^y=y^{2x}`$ then either

(a) $`x=\left(2+\frac{2}{n}\right)^n`$ and $`y=\left(2+\frac{2}{n}\right)^{1+n}`$, and $`n`$, $`n0,1`$;

or

(b) $`(x,y)=(1,1),(2,16)`$ or $`((\frac{4}{5})^{128},(\frac{4}{5})^{125})`$.

and

###### Theorem 3.

If $`x`$ and $`y`$ are positive rational numbers for which $`x^y=y^{3x}`$ then either

(a) $`x=\left(3+\frac{3}{n}\right)^n`$ and $`y=\left(3+\frac{3}{n}\right)^{1+n}`$, and $`n`$, $`n0,1`$;

or

(b) $`x=(\frac{3w_n}{v_n})^{v_n^2}`$ and $`y=(\frac{3w_n}{v_n})^{3w_n^2}`$, $`0n`$;

or

(c) $`x=(\frac{w_n}{v_n})^{3w_n^2}`$ and $`y=(\frac{w_n}{v_n})^{v_n^2}`$, $`0n`$.

In (ii), $`v_n`$ and $`w_n`$ are the integers defined by $`v_n+w_n\sqrt{3}=(1+\sqrt{3})(2+\sqrt{3})^n.`$

Our proofs of Theorems 2 and 3 (perhaps rather surprisingly) involve techniques from transcendental number theory and Diophantine approximation. Moreover, they provide us with an opportunity to illustrate a fairly diverse grab-bag of methods for solving Diophantine problems.

## 4 The cases $`k=1`$ and $`k=2`$

The set $`S_1\left(m\right)`$ is easily computed. Indeed, we immediately find that $`k=1`$ implies that either

$$x=\left(m+\frac{m}{n}\right)^n,y=\left(m+\frac{m}{n}\right)^{1+n}$$

or

$$x=\left(m\frac{m}{n}\right)^n,y=\left(m\frac{m}{n}\right)^{1n}$$

for $`n`$ a positive integer (with $`n2`$ in the latter case). We note that these solutions correspond to the parametrized family (6) of solutions to Euler’s original equation and provide us with part (a) of Theorems 2 and 3.

If $`k2`$, the situation becomes rather more interesting, though the set $`S_2\left(m\right)`$ is also not too difficult to describe: it is either empty or infinite. The following lemmata provide sufficient conditions for the former to occur. As usual, for $`x`$, let $`\nu _2\left(x\right)`$ be the largest integer such that $`2^{\nu _2\left(x\right)}`$ divides $`x`$.

###### Lemma 4.

If $`\nu _2(m)`$ is odd and $`k`$ is even, then $`S_k(m)`$ is empty.

###### Proof.

We have $`\mathrm{gcd}(a,b)=1`$, and since $`k=\left|ab\right|`$ is even, $`a`$ and $`b`$ must both be odd. It follows that $`d=\mathrm{gcd}(b,m)`$ is odd as well. Write, as before, $`b=db^{}`$, $`m=dm^{}`$, $`am^{}=u^k`$, $`b^{}=v^k`$. Then $`\nu _2\left(m\right)=\nu _2\left(d\right)+\nu _2\left(m^{}\right)`$ is odd, hence so is $`\nu _2\left(am^{}\right)=\nu _2\left(u^k\right)=k\nu _2\left(u\right)`$, which contradicts $`k`$ being even. ∎

###### Lemma 5.

If $`m=2^{\nu _2(m)}m_1`$ is a positive integer for which $`m_11\text{ (mod }4\text{)}`$, then $`S_2(m)`$ is empty.

###### Proof.

If $`\nu _2\left(m\right)`$ is odd, this is a special case of Lemma 4. Since $`k=2`$ is even, we may conclude as before that $`a`$ and $`b`$ are both odd. Suppose $`\nu _2\left(m\right)`$ is even, say $`\nu _2\left(m\right)=2t`$. It follows that

$$u^22^{2t}m_1v^2=\pm 2^{2t+1}\frac{m_1}{d}$$

where $`u`$ and $`v`$ are coprime and $`d`$ divides $`m`$. Since $`d|b`$, it follows that $`d`$ is odd, and so $`d|m_1`$. We thus have that $`2^t|u`$, say $`u=2^tu_1`$, whereby

$$u_1^2m_1v^2=\pm 2\frac{m_1}{d}.$$
(9)
Since the right hand side of this equation is even, $`u_1`$ and $`v`$ have the same parity, and, since $`u`$ and $`v`$ are coprime, are both necessarily odd. Since $`m_11\text{ (mod }4\text{)}`$, this implies that the left hand side of (9) is divisible by $`4`$. Since the right hand side of this equation is congruent to $`2`$ modulo $`4`$, this yields the desired contradiction. ∎

It is a well-known fact that if $`m`$ is a positive nonsquare integer and $`c`$ is a nonzero integer, then a single solution in positive integers $`x`$ and $`y`$ to the equation $`x^2my^2=c`$ implies the existence of infinitely many such solutions. In fact, one can find a finite collection of pairs of positive integers, say

$$(x_1,y_1),(x_2,y_2),\mathrm{},(x_r,y_r),$$

with $`x_i^2my_i^2=c`$ for $`i=1,2,\mathrm{},r`$, such that every solution in positive integers $`(x,y)`$ to the equation $`x^2my^2=c`$ satisfies

$$x+y\sqrt{m}=\left(x_i+y_i\sqrt{m}\right)\left(u_1+v_1\sqrt{m}\right)^k$$
(10)
where $`k`$ is a nonnegative integer, $`i\{1,2,\mathrm{},r\}`$ and $`(u_1,v_1)`$ is the smallest positive integer solution to the equation $`u^2mv^2=1`$. The integer $`r`$ here depends upon $`c`$ and, potentially, upon $`m`$.

It follows that if $`S_2\left(m\right)`$ is not empty, then it is infinite. From the theory of Pell equations (or Fermat-Pell equations if one likes) these correspond to elements of finitely many recurrence sequences (see Nagell  for a nice, affordable exposition of such matters). This fact is also a pretty easy consequence of (10). We note that consideration of the case $`m=39`$ demonstrates that the above lemmata do not in fact provide necessary conditions for $`S_2\left(m\right)`$ to be empty. A routine computation shows that for $`2m50`$, $`S_2\left(m\right)`$ is infinite precisely for

$$m\{3,7,11,12,15,19,23,27,28,31,35,43,44,47,48\}.$$

It is also not hard to provide sufficient conditions for $`S_2\left(m\right)`$ to be nonempty. An example of such a result is as follows:

###### Lemma 6.

If $`p3\text{ (mod }4\text{)}`$ is prime then $`S_2(p)`$ is infinite.

###### Proof.

From a venerable (after citing Euler, we can hardly call this old!) result of Petr , precisely one of

$$x^2py^2=2\text{ or }x^2py^2=2$$

is solvable in integers $`x`$ and $`y`$. This implies, in either case, that $`S_2\left(p\right)`$ is infinite. ∎

We leave to our gentle reader the (rather nontrivial) task of deriving necessary conditions for $`S_2\left(m\right)`$ to be nonempty. As we shall see in the next section, the behavior of $`S_k\left(m\right)`$ is dramatically different when $`k3`$.

## 5 Thue equations

A famous theorem of the Norwegian mathematician Axel Thue  asserts, if $`\theta `$ is an algebraic number of degree $`n3`$ and $`ϵ>0`$, that the inequality

$$\left|\theta \frac{x}{y}\right|<\frac{1}{\left|y\right|^{n/2+1+ϵ}}$$

has at most finitely many solutions in integers $`x`$ and $`y`$ with $`y0`$. Note that, if $`\theta =\sqrt[k]{m}`$, then we have the algebraic identity

$$\left|x^kmy^k\right|=y^k\left|\theta \frac{x}{y}\right|\left(\left(\frac{x}{y}\right)^{k1}+\theta \left(\frac{x}{y}\right)^{k2}+\mathrm{}+\theta ^{k1}\right).$$

If $`m`$ is not a perfect $`k`$-th power and $`c0`$ is an integer, then it follows that the equation

$$x^kmy^k=c$$
(11)
has at most finitely many solutions in integers $`x`$ and $`y`$. Such equations are nowadays termed Thue equations (a nice exposition of this active area of research may be found in the book of Fel’dman and Nesterenko ).

In our situation, this immediately implies that $`S_k\left(m\right)`$ is finite for each fixed integer $`k3`$. As we shall see in the next two sections, it is possible to “effectively” determine each such $`S_k\left(m\right)`$, (and to derive an upper bound upon $`k`$ that depends solely upon $`m`$.

## 6 Linear forms in logarithms

The sets $`S_k\left(m\right)`$ are determined by equalities of the shape

$$u^kmv^k=\pm \frac{m}{d}k$$

for $`u`$ and $`v`$ positive integers. In all cases, we thus have

$$\left|u^kmv^k\right|mk.$$

It follows that

$$\left|m^1\left(\frac{u}{v}\right)^k1\right|\frac{k}{v^k}$$

and hence the linear form in logarithms $`\left|k\mathrm{log}\left(u/v\right)\mathrm{log}\left(m\right)\right|`$ is “small”. A classical result of Gel’fond  (extending his work on Hilbert’s 7th problem) indicates that, for any given nonzero algebraic numbers $`\alpha _1,\alpha _2,\beta _1`$ and $`\beta _2`$, with $`\mathrm{log}\alpha _1`$ and $`\mathrm{log}\alpha _2`$ linearly independent over the rationals, we have

$$\left|\beta _1\mathrm{log}\alpha _1\beta _2\mathrm{log}\alpha _2\right|0$$

and, moreover, provides lower bounds for such a form. Applying state-of-the-art versions of these bounds, say those due to Laurent, Mignotte and Nesterenko , we may conclude, as in Theorem 3 of Mignotte , that, in our situation, $`k<600`$ (if $`m=2`$) and, more generally, as in Theorem 2 of , that

$$k<10676\mathrm{log}m.$$

It follows that

$$\underset{k=3}{\overset{\mathrm{}}{}}S_k\left(m\right)$$

is finite. In the remaining sections, we will describe a strategy for explicitly determining this set and illustrate it in the cases $`m=2`$ and $`m=3`$.

## 7 Solving the remaining equations

For small values of $`k`$, it is possible to use standard computational techniques based upon lower bounds for linear forms in logarithms, combined with lattice basis reduction, to solve the Thue equations that occur (a reasonably accessible book which covers this field is that of Smart ). If, however, $`k`$ is moderately large, this becomes a computational problem due to the difficulty in finding systems of independent units in the ring of integers of $`\left(\sqrt[k]{m}\right)`$. (If one has cycles to burn, try, for instance, to find the fundamental units in, say, $`\left(\sqrt[41]{2}\right)`$.)

In our situation, though, we are able to find local obstructions to solvability for virtually all values of $`k`$ under consideration, obviating the need for extensive computations. For the sake of simplicity, let us restrict our attention to the case $`m=2`$ (where, as mentioned previously, we may assume $`k<600`$). Here the equations to be studied are of the shape

$$x^k2y^k=\pm 2^\delta k$$

where $`\delta \{0,1\}`$. Lemma 4 also allows us to suppose that $`k`$ is odd. For each such $`k`$, we consider primes of the form $`p=2nk+1`$ for $`n`$. For these $`p`$, there are at most $`\left(2n+1\right)^2`$ values of $`x^k2y^k`$ modulo $`p`$. If none of these are congruent to $`\pm k`$ or $`\pm 2k`$ modulo $`p`$, we have found a local obstruction to solvability and can thus conclude that $`S_k\left(2\right)`$ is empty. If, for example, $`k=13`$, we consider the above equation(s) modulo $`53`$. Noting that $`x^{13}0,\pm 1,\pm 23\text{ (mod }53\text{)}`$, it follows that

$$x^{13}2y^{13}0,\pm 1,\pm 2,\pm 3,\pm 6,\pm 7,\pm 8,\pm 16,\pm 21,\pm 23,\pm 25\text{ (mod }53\text{)}$$

and hence may conclude that $`S_{13}\left(2\right)`$ is empty. Similar arguments suffice to eliminate all equations (for $`k<600`$) but

$$x^32y^3=\pm 3,x^32y^3=\pm 6,x^52y^5=\pm 10,$$

$$x^72y^7=\pm 7,x^72y^7=\pm 14\text{ and }\text{ }x^{11}2y^{11}=\pm 22$$

(to verify this, the alert reader may wish to write her own piece of code; the cognizant author used Maple ™, being careful to avoid our silico-aceric chum’s rather dubious “msolve” routine).

Similarly, when $`m=3`$, we easily deal with all the equations encountered, with the exceptions of

$$x^33y^3=\pm 3,x^53y^5=\pm 15\text{ and }\text{ }x^73y^7=\pm 21$$

(though we have many more $`k`$ to treat – up to $`k=11728`$, including even values).

## 8 Endgame

As mentioned in the last section, there are various techniques from Diophantine approximation that may be used to solve these remaining equations. The symbolics package Kash , for example, has a built-in Thue solver that can handle all the equations under consideration in a matter of minutes on a Sun Ultra. In any case, nowadays it is a routine matter to verify that of the equations that have so far evaded our net, only those with $`k=3`$ possess solutions, corresponding to $`5^324^3=3`$, $`2^321^3=6`$ and $`3^332^3=3`$. By way of example, to solve the Diophantine equations $`\left|x^32y^3\right|\{3,6\}`$, one can appeal to the inequality

$$\left|x^32y^3\right|\sqrt{\left|x\right|},$$

valid for all integers $`x`$ and $`y`$ (see e.g. ).

The equation $`5^324^3=3`$ maps back to the solution $`(x,y)`$ to $`x^y=y^{2x}`$ given in (3). Similarly, $`2^321^3=6`$ yields the solution $`(x,y)=(2,16)`$, also to $`x^y=y^{2x}`$. The equation $`3^332^3=3`$ which potentially yields a solution to $`x^y=y^{3x}`$, leads to $`a=27,b=24`$, contradicting the coprimality of $`a`$ and $`b`$.

Similar arguments suffice to completely solve $`x^y=y^{mx}`$ for, with a modicum of computation, all values of $`m`$ up to $`40`$ or so. As $`m`$ increases, we are faced with the prospect of solving Thue equations of higher and higher degree, an apparently formidable task. We note that $`S_k\left(m\right)`$ can be nonempty for arbitrarily large $`k`$ (for example, this is always the case for $`S_k\left(2^kk\right)`$ if $`k`$ is odd).

## 9 Postscript : Mutterings on local-global principles

A standard heuristic employed in the field of Diophantine problems is that an equation should be solvable over $``$ precisely when it is solvable over $``$ and over $`p`$-adic fields $`_p`$, for all primes $`p`$. Such local-global or Hasse principles are known to be true in various settings, but false in general.

One of the first instances where local-global principles were shown to fail occurs in work of Skolem  of 1937, where he demonstrated that the equation $`x^32y^3=47`$ can be solved modulo $`p^\alpha `$ for every prime power $`p^\alpha `$, but has no solution over the integers.

In our analysis, we encounter numerous equations for which a like conclusion holds, of degree up to (at least) $`19`$. Such is the case, for instance, for the equations

$$x^35y^3=15,x^52y^5=10,x^53y^5=15,x^72y^7=7,$$

$$x^72y^7=14,x^73y^7=21,x^{11}2y^{11}=22.$$

To actually prove that the Hasse principle fails in these instances is an interesting exercise which we leave to our loyal (and, at this stage, perhaps rather fatigued) reader!

## 10 Acknowledgments

We thank Gerry Myerson for his comments on an earlier version of this paper. Some computations in this paper were performed by using Maple™. Maple and Maple 8 are registered trademarks of Waterloo Maple Inc.


format: text

text:
# Allowed and forbidden transitions in artificial hydrogen and helium atoms

\[

\]

The strength of radiative transitions in atoms is governed by selection rules . Spectroscopic studies of allowed transitions in hydrogen and helium provided crucial evidence for the Bohr’s model of an atom. “Forbidden” transitions, which are actually allowed by higher-order processes or other mechanisms, indicate how well the quantum numbers describe the system. We apply these tests to the quantum states in semiconductor quantum dots (QDs), which are regarded as artificial atoms. Electrons in a QD occupy quantized states in the same manner as electrons in real atoms . However, unlike real atoms, the confinement potential of the QD is anisotropic, and the electrons can easily couple with phonons of the material. Understanding the selection rules for such QDs is an important issue for the manipulation of quantum states. Here we investigate allowed and forbidden transitions for phonon emission in one- and two-electron QDs (artificial hydrogen and helium atoms) by electrical pump-and-probe experiments, and find that the total spin is an excellent quantum number in artificial atoms. This is attractive for potential applications to spin based information storage.

The QD we study is located in a circular pillar (diameter of 0.5 $`\mu `$m) fabricated from an AlGaAs/InGaAs heterostructure (see Figs. 1a and 1b) . Electrons can be injected from, and collected by, the source (s) and drain (d) electrodes through asymmetric tunneling barriers, whose tunneling rates are $`\mathrm{\Gamma }_s`$ (3 ns)-1 and $`\mathrm{\Gamma }_d`$ (100 ns)-1 for the barrier nearest to the source and the drain, respectively. The number of electrons in the dot, $`N`$, can be controlled exactly by applying a voltage, $`V_g`$, to the surrounding gate electrode (g). Electrons are confined in an In0.05Ga0.95As quantum well (thickness $`a`$ = 12 nm) in the vertical ($`z`$) direction, and by a two-dimensional harmonic potential in lateral ($`x`$ and $`y`$) direction appropriate for small $`N`$. The first few electrons occupy the 1s and 2p orbitals associated with the lateral confinement. Since our QD does not have circular symmetry , orbital degeneracy is lifted even at zero magnetic field, and only a two-fold spin degeneracy is expected. For this non-circular (approximately elliptic) QD, we still use 1s, 2p, … to label the orbitals for convenience. This non-circularity does not affect much our discussion or the underlying physics. The following experiments are carried out at a temperature, $`T`$, of $``$100 mK, unless otherwise stated.

First we investigate the $`N=`$ 1 QD (artificial hydrogen), in which single electron occupies the 1s orbital (the ground state) or the 2p orbital (the first excited state). The energy spectrum of these states can be obtained by tunneling spectroscopy . The current, $`I_{sd}`$, increases stepwise with increasing $`V_g`$ (giving peaks in $`dI_{sd}/dV_g`$) each time an empty dot state enters the transport window. Color plot of $`dI_{sd}/dV_g`$ vs $`V_g`$ traces taken as a function of magnetic field ($`B`$) applied in the $`z`$ direction are shown in Fig. 1c. The peak spacings within a given stripe can be related to the energy spacings between corresponding states. The energy spacing between the 1s and 2p states, $`\epsilon _{1s2p}`$, of the $`N`$ = 1 QD deduced from the first current stripe is plotted in Fig. 2e. The lateral confinement of our QD can be approximated by two orthogonal harmonic potentials in the $`x`$ and $`y`$ plane . The characteristic confinement energies, $`\mathrm{}\omega _x`$ 2.5 meV and $`\mathrm{}\omega _y`$ 5.5 meV ($`\mathrm{}h/2\pi `$ is the Planck’s constant), are obtained by fitting $`\epsilon _{1s2p}`$ (see solid line in Fig. 2e). The Zeeman splitting of the electron spin is small ($`E_Z`$ 0.1 meV at $`B=`$ 5 T), and is unresolved in our experiments.

We now focus on the energy relaxation from the 2p state to the 1s state in the $`N=1`$ QD, which changes the electron’s orbital momentum but preserves the spin. Electrical pump-and-probe experiments are performed by applying a gate voltage of rectangular-shaped time dependence, $`V_g(t)`$, which switches between $`V_l`$ and $`V_h`$ as shown in Fig. 2a, to the gate. Experimental details are given in Refs. 9 & 10. First, the $`N`$ = 0 QD is prepared during the low-phase of the pulse ($`V_g=V_l`$, see Fig. 2b). The period, $`t_l`$ = 100 - 200 ns, is made sufficiently long to ensure that both the 1s and 2p states are empty. When the pulse is switched on ($`V_g=V_h`$, see Fig. 2c) such that only the 2p state is located in the transport window, an electron can be injected into the 2p state from the source (pump) with a time constant, $`\mathrm{\Gamma }_s^1`$ 3 ns. The electron can only escape to the drain (probe) more slowly, with a time constant, $`\mathrm{\Gamma }_d^1`$ 100 ns. However, this escape process can be interrupted by the relaxation into the 1s ground state. Thus, the current contains information about the relaxation lifetime, $`\tau _{1s2p}`$. We measure the averaged dc current, $`I_p`$, under the application of the pulse train. Figure 2d shows how this current changes with the pulse length, $`t_h`$. $`I_p`$ is then converted into an average number of tunneling electrons per pulse, $`n_t=I_p(t_h+t_l)/e`$ ($`e`$ is the elementary charge). From a detailed analysis of the rate equations including all possible tunneling processes for the $`N=1`$ QD, we find $`n_t\mathrm{\Gamma }_d\tau _{1s2p}[1\mathrm{exp}(t_h/\tau _{1s2p})]`$ under the condition, $`\mathrm{\Gamma }_s^1\tau _{1s2p}<\mathrm{\Gamma }_d^1`$, required for the relaxation time measurement . We made sure that this condition is satisfied in all measurements for the $`N=1`$ QD. The relaxation time thus estimated from the rise time of $`n_t`$ is $`\tau _{1s2p}`$ 10 ns for the case in Fig. 2d. The small saturation value of $`n_t`$ 0.02 indicates very efficient relaxation.

We discuss the origin of the relaxation in artificial hydrogen. For an energy spacing in the meV region, and at low temperature (thermal energy $``$10 $`\mu `$eV $`\epsilon _{1s2p}`$), spontaneous emission of a phonon, rather than a photon, dominates the relaxation process . The energy spacing coincides with the acoustic-phonon energy in the linear dispersion regime. Because of the discrete energy of the states, the relaxation involves emission of a phonon of energy equal to $`\epsilon _{1s2p}`$ (with corresponding phonon wavelength, $`\lambda _{1s2p}`$). In our experiment, $`\epsilon _{1s2p}`$, and hence $`\lambda _{1s2p}`$, vary with $`B`$ as shown in Figs. 2e and 2f. Note also that the characteristic sizes $`l_x`$ and $`l_y`$ for the lateral dimensions of the QD decrease with increasing $`B`$. The strength of the electron-phonon interaction is expected to be suppressed for $`\lambda _{1s2p}`$ smaller than the characteristics size of the QD (phonon bottleneck effect ). Therefore, the $`B`$ dependence of $`\tau _{1s2p}`$ in Fig. 2g is because the phonon emission is suppressed (i.e. $`\tau _{1s2p}`$ increases) with decreasing $`B`$ when $`\lambda _{1s2p}`$ becomes shorter than $`a`$, $`l_x`$, and $`l_y`$.

In order to be quantitative, we calculate the phonon emission rate from Fermi’s golden rule including both deformation and piezoelectric types of coupling with standard GaAs material parameters . For simplicity, the calculation is done for a circular dot, whose effective confinement energy is $`\mathrm{}\omega _{eff}=\mathrm{}\sqrt{\omega _x\omega _y(1+\omega _c^2/(\omega _x+\omega _y)^2)}`$, where $`\omega _c`$ is the cyclotron frequency. This assumption is reasonable because Coulomb interactions in an elliptic QD just scale with $`\mathrm{}\omega _{eff}`$ . As shown by the solid line in Fig. 2g, we find good agreement with the data. The difference by about a factor of 2 or 3 might come from the assumptions about the confinement potential and uncertainties in the material parameters. Thus, the fast energy relaxation in $`N=`$ 1 QD can be understood by spontaneous emission of a phonon. It should be noted that no clear selection rule for orbital momentum is expected, because the electron-phonon interaction cannot be approximated by a dipole interaction, and because of the anisotropic confinement potential.

In contrast, the relaxation time is remarkably different for a $`N=`$ 2 QD (artificial helium). At low magnetic fields (see second stripe in Fig. 1c for $`B`$ $`<`$ 2.5 T), the many-body ground state is a spin-singlet (labeled S) with two antiparallel-spin electrons occupying the 1s orbital, while the first excited state is a spin-triplet (labeled T) with two parallel-spin electrons, one each occupying the 1s and 2p orbitals . Because of direct Coulomb and exchange interactions, the energy spacing between the two states, $`\epsilon _{ST}`$ ($``$0.6 meV at $`B=`$ 0 T), is smaller than $`\epsilon _{1s2p}`$. Energy relaxation from the first excited state (T) to the ground state (S) not only involves the same change in orbital momentum as that in the $`N=`$ 1 QD, but also requires a spin-flip because of Pauli exclusion (see inset of Fig. 3e). A simple phonon-emission transition from the spin-triplet to the spin-singlet is forbidden by spin conservation.

We now investigate to what degree this transition is “forbidden”. The simple rectangular pulse technique employed for the $`N=`$ 1 QD is not useful for this $`N=`$ 2 QD transition, because the relaxation lifetime, $`\tau _{ST}`$, is always beyond the measurable range ($`\tau _{ST}>\mathrm{\Gamma }_d^1`$ 100 ns) . Instead, we subject the QD to a double-step voltage-pulse, in which $`V_g`$ is switched between three voltages, $`V_l`$, $`V_h`$ and $`V_m`$ as shown in Fig. 3a. First, when $`V_g=V_l`$ (Fig. 3b), the $`N=`$ 1 QD is prepared during a sufficiently long period, $`t_l`$ = 100 ns. When $`V_g`$ is suddenly increased to $`V_h`$ (Fig. 3c), an electron can enter to create the $`N=`$ 2 triplet state within the interval $`\mathrm{\Gamma }_s^1=`$ 3 - 7 ns. $`V_g=V_h`$ for the duration $`t_h=`$ 100 ns \- 100 $`\mu `$s, which is much longer than $`\mathrm{\Gamma }_s^1`$. The triplet state may suffer a relaxation process during this time. When $`V_g`$ is changed to $`V_m`$ (Fig. 3d), an electron in the triplet state can tunnel out to the drain, if the triplet state is not yet relaxed to the singlet state. We set the period, $`t_m`$ = 300 ns, to be longer than $`\mathrm{\Gamma }_d^1`$ so we can read out the signal. We repeatedly apply the double-step pulse (effectively $``$107 times) to obtain a reliable current $`I_p`$, and evaluate the average number of tunneling electrons, $`n_t=I_p(t_l+t_h+t_m)/e`$. Since the current measures the unrelaxed electron number, $`n_t=A\mathrm{exp}(t_h/\tau _{ST})`$ for the condition $`\mathrm{\Gamma }_s^1\tau _{ST}`$ (no upper limit in principle), where $`A1`$ is approximately the ratio of $`\mathrm{\Gamma }_s^1`$ for the triplet to that for the singlet. Figure 3e shows a typical measurement of $`n_t`$ at 0 T, indicating a relaxation time of $`\tau _{ST}`$ 200 $`\mu `$s. This relaxation time is 4 to 5 orders of magnitude longer than that observed in the $`N=`$ 1 QD.

We find no clear $`B`$ dependence of $`\tau _{ST}`$ (always longer than 100 $`\mu `$s), at least for the energy spacings between 0.6 meV at $`B=`$ 0 T and 0.24 meV at $`B=`$ 2 T (not shown). We also investigate the temperature dependence of $`\tau _{ST}`$ (Fig. 3f). No clear change is observed up to 0.5 K. $`\tau _{ST}`$ decreases above 0.5 K, where thermal excitation from the QD to the electrodes becomes important.

On the other hand, we do find that $`\tau _{ST}`$ strongly depends on the high gate-pulse voltage $`V_h`$, during which relaxation takes place, as shown in Fig. 3g. Although $`V_h`$ is swept deep into the $`N=`$ 2 Coulomb blockade region (-1.27 V $`<V_g<`$ -1.16 V), $`\tau _{ST}`$ decreases rapidly at $`V_h`$ -1.18 V. This $`V_h`$ dependence implies a strong influence of the source and drain electrodes. Even though the Coulomb blockade is robust in the suppression of transport, higher-order tunneling processes can contribute to the relaxation and change $`\tau _{ST}`$. An electron in the dot can be replaced with an electron of opposite spin from the electrodes (see inset of Fig. 3g). This results in energy loss in the QD, whereas the electrode gains the same energy. This inelastic cotunneling rate, $`\tau _{co}^1`$, is estimated by considering second-order tunneling processes . For the relaxation mechanisms considered here, the $`N=`$ 2 QD can relax virtually through $`N=`$ 1 or $`N=`$ 3 intermediate states. Note that this process does not cause a net current even at a finite voltage of $`e|V_{sd}|<\epsilon _{ST}`$ . Assuming $`V_{sd}=`$ 0 V and zero temperature for simplicity, we obtain $`\tau _{co}^1=\epsilon _{ST}(\mathrm{}\mathrm{\Gamma }_s+\mathrm{}\mathrm{\Gamma }_d)^2(\mathrm{\Delta }_1^1+\mathrm{\Delta }_3^1)^2/h`$. Here, $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_3`$, respectively are the energies required to excite the initial $`N=`$ 2 triplet state to the $`N=`$ 1 and 3 intermediate states, as indicated in Fig. 3c. We can extract $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_3`$ from $`V_h`$, and the values are shown in Fig. 3g. The solid line shows $`\tau _{co}`$, the relaxation time due to cotunneling, calculated with experimentally deduced parameters ($`\epsilon _{ST}=`$ 0.6 meV, and $`(\mathrm{\Gamma }_s+\mathrm{\Gamma }_d)^1=`$ 7 ns). Clearly the observed relaxation time can be well understood by inelastic cotunneling.

Our observations for $`N=`$ 1 and 2 QDs can be compared with real atoms . The transition from the 2p state to the 1s state in atomic hydrogen is allowed by photon emission (the Lyman $`\alpha `$ transition line), while that in artificial hydrogen is allowed by phonon emission. The transition from the spin-triplet state to the spin-singlet state is forbidden by conservation of the total spin for both atomic helium and artificial helium. Moreover, these “forbidden” transitions can occur by collisions with electrons for the helium atom, and by cotunneling for the $`N=`$ 2 QD. Spin conservation also applies to any transition in light atoms, and probably likewise to few-electron QDs. The difference between the allowed and forbidden transitions leads to more than 11 orders of magnitude difference in the relaxation times for real hydrogen and helium atoms. Our observation of 4 to 5 orders of magnitude difference in our $`N=1`$ and 2 artificial atoms is not as high, but is still surprisingly large. Note that this difference would become larger if cotunneling in QDs can be suppressed by using thicker tunneling barriers. Very importantly, the large difference between $`\tau _{1s2p}`$ and $`\tau _{ST}`$ originates from the fact that other effects, such as spin-orbit and hyperfine interactions , must have only a weak effect on the breaking of the “forbidden” symmetries. We now discuss how small these hidden contributions are by focusing on the spin-orbit interactions.

Spin-orbit interactions are predicted to give the dominant contribution to spin relaxation in GaAs QD systems , although this is still an extremely small effect. For simplicity, we consider the spin-orbit interaction energy, $`\mathrm{\Delta }_{so}`$, only for coupling between the 1s and 2p orbitals, but including all effects which mix spin and orbital degrees of freedom. Simple perturbation theory predicts that the relaxation time from the triplet to the singlet is given by $`\tau _{ST,so}(\epsilon _{ST}/\mathrm{\Delta }_{so})^2\tau _{phonon}(\epsilon _{ST})`$. Here, $`\tau _{phonon}^1(\epsilon _{ST})`$ is the phonon emission rate at the phonon energy, $`\epsilon _{ST}`$, and we know that $`\tau _{phonon}`$ is well accounted for by the electron-phonon interaction. Therefore, we can deduce an upper bound of $`\mathrm{\Delta }_{so}<`$ 4 $`\mu `$eV from our observations ($`\tau _{ST}>`$ 200 $`\mu `$s). This value is close to the spin-orbit induced spin splitting energy ($``$ 2.5 $`\mu `$eV) observed in a GaAs two-dimensional electron gas system . Note that spin-orbit interactions are significantly enhanced in nanoparticles, e.g. copper, probably because of impurities or interfaces , but our vertical semiconductor QD are largely free of these undesirable factors.

Our experiments indicate that the spin degree of freedom in QDs is well isolated from the orbital degree of freedom. This is particularly attractive for applications to spin memories and spin quantum bits (qubits) . For a simple scheme involving just a single-electron spin in a magnetic field, the spin-orbit interactions can degrade the energy relaxation time ($`T_1`$) of a spin qubit. We estimate the dominant contribution, $`T_{1,so}`$, using a perturbative approach, $`T_{1,so}(\epsilon _{1s2p}/\mathrm{\Delta }_{so})^2\tau _{phonon}(\epsilon _Z)`$. Since $`\mathrm{\Delta }_{so}<`$ 4 $`\mu `$eV, this yields $`T_{1,so}>`$ 1 ms for a Zeeman splitting $`\epsilon _Z`$ 0.1 meV and $`\epsilon _{1s2p}`$ 1.2 meV at $`B`$ $``$ 5 T ($`T_{1,so}>`$ 100 $`\mu `$s at $`B`$ $``$ 9 T). This $`T_{1,so}`$ is thus comparable to that obtained by electron-spin-resonance for donor states in GaAs , and is much longer than the time required for typical one- and two-qubit operations . Note that small spin-orbit interactions are also desirable with respect to the dephasing time ($`T_2`$) of a spin qubit. Our results therefore encourage further research in the use of the spin degree of freedom in QDs.

Acknowledgement

We thank G. E. W. Bauer, T. Honda, T. Inoshita, A. V. Khaetskii, L. P. Kouwenhoven for discussions and help.

Correspondence and requests for materials should be addressed to T.F. (e-mail: fujisawa@will.brl.ntt.co.jp).


text:
# Robust Ising Gates for Practical Quantum Computation

## Abstract

I describe the use of techniques based on composite rotations to combat systematic errors in controlled phase gates, which form the basis of two qubit quantum logic gates. Although developed and described within the context of Nuclear Magnetic Resonanace (NMR) quantum computing these sequences should be applicable to any implementation of quantum computation based on Ising couplings. In combination with existing single qubit gates this provides a universal set of robust quantum logic gates.

Quantum computers bennett00  use quantum mechanical effects to implement algorithms which are not accessible to classical computers, and thus are able to tackle otherwise intractable problems shor99 . They are extremely vulnerable to the effects of errors, and considerable effort has been expended on correcting random errors arising from decoherence processes shor95 ; steane96 ; steane99 . It is, however, also important to consider the effects of systematic errors, which arise from reproducible imperfections in the experimental apparatus used to implement quantum computations.

It makes sense to consider systematic errors as some of them can be tackled relatively easily. In the context of Nuclear Magnetic Resonance (NMR) quantum computation cory96 ; cory97 ; gershenfeld97 ; jones98c ; jones01a ; vandersypen01  the concept of composite rotations (also called composite pulses) ernst87 ; freeman97b ; levitt86 ; wimperis94  has been used to tackle both off-resonance effects cummins00  and more recently pulse length errors cummins02 . All these results, however, have been concerned with single qubit gates, and while these play a central role in quantum logic circuits it is also vital to consider two qubit controlled logic gates, such as the controlled-not gate.

The combination of a complete set of single qubit gates and the controlled-not gate is *universal* for quantum computation barenco95 , that is any quantum logic circuit can be built entirely from single qubit gates and controlled-not gates. There are, however, other two qubit gates which can be used to form a universal set, most notably the controlled phase-shift (controlled-$`\sigma _z`$) gate jones01a . This performs the transformation

$$|1|1|1|1$$
(1)
while leaving other eigenstates unchanged, and so is described by the matrix

$$\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right).$$
(2)
Applying Hadamard gates to the target spin before and after the controlled-$`\sigma _z`$ gate converts it to a controlled-$`\sigma _x`$ gate, that is controlled-not jones01a .

The controlled phase-shift gate can itself be considered as a derivative of the Ising ($`zz`$) coupling gate

$$\left(\begin{array}{cccc}e^{i\varphi /2}& 0& 0& 0\\ 0& e^{i\varphi /2}& 0& 0\\ 0& 0& e^{i\varphi /2}& 0\\ 0& 0& 0& e^{i\varphi /2}\end{array}\right)$$
(3)
which corresponds to evolution under the Ising ($`\sigma _{1z}\sigma _{2z}`$) Hamiltonian; Eqns. 2 and 3 are related by single qubit $`z`$ rotations and a global phase shift jones01a . The $`z`$ rotations can be adsorbed into abstract reference frames knill00  and the global phase shift has no physical significance, and thus the two equations are essentially equivalent.

The Ising coupling is of great experimental importance, as it provides the basic quantum logic gate for many proposed implementations of quantum computing cory96 ; lloyd93 ; ioffe99 ; cirac00 ; briegel01 ; raussendorf01 . Its role is particularly clear in NMR systems, where the scalar coupling ($`J`$ coupling) takes the Ising form in the weak coupling limit ernst87 , and from this point I will adopt the product operator notation widely used in NMR ernst87 ; jones01a ; sorensen83 . However the underlying principles remain applicable in any system based on Ising couplings.

Consider a system of two spin-$`\frac{1}{2}`$ nuclei, $`I`$ and $`S`$. The Ising gate (Eq. 3) is implemented by free evolution under the $`J`$ coupling Hamiltonian

$$_{IS}=\pi J\mathrm{\hspace{0.17em}2}I_zS_z$$
(4)
for a time $`\tau =\varphi /\pi J`$, where $`J`$ is the coupling strength and $`\varphi `$ is the desired evolution angle. In order to implement accurate controlled phase-shift gates by this means it is clearly necessary to know $`J`$ with corresponding accuracy.

Obtaining accurate values of $`J`$ for individual spin–spin couplings is easy in small molecules, but is much more difficult in larger systems. In particular, consider an array of qubits coupled by Ising interactions, with $`J`$ couplings that are nominally identical but in fact differ from one another as a result of imperfections in the lattice. In such a system it is desirable to be able to perform some accurately known $`zz`$ evolution over a *range* of values of $`J`$. While this might appear tricky, it is in fact easy to achieve using composite pulse techniques.

The problem of performing accurate $`zz`$ rotations is conceptually similar to that of correcting for pulse length errors in single qubit gates cummins02 , and the solutions are closely related. I begin by describing existing techniques for tackling pulse length errors. While a wide variety of composite pulse sequences have been developed within the context of NMR experiments levitt86 , most of these are not suitable for use in quantum computers as they make assumptions about the form of the initial state. Instead it is necessary to use Class A composite pulses levitt86  which work well for any initial state. Two suitable families of composite pulses have been described: the scrofulous pulse sequences of Cummins *et al.* cummins02 , and the BB1 family, initially described by Wimperis wimperis94 , and subsequently developed by Cummins *et al.* cummins02 . As the performance of the BB1 sequence is much better than that of scrofulous I will concentrate on the time symmetrized variant of BB1.

Consider a $`\theta _0`$ pulse, which implements a single qubit rotation around the $`x`$ axis on the Bloch sphere ernst87 . In an ideal system the propagator will be

$$U=\mathrm{exp}[i\theta I_x]$$
(5)
but in the presence of a fractional error $`ϵ`$ in the strength of the driving field the actual propagator will be

$$V=\mathrm{exp}[i\theta (1+ϵ)I_x].$$
(6)
The accuracy of this experimental implementation can be assessed using the propagator fidelity

$$=\frac{|\text{Tr}(VU^{})|}{\text{Tr}(UU^{})}$$
(7)
(note that it is necessary to take the absolute value of the numerator as $`U`$ and $`V`$ could in principle differ by an irrelevant global phase shift). In this case

$$=\mathrm{cos}\left(\frac{ϵ\theta }{2}\right)1\frac{ϵ^2\theta ^2}{8}.$$
(8)
A better approach is to replace the naive implementation, Eq. 6, with the composite pulse sequence

$$(\theta /2)_0\pi _{\varphi _1}\mathrm{\hspace{0.17em}2}\pi _{\varphi _2}\pi _{\varphi _1}(\theta /2)_0$$
(9)
with $`\varphi _2=3\varphi _1`$ and $`\varphi _1=\mathrm{arccos}(\theta /4\pi )`$. This gives a fidelity expression

$$1\frac{ϵ^6(\theta ^614\pi ^2\theta ^432\pi ^4\theta ^2)}{9216}$$
(10)
in which the second *and* fourth order error terms have been completely removed cummins02 .

A very similar approach can be used to tackle systematic errors in Ising coupling gates; in effect Ising coupling corresponds to rotation about the $`2I_zS_z`$ axis, and errors in $`J`$ correspond to errors in a rotation angle about this axis. These can be parameterised by the fractional error in $`J`$:

$$ϵ=\frac{J_{\text{real}}}{J_{\text{nominal}}}1.$$
(11)
Such errors can be overcome by rotating about a sequence of axes tilted from $`2I_zS_z`$ towards another axis, such as $`2I_zS_x`$. Defining

$$\theta _\varphi \mathrm{exp}[i\times \theta \times (2I_zS_z\mathrm{cos}\varphi +2I_zS_x\mathrm{sin}\varphi )]$$
(12)
allows the naive sequence $`\theta _0`$ to be replaced by Eq. 9 as before. The tilted evolutions can be achieved by sandwiching a $`2I_zS_z`$ rotation (that is, free evolution under the Ising Hamiltonian) between $`\varphi _y`$ pulses applied to spin $`S`$ ernst87 . For the case that $`\theta =\pi /2`$ (which forms the basis of the controlled-not gate) the final sequence takes the form shown in Fig. 1.

Clearly it is vital that any robust implementation of a quantum gate must be built from components that can themselves be implemented robustly. The robust Ising gate uses only two components: single qubit rotations around the $`\pm y`$ axes, for which robust versions are already known cummins02 , and periods of evolution under the Ising coupling. Note that the five time periods must have lengths in the integer ratios $`1:4:8:4:1`$, but it is not necessary to accurately control the absolute length of the periods, as errors in absolute lengths are equivalent to errors in the value of $`J`$.

The fidelity gain achieved for coupling gates by this approach is *identical* to that achieved for single qubit rotations. In particular for the case that $`\theta =\pi /2`$ the naive fidelity is

$$1\frac{ϵ^2\pi ^2}{32}$$
(13)
while the BB1 approach gives

$$1\frac{ϵ^663\pi ^6}{65536}.$$
(14)
Fidelity plots are shown in Fig. 2.

Clearly the BB1 approach compensates extremely well for small errors in $`J`$ values, especially within the range $`\pm 10\%`$. Over this range the infidelity of the BB1 sequence is always less that one part in $`10^6`$. To achieve comparable fidelity with a naive gate it is necessary to determine the value of $`J`$ to better than $`0.2\%`$. Thus the robust gate can achieve an infidelity of $`10^6`$ over a range of $`J`$ values more than 50 times wider than the naive gate. If even higher fidelities are desired the improvement provided by the robust gate is even greater. In combination with existing robust single qubit gates, the robust Ising gate provides a complete set of robust gates for quantum computation within the Ising model.

###### Acknowledgements.

I thank the Royal Society of London for financial support. I am grateful to S. Benjamin, L. Hardy, M. Mosca and A. M. Steane for helpful discussions.


text:
# Using Approximating Polynomials in Partial-Global Dynamical Simulations

## 1 MOTIVATION

Dynamical fermion simulations with smeared links are used more and more frequently. Smeared links improve the chiral and flavor symmetry of Wilson and staggered fermionic actions (see references in ). The main obstacle in using these actions is finding efficient algorithms to simulate them. In this paper we present some techniques that are used to improve the efficiency of these simulations. We will focus here on the HYP action , but we believe that these methods are more general and can be used for different actions as well.

The HYP action  is a modified staggered fermion action

$$S=S_G\left(U\right)+\overline{S}_G(V)+S_F\left(V\right)$$
(1)
where $`S_G`$ is the traditional pure gauge action defined in terms of thin links $`U`$, $`\overline{S}_G`$ is a gauge action defined in terms of smeared links $`V`$ and

$$S_F=\frac{n_f}{4}\mathrm{ln}det\mathrm{\Omega },$$
(2)
where $`\mathrm{\Omega }=M^{}M`$ restricted to even sites, is the staggered fermion action defined in terms of HYP links $`V`$ rather than thin links. The role of $`\overline{S}_G`$ is discussed in .

To simulate this action we use a two step approach. We first update a part of the thin links, $`UU^{}`$, using a heath-bath and/or an over-relaxation step, then we smear the thin links using the hyper-cubic (HYP) blocking  and compute the fat link action, $`S_G\left(V^{}\right)+S_F\left(V^{}\right)`$. We accept this new configuration with probability

$$P_{acc}=\mathrm{min}\{1,e^{\overline{S}_G\left(V^{}\right)+\overline{S}_G\left(V\right)S_F\left(V^{}\right)+S_F\left(V\right)}\}.$$
(3)
The most difficult step in the algorithm is computing the determinant ratio

$`e^{S_F\left(V^{}\right)+S_F\left(V\right)}=\left({\displaystyle \frac{det\mathrm{\Omega }\left(V^{}\right)}{det\mathrm{\Omega }\left(V\right)}}\right)^{n_f/4}`$
$`={\displaystyle \frac{𝑑\xi 𝑑\xi ^{}e^{\xi ^{}\left(\mathrm{\Omega }^{1/2}\mathrm{\Omega }\mathrm{\Omega }^{1/2}\right)^{n_f/4}\xi }}{𝑑\xi 𝑑\xi ^{}e^{\xi ^{}\xi }}}.`$ (4)
Instead of computing the ratio itself we use a stochastic estimator suggested by the last part of the relation above

$$E_A\left[\xi \right]=e^{\xi ^{}\left(A1\right)\xi },$$
(5)
where $`A=\left(\mathrm{\Omega }^{1/2}\mathrm{\Omega }\mathrm{\Omega }^{1/2}\right)^{n_f/4}`$. This estimator averages to the determinant ratio since $`A`$ is hermitian and positive definite. Furthermore, the variance of this estimator is

$$\sigma ^2\left(E_A\right)=\frac{1}{det\left(2A1\right)}\frac{1}{\left(detA\right)^2}.$$
(6)
The relation above holds only if the spectrum of $`A`$ is bounded from below by 1/2 otherwise the variance diverges. This condition is not necessarily fulfilled by our matrix. To address this problem we modify our stochastic estimator

$$E_A^{}\left[\xi \right]=e^{_{j=1}^n\xi _j^{}\left(A^{1/n}1\right)\xi _j}.$$
(7)
This estimator has a finite variance if the lowest eigenvalue of $`A`$ is greater than $`1/2^n`$. This condition can be fulfilled if we choose $`n`$ properly. However, this forces us to compute the $`n^{th}`$ root of fermionic matrices. This is why we are forced to use the polynomial approximation.

In order to make the algorithm more efficient we employed a reduced matrix : $`\mathrm{\Omega }_r=\mathrm{\Omega }e^{2f\left(\mathrm{\Omega }\right)}`$, where $`f`$ is a cubic function. The fermionic part will be $`S_F=\frac{n_f}{4}\mathrm{ln}det\mathrm{\Omega }_r`$ and the reminder is absorbed in $`\overline{S}_G`$.

Lastly, we mention that it can be proved that the algorithm satisfies the detailed balance condition for either of the estimators presented above.

## 2 POLYNOMIAL APPROXIMATION

To compute the estimator (7), where $`A=\left(\mathrm{\Omega }_r^{1/2}\mathrm{\Omega }_r\mathrm{\Omega }_r^{1/2}\right)^{n_f/4}`$, we need to find polynomial approximations for the matrix functions involved in the estimator

$`P_k^{\left(n\right)}\left(\mathrm{\Omega }\right)`$ $``$ $`F\left(\mathrm{\Omega }\right)^{1/2n},`$
$`Q_k^{\left(n\right)}\left(\mathrm{\Omega }\right)`$ $``$ $`F\left(\mathrm{\Omega }\right)^{1/n},`$ (8)
where $`F\left(\mathrm{\Omega }\right)=\left(\mathrm{\Omega }\mathrm{exp}\left[2f\left(\mathrm{\Omega }\right)\right]\right)^{n_f/4}`$, $`n`$ is the break-up level and $`k`$ is the polynomial order.

For a polynomial to approximate a matrix function it needs to approximate the function well on the entire spectrum of the matrix. To determine these polynomials we employed a least square method . The coefficients $`t_i`$ of the polynomial $`T`$ approximating the function $`g`$ are determined by minimizing the quadratic “distance”

$$\delta ^2\left(t_i\right)=_{\lambda _0}^{\lambda _1}𝑑x\rho (x)\left(T\left(x\right)g\left(x\right)\right)^2,$$
(9)
where $`\lambda _0`$, $`\lambda _1`$ are the spectral bounds of $`\mathrm{\Omega }`$ and $`\rho `$ is a weight function, ideally $`\mathrm{\Omega }`$’s spectral density.

This method has the advantage that it turns into a linear problem and the solution is

$$t_i=\underset{j=0}{\overset{k}{}}\left(K^1\right)_{ij}g_j,$$
(10)
where $`k`$ is the polynomial order and

$`K_{ij}`$ $`=`$ $`{\displaystyle _{\lambda _0}^{\lambda _1}}𝑑x\rho \left(x\right)x^{i+j},`$
$`g_i`$ $`=`$ $`{\displaystyle _{\lambda _0}^{\lambda _1}}𝑑x\rho \left(x\right)x^ig\left(x\right).`$ (11)
For our polynomials we used $`\lambda _0=4m^2`$, the lowest eigenvalue for the staggered fermions matrix. We set the upper bound to $`\lambda _1=16.4`$ which we found to be a good limit for the range of parameters we used in our simulations. It should be mentioned that the only hard upper bound for the spectrum of $`\mathrm{\Omega }`$ is $`64+4m^2`$ but this limit is reduced dynamically. The smearing reduces the value of the highest eigenvalue even further.

The weight function $`\rho \left(x\right)`$ should approximate the spectral density for $`\mathrm{\Omega }`$. However, we found that for reasonable weight functions the coefficients depend very little on the detailed form. Thus, we used a weight function of the form $`\rho \left(x\right)=x^\alpha `$ which is more convenient for our computation. The advantage of this form is that we can compute the $`K_{ij}`$ matrix analytically and also it generates a recursion relation for the coefficients $`g_i`$.

To measure the accuracy of the approximation we used the distance function (9) where we set $`\rho \left(x\right)=1`$. In Fig.1 we plot the distance between $`\left[P^{\left(n\right)}\right]^{2n}`$, $`\left[Q^{\left(n\right)}\right]^n`$ and $`F\left(x\right)^1`$ as a function of the polynomial order. We notice first that $`Q^{\left(1\right)}`$ is a much better approximation than $`Q^{\left(n\right)}`$ for $`n>1`$. This is due to the fact that $`F\left(x\right)`$ (which is the function that $`Q^{\left(1\right)}`$ approximates) has an almost polynomial behavior close to $`\lambda _0`$ where the function is the hardest to approximate. Another interesting feature that we notice in the graphs is that the order of the approximation doesn’t improve with increasing break-up level. This is puzzling, at first, since as $`n`$ increases the equivalent order of the polynomial approximating $`F\left(x\right)^{\pm 1}`$ is $`n\times k`$. However, the number of variables varied to minimize $`\delta `$ stays the same and it seems that this is the determining factor for the order of approximation.

We found that when $`\delta 10^8`$ the numerical roundoff error becomes dominant and the higher order polynomials do not improve the precision of the calculation. Thus, we will be using the polynomials with $`\delta 10^8`$ as “exact” ones. We see in Fig. 1 that these polynomials have an order of $`k100`$. The only exceptions are the polynomials for $`n=1`$ but they are of no real interest since they don’t satisfy the minimum eigenvalue requirement.

To further reduce the cost of computing the estimator we implemented a two step approach. Since we only use the value of the estimator in the accept reject step we only need to know whether the value of the estimator is larger or smaller that $`\mathrm{ln}r`$, where $`r`$ is the random number. We will first compute the estimator approximatively, using small order polynomials; if the distance between the approximate value and $`\mathrm{ln}r`$ is greater than $`ϵ`$ we use this value in the accept reject step. If not, we “fall-back”, i.e. we compute the “exact” value of the estimator using the large polynomials. The value of $`ϵ`$ is determined in a small run as the maximum difference between the approximate estimator and its “exact” value. The rate of fall-back is usually around 10% and this reduces the computational cost substantially. Having fixed the large polynomials by the $`\delta `$ condition we determine the optimum small ones that minimize the computational cost. The cost is given by the number of $`M^{}M`$ multiplies and it is determined by the order of the polynomials used. We found that the order of the optimum small polynomial is independent of the break-up level (see Fig.2) and it increases as we decrease the mass. We find that for $`am=0.10`$ we need $`k=24`$, for $`am=0.06`$ $`k=32`$ and for $`am=0.04`$ $`k=48`$.

## 3 CONCLUSIONS

We have showed how to use the polynomial approximation for simulating dynamical fermions. First we determine the minimum break-up level $`n`$ using the minimum eigenvalue requirement. Then we determine the large polynomials using the $`\delta 10^8`$ conditions. We then determine the optimum small polynomials using the minimum cost condition. For a more complete discussion see .


text:
# A 1'' Telescope: The Optimal Approach to Bright-Star Planetary Transits

## 1 Introduction

In Pepper, Gould & DePoy (astro-ph/0208042), we argued that all-sky surveys are the best way to find planets transiting bright $`V<10`$ stars. Such transits offer the best opporunity for detailed studies of planets. We showed that the number of systems probed is given by

$$N_p=\frac{4}{3}\pi n\eta d_0^3\left(\frac{R_0}{a_0}\right)\left(\frac{L}{L_0}\right)^{3/2}\left(\frac{R}{R_0}\right)^{7/2}\left(\frac{a}{a_0}\right)^{5/2}\left(\frac{r}{r_0}\right)^6,$$
(1)
where $`d_0`$ is the maximum distance at which an $`i=90^\mathrm{o}`$ transit can be detected for a star of luminosity $`L_0`$, radius $`R_0`$, with a planet at semi-major axis $`a_0`$ and radius $`r_0`$, and where $`n`$ is the local number density of such stars and $`\eta =0.719`$ is a numerical factor.

We then showed that the sensitivity of a given survey will depend almost entirely on $`\gamma `$, the number of photons collected from a fiducial star of some fixed designated magnitude, which we arbitrarily chose to be $`V=10`$. That is, it will not depend on the details of the all-sky observing program. We then normalized equation (1) in terms of $`\gamma `$,

$$N_p=600F(M_V)\left(\frac{a}{a_0}\right)^{5/2}\left(\frac{r}{r_0}\right)^6\left(\frac{\gamma }{\gamma _0}\right)^{3/2}\left(\frac{\mathrm{\Delta }\chi _{\mathrm{min}}^2}{36}\right)^{3/2}$$
(2)
where we adopted $`\gamma _0=2\times 10^7`$, $`a_0=20R_{}`$, $`r_0=0.10R_{}`$, and where we have made our evaluation at $`M_V=5`$ (i.e. $`R=0.97R_{}`$, $`V_{\mathrm{max}}=10`$, $`d_0=100`$pc, and $`n=0.004\mathrm{pc}^3`$). Here, $`F(M_V)`$ is a function, which is shown in Figure 1, and $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ is the minimum $`\chi ^2`$ improvement of a transit-model relative to a constant-flux model, which is set to avoid excessive false positives. Note that $`\gamma _0=2\times 10^7`$ corresponds approximately to 1000 20-second exposures with a 5 cm telescope and a broadened $`(V+R)`$ type filter for one $`V=10`$ mag fiducial star.

Our goal was to apply this formula to the problem of telescope design in a subsequent paper. However, in light of the referee report, we decided to append our work on telescope design as an additional chapter of the original paper. The following is that chapter.

## 2 Implications for Telescope Design

We now apply the general analysis of astro-ph/0208042 to the problem of optimizing telescope design for quickly locating a “large” number of bright transiting systems. Since only one such system is now known, we define “large” as $`𝒪(10)`$. From equation (2) and the 0.75% frequency of hot jupiters measured from RV surveys, there are roughly five such systems to be discovered over the whole sky per magnitude for $`V_{\mathrm{max}}=10`$. Hence, from Figure 1, of order 15 are to be discovered from all spectral types. It would, of course, be possible to discover even more by going fainter, but setting this relatively bright limit is advisable for three reasons. First, as we argued in the introduction, the brightest transits are the most interesting scientifically, and most of the transits detected in any survey will be close to the magnitude limit. Second, as we discuss below, a wide range, $`\mathrm{\Delta }V=V_{\mathrm{max}}V_{\mathrm{min}}`$, of survey sensitivity can only be achieved at considerable cost to the observing efficiency. Hence, if high efficiency is to be maintained, setting $`V_{\mathrm{max}}`$ fainter means eliminating the brightest (most interesting) systems. Third, at $`V_{\mathrm{max}}=10`$, we are already reaching distances of 100 pc for G stars. Hence the number of transits observed in fainter surveys will not continue to grow as $`d_0^3`$ as in equation (1).

In previous sections, we ignored the loss of sensitivity to systems that are brighter than $`V_{\mathrm{min}}`$, which is set by saturation of the detector (or more precisely, by the flux at which detector non-linearities can no longer be accurately calibrated). This fraction is $`10^{0.6\mathrm{\Delta }V}`$, or 6% for $`\mathrm{\Delta }V=2`$, which we therefore adopt as a sensible goal. That is, we wish to optimize the telescope design for,

$$8=V_{\mathrm{min}}<V<V_{\mathrm{max}}=10.$$
(3)
(In any event, essentially all stars $`V<8`$ have already been surveyed for XSPs using RV, and the problem of determining which among the planet-bearers have transits is trivial compared to the problem of conducting an all-sky photometric variability survey.)

Optimization means maximizing the photon collection rate, $`\gamma /T`$, where $`T`$ is the duration of the experiment and, again, $`\gamma `$ is the total number of photons collected from a fiducial $`V=10`$ mag star. Explicitly,

$$\gamma =KTD^2\frac{(\mathrm{\Delta }\theta )^2}{4\pi },$$
(4)
where $`\mathrm{\Delta }\theta `$ is the angular size of the detector, $`D`$ is the diameter of the primary-optic, $``$ is the fraction of the time actually spent exposing, and $`K`$ is a constant that depends on the telescope, filter, and detector throughput. For our calculations, we assume $`K=K_040e^{}\mathrm{cm}^2\mathrm{s}`$, which is appropriate for a broad $`(V+R)`$ filter and the fiducial $`V=10`$ mag star. The design problems are brought into sharper relief by noting that $`\mathrm{\Delta }\theta =/D`$, where $``$ is the linear size of the detector and $``$ is the focal ratio, or $`f`$/#, of the optics. Equation (4) then becomes

$$\gamma =\frac{K^2T}{4\pi ^2}.$$
(5)
That is, almost regardless of other characteristics of the system, the camera should be made as fast as possible. We will adopt $`=2`$, beyond which it becomes substantially more difficult to design the optics. A more remarkable feature of equation (5) is that all explicit dependence on the size of the primary optic has vanished: a 1” telescope and an 8m telescope would appear equally good! Actually, as we now show, there is a hidden dependence of $``$ on $`D`$, which favors small telescopes.

The global efficiency $``$ can be broken down into two factors, $`=_0_S`$, where $`_0`$ is the fraction of time available for observing (i.e., during which the sky is dark, the weather is good, etc.), and $`_S`$ is the fraction of this available observing time that the shutter is actually open. The first factor is not affected by telescope design and so will be ignored for the moment. The second factor should be maximized. Since the readout time is fixed, the smaller is the telescope aperture, the longer can be the exposures before a $`V_{\mathrm{min}}=8`$ mag star saturates, and so the smaller fraction of time is lost to read-out. For large-format detectors, the pixel size is typically $`\mathrm{\Delta }x_p=15\mu `$m, for which significant non-linearities set in at about $`60,000e^{}`$. Hence, exposure times are

$$T_{\mathrm{exp}}=38\mathrm{s}\frac{K_0}{K}\left(\frac{D}{2.5\mathrm{cm}}\right)^2$$
(6)
That is, the exposure time for a $`1^{\prime \prime }`$ “telescope” is already of the order of typical read-out times for large-format detectors. Clearly, smaller is better, but are there constraints from going too small?

One potential constraint comes from sky noise. To stay within the photon-noise limited regime, which has been assumed in all of our calculations, the sky inside one pixel should be at least one magnitude fainter than $`V_{\mathrm{max}}`$, and preferably two mags fainter. Assuming the sky brightness in our broadening passband reaches a maximum of $`V=19.7\mathrm{arcsec}^2`$, and again assuming $`\mathrm{\Delta }x_p=15\mu `$m pixels, the sky in one pixel is

$$V_{\mathrm{sky}}=10.8+5\mathrm{log}\left(\frac{D}{2.5\mathrm{cm}}\frac{}{2}\right).$$
(7)
Hence, the sky is a bit bright for a $`1^{\prime \prime }`$ telescope, but appears quite satisfactory for a $`2^{\prime \prime }`$.

Finally, one must be careful that the field of view is not too large, or the focal-plane distortions at its edges will be difficult (and expensive) to correct. For example, for a $`4k\times 4k`$ detector with $`\mathrm{\Delta }x_p=15\mu `$m pixels, $`\mathrm{\Delta }\theta =68^{}(D/2.5\mathrm{cm})^1(/2)^1`$. At $`1^{\prime \prime }`$, this is probably too large to correct at reasonable expense. With this size detector, both sky-noise considerations, and problems in optics design argue for a $`1.5^{\prime \prime }`$ telescope. However, from equation (6), such a “large” telescope will most likely be dominated by read-out time. In summary, we conclude that a $`1^{\prime \prime }`$ to $`2^{\prime \prime }`$ telescope equipped with a $`4x\times 4k`$, $`=6`$cm detector and with $`=2`$ focal ratio is optimal for this observing problem. Among all existing transit programs of which we are aware, the WASP telescope ($`D=2.5^{\prime \prime }`$ lens, $`=2.8`$ focal ratio, $`2k\times 2k`$, $`=3`$cm detector, Street et al. 2002) comes closest to meeting these design specifications.

Adopting $`\gamma =\gamma _0`$, $`K=K_0`$, $`_0=20\%`$, $`_S=50\%`$, $`=6`$cm, and $`=2`$, we find from equation (5) that the required duration of the experiment using our optimally-designed telescope is $`T=3\mathrm{months}`$. Clearly, roughly a year is required just to get around the sky. The S/N acquired during such a year-long search would therefore roughly double the minimum requirements calculated here.

This work was supported in part by grant AST 02-01266 from the NSF.


text:
# Pareto’s law: a model of human sharing and creativity

\[

## Abstract

A computational model for the distribution of wealth among the members of an ideal society is presented. It is determined that a realistic distribution of wealth depends upon two mechanisms: an asymmetric flux of wealth in trading transactions that advantages the poorer of the two traders and a non-stationary creation and destruction of individual wealth. The former mechanism redistributes wealth by reducing the gap between the rich and poor, leading to the emergence of a middle class. The latter mechanism, together with the former one, generates a distribution of wealth having a power-law tail that is compatible with Pareto’s law.

\]

More than a century ago the Italian sociologist and economist, Vilfredo Pareto, studied the distribution of income among people of different western countries and found an inverse power law. The cumulative probability $`P(w)`$ of people whose income is at least $`w`$ seemed to vary as $`P(w)w^\alpha `$. Pareto mistakenly believed that the exponent $`\alpha `$ was universal constant with an approximate value of $`\alpha =1.5`$. In our discussion we use the terms distribution of income and distribution of wealth interchangeably, and although the two distributions may not exactly coincide, they are in fact strongly dependent on one another. The distinction between the two will not influence the conclusions drawn from the dynamical models of economic interactions we present here.

Societies in the West have historically been partitioned into three classes; the poor, the middle and the rich. The relative size of each class is determined by somewhat arbitrarily assigned levels of income, but it is safe to say that the smallest class is the rich and for a stable society the largest is the middle class. This partitioning must be incorporated into any mathematical model describing how wealth is distributed within a society. The inverse power law of Pareto does not have this characteristic partitioning, because the derivative of a power law is still a power law. According to data, across the full range of income, we should expect the probability density function (pdf) $`p(w)`$ to increase at low income, reach a maximum and, finally, decrease with increasing wealth. Moreover, fitting the cumulative probability $`P(w)`$ instead of the pdf $`p(w)`$ will certainly mask important properties of the dynamics of the economy. Therefore, herein we model the density function $`p(w)`$. The cumulative probability $`P(w)`$ is related to the pdf $`p(w)`$ via the relation $`P(w)=_w^{\mathrm{}}p(x)𝑑x`$.

We investigate a computational model of the distribution of wealth, which is based on the assumption that the pdf of wealth is the result of a non-stationary total wealth and asymmetric trading mechanisms in the economy. For a non-stationary economy we postulate a situation in which the total wealth of a society is not conserved, since an individual may create or consume wealth. For an asymmetric trading economy we postulate a mechanism by which, through trades, wealth randomly moves from one individual to another, but in such a way as to, on average, slightly advantage the poorer trader in the exchange.

Wealth decays because of human consumption and needs to be continuously recreated through work and human creativity. People must work to live! The mechanisms that regulate the creation and consumption of wealth are neither identical nor symmetric, and they do not compensate exactly for one another. This lack of compensation implies that the wealth of an individual, even in the absence of any type of trade, fluctuates in time. Individual wealth may increase or decrease according to the situation. Moreover, also the assumption that all transactions occur at precisely the right price, a fact that implies the absence of any flux of wealth between traders, is unrealistic. Some sellers are more persuasive and some buyers are more discerning than others. There are good and bad deals. Therefore, the proper value of an asset may be lower or higher than the price paid for it. This dynamic of a transaction continuously generates a random flux of wealth from one trader to the other in every transaction.

Finally, we assume that this random flux of wealth during a transaction depends upon the difference in wealth between the two traders and may be statistically biased toward the poorer partner. This mechanism seems to be needed to redistribute the wealth among people. More importantly, without this mechanism, it would be true that the rich would get richer and the poor would get poorer. But this latter flux of wealth from the poor to the rich would be a consequence of the excess wealth of the rich and not the result of any hypothetical economic abuses of the poor by the rich. In a fair game in which the rich and the poor have the same probability to win or lose in any transaction involving one of each class, the rich ultimately win because of their larger resources . So for a society to avoid the rich getting too rich and the poor getting too poor with the fallout of a social collapse, the poor must be slightly advantaged in their transactions with the rich.

The symmetry-breaking mechanism that advantages the poor, in their transactions with the rich, should be thought to be both a necessity for, and a consequence of, an efficient and stable society. In fact, the necessity to slightly advantage the poor requires a policy of rights and freedom, balanced by duties and responsibilities, that has been typical of western societies. For example, in a modern democratic society the right of the employed worker to be protected by a union is recognized, as is their right to strike against their employers in particular cases. These rights allow the workers to obtain a salary statistically higher than the real value of their job, compatible to the actual economy. In the same way, the freedom to decide the price of their own goods is not only an incentive to a larger and better production of goods, but allows the self-employed of the working class to statistically increase the price of their own products. In fact, many such workers are farmers, fishermen and small scale artisans, who produce the basic needs of a society, for example, food. Even the richest person needs to eat, therefore, if he/she cannot enslave the farmer and wants to eat, he/she must buy food at the price the farmer decides. We understand that the economic situation is much more complex than we are able to detail here, but the policy of rights and freedom in an economic context contributes the symmetry-breaking mechanism in our model as we describe.

In the physics tradition we hypothesize a mathematical model describing the economic mechanisms in an idealized situation, predict the resulting distribution of wealth and then compare the predicted distribution function with data . The correctness of our assumptions are ultimately determined by data and not logical discourse. In all our simulations we suppose an economic society of 10,000 agents, who at the onset have the same degree of wealth. We calculate the wealth pdf, $`p(w)`$, of our hypothetical society after 100 million random transactions between all agents. We proceed by steps, introducing models of increasing complexity and discussing the consequences of each of them in turn. The first three models assume the conservation of total wealth, while the last model assumes that the total wealth of a society may change.

Perfect trade price. The simplest hypothesis is to suppose that wealth cannot be created or destroyed and that all transactions between agents occur at the right price. This complete symmetry means that each agent acquires the exact same amount of wealth that is given up to the other trader. The model is based upon the assumption that the price of an asset is defined by its intrinsic value. Consequently the wealth of each agent remains constant and the final distribution of a closed society’s wealth remains identical to the initial uniform distribution. This model is clearly at odds with our experience of the economic world.

Random symmetric trade price. Again we suppose that wealth cannot be created or destroyed, however we now assume that the transactions between agents may occur at a price that randomly fluctuates around a hypothetical ideal value. Therefore, the wealth of each agent may increase or decrease according to whether he/she makes good or bad deals. We still assume symmetry, that is, both agents have the same chance to win or lose, such as in a coin toss. We also suppose that the maximum wealth that may move from one agent to the other, in a single transaction, is limited to some fraction of the wealth of the poorer of the two traders. This restrictive condition is required because only in a robbery is it realistic to assume that a person will give over to another more wealth than he/she possesses.

Fig. 1 shows the pdf of wealth among the agents after 100 million random transactions. It is a very wide inverse power law of the type $`p(w)1/w`$. The fraction of wealth that moves from one agent to the other in a single transaction is assumed to vary between 0% and 50% of the poorer agent’s wealth. Fig. 1 shows that the distance in wealth between the richest and the poorest is huge. Practically, the entire wealth of society concentrates in the hand of very few people. In fact, even if both agents in a transaction have the same chance to win or lose, the risk for the rich trader is smaller because if he/she loses, the loss is a smaller fraction of his/her own wealth than that which the poorer agent may lose. Consequently, there is a high probability in this model that a very few people accumulate almost the entire wealth available and the others become extremely poor, as Fig. 1 shows. An ideal society that adopts such a trading policy will ultimately collapse.

Random asymmetric trade price. The third model assumes that the poorer of the two traders is slightly advantaged in the trade, and has a greater chance to formulate a good deal. In practice, we implement this ability to make a good deal by assuming that if $`w_p`$ and $`w_r`$ are the wealth of the poorer and the richer of the two traders, respectively, the probability for the poorer trader to profit is

$$\mathrm{\Pi }=0.5+f\frac{w_rw_p}{w_r+w_p},$$
(1)
where $`0f0.5`$ is the asymmetry flux index. Eq. (1) assumes that if the two agents have the same degree of wealth, they have the same chance to get a good deal, while if one trader is much richer than the other, the poorer is advantaged with the maximum probability of $`\mathrm{\Pi }=0.5+f`$.

In the computer simulations shown in Fig. 2, we calculate three different situations with $`f=0.10`$, $`f=0.15`$ and $`f=0.25`$, and the fraction of wealth that may move in the transaction is assumed to vary between 0% and 10% of the poorer agent’s wealth. We determine that these computer-generated wealth pdf, as Fig. 2 shows, are very well fit by curves of the form

$$p(w)=aw^\gamma \mathrm{exp}(bw^\delta ).$$
(2)
Fig. 2 shows that with this model the economy of a society does not collapse as it did in the previous case. The random asymmetric trade price model yields a stable distribution of wealth, $`p(w)`$, well fitted by (2). The present model leads to a partitioning of society into three parts: poor, middle and rich classes, where the middle is the largest class followed by the poor and finally the rich. The economic gap between richest and poorest is not unrealistically wide, as it was in the random symmetric model. If the curves of Fig. 2 seem unrealistic, it is because the gap in the wealth between the richest and poorest is too small! However, the figure shows that the separation in wealth between the rich and poor increases by decreasing the asymmetric flux index $`f`$ and a more realistic distribution with an exponential-like tail may emerge by decreasing this parameter.

In summary, the random asymmetric trade price model produces a stable wealth pdf that may be fitted by functions of the type of Eq. (2). The only difficulty with this model is that it does not reproduce many experimental data, like those analyzed by Pareto, which suggest an inverse power-law distribution that is not fitted by Eq. (2). We shall show that this failure is due to the assumption that society’s total wealth remains constant.

Random asymmetric trade price with nonstationary wealth. Our last computational model implements the mechanism of the previous model with the additional property that the total wealth is not constant. The wealth of each individual trader is subject to a random destruction-creation dynamics, due to the natural ebb and flow of wealth through human creativity, work and investment. We introduce the non-stationary mechanism through a risk index $`r>0`$, that measures the standard deviation of the destruction-creation wealth process. In our computer simulation we assume that after every 10,000 transactions the wealth of each trader is reinitialized by the following expression

$$w_i(t+1)=(1+r\xi )w_i(t),$$
(3)
where $`w_i(t)`$ is the wealth of the i-$`th`$ agent after the t-$`th`$ epoch and $`w_i(t+1)`$ is that agent’s wealth at the start of the new epoch, $`t+1`$. Finally, $`\xi `$ is a zero-centered Gaussian random variable with unit variance. The rare instances of too negative a random value of the variable $`\xi `$, that makes the wealth $`w(i)`$ negative, are neglected.

Fig. 3a shows the distribution of wealth for three computer simulations obtained by assuming a fixed risk index $`r=0.1`$ and varying the asymmetry flux index $`f`$: $`f=0.10`$, $`f=0.15`$ and $`f=0.30`$. Fig. 3b shows another three computer simulations obtained by assuming a fixed asymmetry flux index $`f=0.1`$ and varying the risk index, $`r=0.05`$, $`r=0.10`$ and $`r=0.15`$. The fraction of wealth that may move between agents in a single transaction is assumed again to vary between 0% and 10% of the poorer agent’s wealth. Finally, Figs. 3a and 3b show that the computer-generated wealth pdf can be well fitted for a very large region by curves of the form

$$p(w)=aw^\gamma /(1+bw)^{\gamma +\delta }.$$
(4)
As in the previous model, we obtain a stable wealth distribution $`p(w)`$ that partitions society into the three classes. The wealth gap between the richest and poorest increases by decreasing the asymmetric flux index $`f`$ and/or by increasing the risk index $`r`$.

We stress that Eq. (4) is characterized by an inverse power-law tail of the type $`p(w)1/w^\delta `$, where $`\delta `$ is the inverse power-law exponent. This inverse power-law distribution is compatible with Pareto’s measurements. In fact, Fig. 4 shows the cumulative probability $`P(w)`$ for an economy with $`f=0.15`$ and $`r=0.075`$, fitted with the Pareto distribution and having a Pareto exponent of $`\alpha =\delta 1=1.48`$. This figure is consistent with the fit to U.S. income data made in Figure 35 of Ref..

A society that allows a small class of people to absorb all its wealth is unstable and will either collapse economically, or be destroyed by revolution. This may also happen in the absence of economic exploitation of the poor by the rich, as we explained above. A stable society requires that the poor have an advantage in transactions with the wealthy and are protected by particular rights and marketing freedom. On the other hand, in a real human society, wealth does not only move randomly from one individual to another through trades, but is continuously created and destroyed through an individual’s work, creativity and investment. The human aspiration of personal riches requires a positive risk index $`r`$ that generates both a large middle class and a small rich class, when applied across the society. The asymmetry flux index $`f`$ and the risk index $`r`$ express the human need to share and the human need to create, respectively. Therefore, we conclude from the model that both mechanisms must be present to some degree in any human society that functions effectively and remains stable. In fact, when these two mechanisms are present, according to our model, the wealth pdf assumes the inverse power-law distribution of Pareto at the high income end, and, more realistically, still retains a small but finite population at the low income end.

Acknowledgment:  N.S. thanks the Army Research Office for support under grant DAAG5598D0002.