text:
# Diophantine Methods for Exponential Sums, and Exponential Sums for Diophantine Problems

## 1. Introduction

Over the past fifteen years or so, the Hardy-Littlewood method has experienced a renaissance that has left virtually no facet untouched in its application to diophantine problems. Our purpose in this paper is to sketch what might be termed the past, present, and future of these developments, concentrating on aspects associated with diagonal diophantine problems, and stressing modern developments that make increasing use of less traditional diophantine input within ambient analytic methods. We avoid discussion of the Kloosterman method and its important recent variants (see  and ), because the underlying ideas seem inherently constrained to quadratic, and occasionally cubic, diophantine problems. Our account begins with a brief introduction to the Hardy-Littlewood (circle) method, using Waring’s problem as the basic example. The discussion here illustrates well the issues involved in the analysis of systems of diagonal equations over arbitrary algebraic extensions of $``$, and motivates that associated with more general systems of homogeneous equations (see  and ).

Let $`s`$ and $`k`$ be natural numbers with $`s>k2`$, and consider an integer $`n`$ sufficiently large in terms of $`s`$ and $`k`$. The circle method employs Fourier analysis in order to obtain asymptotic information concerning the number, $`R(n)=R_{s,k}(n)`$, of integral solutions of the equation $`x_1^k+\mathrm{}+x_s^k=n`$. Write $`P=n^{1/k}`$ and define the exponential sum $`f(\alpha )=f(\alpha ;P)`$ by

$$f(\alpha )=\underset{1xP}{}e(\alpha x^k),$$

wherein $`e(z)`$ denotes $`e^{2\pi iz}`$. Then it follows from orthogonality that

$$R(n)=_0^1f(\alpha )^se(n\alpha )𝑑\alpha .$$

When $`\alpha `$ is well-approximated by rational numbers with small denominators, one has sharp asymptotic information concerning $`f(\alpha )`$. In order to be precise, let $`Q`$ satisfy $`1Q\frac{1}{2}P^{k/2}`$, and define the major arcs $`𝔐=𝔐(Q)`$ to be the union of the intervals $`𝔐(q,a)=\{\alpha [0,1):|q\alpha a|QP^k\}`$, with $`0aqQ`$ and $`(a,q)=1`$. Also, put

$$S(q,a)=\underset{r=1}{\overset{q}{}}e(ar^k/q)\text{and}v(\beta )=_0^Pe(\beta \gamma ^k)𝑑\gamma ,$$

and define $`f^{}(\alpha )`$ for $`\alpha [0,1)`$ by taking $`f^{}(\alpha )=q^1S(q,a)v(\alpha a/q)`$, when $`\alpha `$ lies in $`𝔐(q,a)𝔐(Q)`$, and otherwise by setting $`f^{}(\alpha )=0`$. Then the sharpest available estimate (see Theorem 4.1 of ) establishes that 111Given a complex-valued function $`f(t)`$ and positive function $`g(t)`$, we use Vinogradov’s notation $`f(t)g(t)`$, or Landau’s notation $`f(t)=O(g(t))`$, to mean that when $`t`$ is large, there is a positive number $`C`$ for which $`f(t)Cg(t)`$. Similarly, we write $`f(t)g(t)`$ when $`g(t)f(t)`$, and $`f(t)g(t)`$ when $`f(t)g(t)f(t)`$. Also, we write $`f(t)=o(g(t))`$ when as $`t\mathrm{}`$, one has $`f(t)/g(t)0`$. Finally, we use the convention that whenever $`ϵ`$ occurs in a formula, then it is asserted that the statement holds for each fixed positive number $`ϵ`$. $`f(\alpha )=f^{}(\alpha )+O(Q^{1/2+ϵ})`$, uniformly for $`\alpha 𝔐(Q)`$. The functions $`S(q,a)`$ and $`v(\beta )`$ are rather well-understood, and thus one deduces that whenever $`s\mathrm{max}\{4,k+1\}`$ and $`QP`$, then

$$_𝔐f(\alpha )^se(n\alpha )𝑑\alpha =\frac{\mathrm{\Gamma }(1+1/k)^s}{\mathrm{\Gamma }(s/k)}𝔖_{s,k}(n)n^{s/k1}+O(n^{s/k1\delta }),$$
$`1.1`$
for a suitable positive number $`\delta `$. Here, the $`\mathrm{\Gamma }`$-function is that familiar from classical analysis, and the singular series $`𝔖_{s,k}(n)`$ is equal to the product of $`p`$-adic densities $`_pv_p(n)`$, where for each prime $`p`$ we write

$$v_p(n)=\underset{h\mathrm{}}{lim}p^{h(1s)}\text{card}\{𝕩(/p^h)^s:x_1^k+\mathrm{}+x_s^kn(modp^h)\}.$$

An asymptotic formula for $`R(n)`$, with leading term determined by the major arc contribution (1.1), now follows provided that the corresponding contribution arising from the minor arcs $`𝔪=[0,1)𝔐`$ is asymptotically smaller. Although such is conjectured to hold as soon as $`s\mathrm{max}\{4,k+1\}`$, this is currently known only for larger values of $`s`$. It is here that energy is focused in current research. One typically estimates the minor arc contribution via an inequality of the type

$$\left|_𝔪f(\alpha )^se(n\alpha )𝑑\alpha \right|\left(\underset{\alpha 𝔪}{sup}|f(\alpha )|\right)^{s2t}_0^1|f(\alpha )|^{2t}𝑑\alpha .$$
$`1.2`$
For suitable choices of $`t`$ and $`Q`$, one now seeks bounds of the shape

$$\underset{\alpha 𝔪}{sup}|f(\alpha )|P^{1\tau +ϵ}\text{and}_0^1|f(\alpha )|^{2t}𝑑\alpha P^{2tk+\delta +ϵ},$$
$`1.3`$
with $`\tau >0`$ and $`\delta `$ small enough that $`(s2t)\tau >\delta `$. The right hand side of (1.2) is then $`o(n^{s/k1})`$, which is smaller than the main term of (1.1) whenever $`𝔖_{s,k}(n)1`$. The latter is assured provided that non-singular $`p`$-adic solutions can be found for each prime $`p`$, and in any case when $`s4k`$. Classically, one has two apparently incompatible approaches toward establishing the estimates (1.3). On one side is the differencing approach introduced by Weyl , and pursued by Hua , that yields an asymptotic formula for $`R(n)`$ whenever $`s2^k+1`$. The ideas introduced by Vinogradov , meanwhile, provide the desired asymptotic formula when $`s>Ck^2\mathrm{log}k`$, for a suitable positive constant $`C`$.

## 2. Efficient differencing and smooth Weyl sums

Since the seminal work of Vaughan , progress on diagonal diophantine problems has been based, almost exclusively, on the use of smooth numbers, by which we mean integers free of large prime factors. In brief, one seeks serviceable substitutes for the estimates (1.3) with the underlying summands restricted to be smooth, the hope being that this restriction might lead to sharper bounds. Before describing the kind of conclusions now available, we must introduce some notation. Let $`𝒜(P,R)`$ denote the set of natural numbers not exceeding $`P`$, all of whose prime divisors are at most $`R`$, and define the associated exponential sum $`h(\alpha )=h(\alpha ;P,R)`$ by

$$h(\alpha ;P,R)=\underset{x𝒜(P,R)}{}e(\alpha x^k).$$

When $`t`$ is a positive integer, we consider the mean value $`S_t(P,R)=_0^1|h(\alpha )|^{2t}𝑑\alpha `$, which, by orthogonality, is equal to the number of solutions of the diophantine equation $`x_1^k+\mathrm{}+x_t^k=y_1^k+\mathrm{}+y_t^k`$, with $`x_i,y_i𝒜(P,R)`$ $`(1it)`$. We take $`RP^\eta `$ in the ensuing discussion, with $`\eta `$ a small positive number222We adopt the convention that whenever $`\eta `$ appears in a statement, implicitly or explicitly, then it is asserted that the statement holds whenever $`\eta >0`$ is sufficiently small in terms of $`ϵ`$.. In these circumstances one has $`\text{card}(𝒜(P,R))c(\eta )P`$, where the positive number $`c(\eta )`$ is given by the Dickman function, and it follows that $`S_t(P,R)P^t+P^{2tk}`$. It is conjectured that in fact $`S_t(P,R)P^ϵ(P^t+P^{2tk})`$. We refer to the exponent $`\lambda _t`$ as being permissible when, for each $`ϵ>0`$, there exists a positive number $`\eta =\eta (t,k,ϵ)`$ with the property that whenever $`RP^\eta `$, one has $`S_t(P,R)P^{\lambda _t+ϵ}`$. One expects that the exponent $`\lambda _t=\mathrm{max}\{t,2tk\}`$ should be permissible, and with this in mind we say that $`\delta _t`$ is an associated exponent when $`\lambda _t=t+\delta _t`$ is permissible, and that $`\mathrm{\Delta }_t`$ is an admissible exponent when $`\lambda _t=2tk+\mathrm{\Delta }_t`$ is permissible.

The computations required to determine sharp permissible exponents for a specific value of $`k`$ are substantial (see ), but for larger $`k`$ one may summarise some general features of these exponents. First, for $`0t2`$ and $`k2`$, it is essentially classical that the exponent $`\delta _t=0`$ is associated, and recent work of Heath-Brown  provides the same conclusion also when $`t=3`$ and $`k238,607,918`$. When $`t=o(\sqrt{k})`$, one finds that associated exponents exhibit quasi-diagonal behaviour, and satisfy the property that $`\delta _t0`$ as $`k\mathrm{}`$. To be precise, Theorem 1.3 of  shows that whenever $`k3`$ and $`2<t2e^1k^{1/2}`$, then the exponent

$$\delta _t=\frac{4k^{1/2}}{et}\mathrm{exp}\left(\frac{4k}{e^2t^2}\right),$$
$`2.1`$
is associated. For larger $`t`$, methods based on repeated efficient differencing yield the sharpest estimates. Thus, the corollary to Theorem 2.1 of  establishes that for $`k4`$, an admissible exponent $`\mathrm{\Delta }_t`$ is given by the positive solution of the equation $`\mathrm{\Delta }_te^{\mathrm{\Delta }_t/k}=ke^{12t/k}`$. The exponent $`\lambda _t=2tk+ke^{12t/k}`$ is therefore always permissible. Previous to repeated efficient differencing, analogues of these permissible exponents had a term of size $`ke^{t/k}`$ in place of $`ke^{12t/k}`$ (see ), so that in a sense, the modern theory is twice as powerful as that available hitherto.

The above discussion provides a useable analogue of the mean-value estimate in (1.3). We turn next to localised minor arc estimates. Take $`Q=P`$, and define $`𝔪`$ as in the introduction. Suppose that $`s`$, $`t`$ and $`w`$ are parameters with $`2sk+1`$ for which $`\mathrm{\Delta }_s`$, $`\mathrm{\Delta }_t`$ and $`\mathrm{\Delta }_w`$ are admissible exponents, and define

$$\sigma (k)=\frac{k\mathrm{\Delta }_t\mathrm{\Delta }_s\mathrm{\Delta }_w}{2(s(k+\mathrm{\Delta }_w\mathrm{\Delta }_t)+tw(1+\mathrm{\Delta }_s))}.$$

Then Corollary 1 to Theorem 4.2 of  shows that $`sup_{\alpha 𝔪}|h(\alpha )|P^{1\sigma (k)+ϵ}`$, and for large $`k`$ this estimate holds with $`\sigma (k)^1=k(\mathrm{log}k+O(\mathrm{log}\mathrm{log}k))`$. Applying an analogue of (1.2) with $`h`$ in place of $`f`$, and taking 333We write $`[z]`$ to denote $`\mathrm{max}\{n:nz\}`$. $`t=[\frac{1}{2}k(\mathrm{log}k+\mathrm{log}\mathrm{log}k+1)]`$ and $`s=2t+k+[Ak\mathrm{log}\mathrm{log}k/\mathrm{log}k]`$, for a suitable $`A>0`$, we deduce from our discussion of permissible exponents that $`_𝔪h(\alpha )^se(n\alpha )𝑑\alpha =o(n^{s/k1})`$. By considering the representations of a given integer $`n`$ with all of the $`k`$th powers $`R`$-smooth, it is now apparent that a modification of the argument sketched in the introduction shows that $`R(n)𝔖_{s,k}(n)n^{s/k1}`$ as soon as one confirms that

$$_𝔐h(\alpha )^se(n\alpha )𝑑\alpha c(\eta )^s\frac{\mathrm{\Gamma }(1+1/k)^s}{\mathrm{\Gamma }(s/k)}𝔖_{s,k}(n)n^{s/k1}.$$
$`2.2`$
Sharp asymptotic information concerning $`h(\alpha )`$ is available throughout $`𝔐(Q)`$ only when $`Q`$ is a small power of $`\mathrm{log}P`$, and so the proof of (2.2) involves pruning technology. Such machinery, in this case designed to estimate the contribution from a set of the shape $`𝔐(P)𝔐((\mathrm{log}P)^\delta )`$, has evolved into a powerful tool. Such issues can be handled these days with a number of variables barely exceeding $`\mathrm{max}\{4,k+1\}`$.

This approach leads to the best known upper bounds on the function $`G(k)`$ in Waring’s problem, defined to be the least integer $`r`$ for which all sufficiently large natural numbers are the sum of at most $`r`$ positive integral $`k`$th powers.

Theorem 2.1.  One has $`G(k)k(\mathrm{log}k+\mathrm{log}\mathrm{log}k+2+O(\mathrm{log}\mathrm{log}k/\mathrm{log}k))`$.

This upper bound (Theorem 1.4 of ) refines an earlier one of asymptotically similar strength (Corollary 1.2.1 of ) that gave the first sizeable improvement of Vinogradov’s celebrated bound $`G(k)(2+o(1))k\mathrm{log}k`$, dating from 1959 (see ). Aside from Linnik’s bound $`G(3)7`$ (see ), all of the sharpest known bounds on $`G(k)`$ for smaller $`k`$ are established using variants of these methods. Thus one has $`G^\mathrm{\#}(4)12`$ (see , and here the $`\mathrm{\#}`$ denotes that there are congruence conditions modulo $`16`$), $`G(5)17`$, $`G(6)24`$, $`G(7)33`$, $`G(8)42`$, $`G(9)50`$, $`G(10)59`$, $`G(11)67`$, $`G(12)76`$, $`G(13)84`$, $`G(14)92`$, $`G(15)100`$, $`G(16)109`$, $`G(17)117`$, $`G(18)125`$, $`G(19)134`$, $`G(20)142`$ (see , , , ).

Unfortunately, shortage of space obstructs any but the crudest account of the ideas underlying the proof of the mean value estimates that supply the above permissible exponents. The use of exponential sums over smooth numbers occurs already in work of Linnik and Karatsuba (see ), but only with Vaughan’s new iterative method  is a flexible homogeneous approach established. An alternative formulation suitable for repeated efficient differencing is introduced by the author in . Suppose that the exponent $`\lambda _s`$ is permissible, and consider a polynomial $`\psi [t]`$ of degree $`d2`$. Given positive numbers $`M`$ and $`T`$ with $`MT`$, and an element $`x𝒜(T,R)`$ with $`x>M`$, there exists an integer $`m`$ with $`m[M,MR]`$ for which $`m|x`$. Consequently, by applying a fundamental lemma of combinatorial flavour, one may bound the number of integral solutions of the equation

$$\psi (z)\psi (w)=\underset{i=1}{\overset{s}{}}(x_i^ky_i^k),$$
$`2.3`$
with $`1z,wP`$ and $`x_i,y_i𝒜(T,R)`$ $`(1is)`$, in terms of the number of integral solutions of the equation

$$\psi (z)\psi (w)=m^k\underset{i=1}{\overset{s}{}}(u_i^kv_i^k),$$
$`2.4`$
with $`1z,wP`$, $`M<mMR`$, $`(\psi ^{}(z)\psi ^{}(w),m)=1`$ and $`u_i,v_i𝒜(T/M,R)`$ $`(1is)`$. The implicit congruence condition $`\psi (z)\psi (w)(modm^k)`$ may be analytically refined to the stronger one $`zw(modm^k)`$, and in this way one is led to replace the expression $`\psi (z)\psi (w)`$ by the difference polynomial $`\psi _1(z;h;m)=m^k(\psi (z+hm^k)\psi (z))`$. Notice that when $`MP^{1/k}`$, one is forced to conclude that $`z=w`$, and then the number of solutions of (2.4) is bounded above by $`PMRS_s(T/M,R)P^{1+ϵ}M(T/M)^{\lambda _s}`$. Otherwise, following an application of Schwarz’s inequality to the associated mean value of exponential sums, one may recover an equation of the shape (2.3) in which $`\psi (z)`$ is replaced by $`\psi _1(z)`$, and $`T`$ is replaced by $`T/M`$, and repeat the process once again. This gives a repeated differencing process that hybridises that of Weyl with the ideas of Vinogradov.

It is now possible to describe a strategy for bounding a permissible exponent $`\lambda _{s+1}`$ in terms of a known permissible exponent $`\lambda _s`$. We initially take $`T=P`$ and $`\psi (z)=z^k`$, and observe that $`S_{s+1}(P,R)`$ is bounded above by the number of solutions of (2.3). We apply the above efficient differencing process successively with appropriate choices for $`M`$ at each stage, say $`M=P^{\varphi _i}`$, with $`0\varphi _i1/k`$, for the $`i`$th differencing operation. After some number of steps, say $`j`$, we take $`\varphi _j=1/k`$ in order to force the above diagonal situation that is easily estimated. One then optimises choices for the $`\varphi _i`$ in order to extract the sharpest upper bound for $`S_{s+1}(P,R)`$, and this in turn yields a permissible exponent $`\lambda _{s+1}`$. It transpires that in this simplified treatment, successive admissible exponents are related by the formula $`\mathrm{\Delta }_{s+1}=\mathrm{\Delta }_s(1\varphi )+k\varphi 1`$, wherein one may take $`\varphi `$ very close to $`1/(k+\mathrm{\Delta }_s)`$. Thus one finds that $`\mathrm{\Delta }_{s+1}`$ is essentially $`\mathrm{\Delta }_s(12/(k+\mathrm{\Delta }_s))`$, an observation that goes some way to explaining how it is that this method is about twice as strong as previous approaches that would correspond to choices of $`\varphi `$ close to $`1/k`$.

Refined versions of this differencing process make use of all known permissible exponents $`\lambda _s`$ in order to estimate a particular exponent $`\lambda _t`$, and in such circumstances the process becomes highly iterative, and entails significant computation. Such variants make use of refined Weyl estimates for difference polynomials, and estimates for the number of integral points on curves and surfaces (see ). Variants of these methods apply also in the situation of Vinogradov’s mean value theorem (see ), smooth Weyl sums with polynomial arguments (see ), and even for sums relevant to counting rational lines on hypersurfaces (see ).

Frequent reference to underlying diophantine equations seems to limit these methods to estimating even moments of smooth Weyl sums, and until recently fractional moments could be estimated only by applying Hölder’s inequality to interpolate linearly between permissible exponents. However, a method  is now available that permits fractional moments to be estimated non-trivially, thereby “breaking classical convexity”, and moreover the number of variables being differenced need not even be an integer. These new estimates can be applied to sharpen permissible exponents (with integral argument), and indeed the associated exponent (2.1) is established in this way. Another consequence  of these developments is the best available lower bound for $`N(X)`$, which we define to be the number of integers not exceeding $`X`$ that are represented as the sum of three positive integral cubes. One has $`N(X)X^{1\xi /3ϵ}`$, where $`\xi =(\sqrt{2833}43)/41=0.24941301\mathrm{}`$ arises from the permissible exponent $`\lambda _3=3+\xi `$ for $`k=3`$. Earlier, Vaughan  obtained an estimate of the latter type with $`13/4`$ in place of $`3+\xi `$.

## 3. Arithmetic variants of Bessel’s inequality

Already in our opening paragraph we alluded to some of the applications accessible to the methods of §2. We now turn to less obvious applications that have experienced recent progress. We illustrate ideas once again with a simple example, and consider the set $`𝒵(N)`$ of integers $`n`$, with $`N/2<nN`$, that are not represented as the sum of $`s`$ positive integral $`k`$th powers. The standard approach to estimating $`Z(N)=\text{card}(𝒵(N))`$ is via Bessel’s inequality. We now take $`P=N^{1/k}`$. When $`𝔅[0,1)`$, write $`R^{}(n;𝔅)=_𝔅h(\alpha )^se(n\alpha )𝑑\alpha `$, and write also $`R^{}(n)=R^{}(n;[0,1))`$. The theory of §2 ensures that when $`Q`$ is a sufficiently small power of $`\mathrm{log}P`$, and $`s4k`$, then $`R^{}(n;𝔐)n^{s/k1}`$. Under such circumstances, an application of Bessel’s inequality reveals that $`Z(N)`$ is bounded above by

$`{\displaystyle \underset{N/2<nN}{}}\left|{\displaystyle \frac{R^{}(n)R^{}(n;𝔐)}{R^{}(n;𝔐)}}\right|^2`$ $`(N^{s/k1})^2{\displaystyle \underset{n}{}}\left|{\displaystyle _𝔪}h(\alpha )^se(n\alpha )𝑑\alpha \right|^2`$ $`3.1`$
$`(N^{s/k1})^2{\displaystyle _𝔪}|h(\alpha )|^{2s}𝑑\alpha .`$
When $`s\frac{1}{2}k(\mathrm{log}k+\mathrm{log}\mathrm{log}k+2+o(1))`$, the minor arc integral in (3.1) is $`o(N^{2s/k1})`$, and thus it follows that $`Z(N)=o(N)`$. Thus one may conclude that almost all integers are sums of $`s(\frac{1}{2}+o(1))k\mathrm{log}k`$ positive integral $`k`$th powers.

The application of Bessel’s inequality in (3.1) makes inefficient use of underlying arithmetic information, and fails, for example, to effectively estimate the number of values of a polynomial sequence not represented in some prescribed form. Suppose instead that we define a Fourier series over the exceptional set itself, namely $`K(\alpha )=_ne(n\alpha )`$, where the summation is over $`n𝒵(N)`$. Since $`R^{}(n)=0`$ for $`n𝒵(N)`$, one has $`R^{}(n;𝔪)=R^{}(n;𝔐)`$, and thus we see that

$$N^{s/k1}Z(N)_𝔐h(\alpha )^sK(\alpha )𝑑\alpha =\left|_𝔪h(\alpha )^sK(\alpha )𝑑\alpha \right|.$$

Applying Schwarz’s inequality in combination with Parseval’s identity, we recover the previous consequence of Bessel’s inequality via the bound

$$\left|_𝔪h(\alpha )^sK(\alpha )𝑑\alpha \right|\left(_0^1|K(\alpha )|^2𝑑\alpha \right)^{1/2}\left(_𝔪|h(\alpha )|^{2s}𝑑\alpha \right)^{1/2}.$$
$`3.2`$
However, this formulation permits alternate applications of Schwarz’s inequality or Hölder’s inequality. For example, the left hand side of (3.2) is bounded above by

$$\left(_0^1|h(\alpha )^{2t}K(\alpha )^2|𝑑\alpha \right)^{1/2}\left(_𝔪|h(\alpha )|^{2s2t}𝑑\alpha \right)^{1/2},$$
$`3.3`$
and also by

$$\left(_0^1|K(\alpha )|^4𝑑\alpha \right)^{1/4}\left(_𝔪|h(\alpha )|^{4s/3}𝑑\alpha \right)^{3/4}.$$
$`3.4`$
In either case, the diophantine equations underlying the integrals on the left hand sides of (3.3) and (3.4) contain arithmetic information that can be effectively exploited whenever the set $`𝒵(N)`$ is reasonably thin.

The strategy sketched above has been exploited by Brüdern, Kawada and Wooley in a series of papers devoted to additive representation of polynomial sequences. Typical of the kind of results now available is the conclusion  that almost all values of a given integral cubic polynomial are the sum of six positive integral cubes. Also, Wooley , , has derived improved (slim) exceptional set estimates in Waring’s problem when excess variables are available. For example, write $`E(N)`$ for the number of integers $`n`$, with $`1nN`$, for which the anticipated asymptotic formula fails to hold for the number of representations of an integer as the sum of a square and five cubes of natural numbers. Then in  it is shown that $`E(N)N^ϵ`$.

As a final illustration of such ideas, we highlight an application to the solubility of pairs of diagonal cubic equations. Fix $`k=3`$, define $`h(\alpha )`$ as in §2, and put $`c(n)=_0^1|h(\alpha )|^5e(n\alpha )𝑑\alpha `$ for each $`n`$. Brüdern and Wooley  have applied the ideas sketched above to estimate the frequency with which large values of $`|c(n)|`$ occur, and thereby have shown that, with $`\xi `$ defined as in the previous section,

$$\underset{x,y𝒜(P,R)}{}|c(x^3y^3)|^2=_0^1_0^1|h(\alpha )^5h(\beta )^5h(\alpha +\beta )^2|𝑑\alpha 𝑑\beta P^{6+\xi +ϵ}.$$

On noting that $`6+\xi <6.25`$, cognoscenti will recognise that this twelfth moment of smooth Weyl sums, in combination with a classical exponential sum equipped with Weyl’s inequality, permits the discussion of pairs of diagonal cubic equations in $`13`$ variables via the circle method. The exponent $`6+\xi `$ improves an exponent $`6+2\xi `$ previously available for a (different) twelfth moment. Brüdern and Wooley  establish the following conclusion.

Theorem 3.1. Suppose that $`s13`$, and that $`a_i,b_i`$ $`(1is)`$ are fixed integers. Then the Hasse principle holds for the pair of equations

$$a_1x_1^3+\mathrm{}+a_sx_s^3=b_1x_1^3+\mathrm{}+b_sx_s^3=0.$$

The condition $`s13`$ improves on the previous bound $`s14`$ due to Brüdern , and achieves the theoretical limit of the circle method for this problem.

## 4. Arithmetic geometry via descent

Let $`F(𝕩)[x_1,\mathrm{},x_s]`$ be a homogeneous polynomial of degree $`d`$, and consider the number, $`N(B)`$, of integral zeros of the equation $`F(𝕩)=0`$, with $`𝕩[B,B]^s`$. When $`s`$ is sufficiently large in terms of $`d`$, the circle method shows under modest geometric conditions that $`N(B)`$ is asymptotic to the expected product of local densities. For fairly general polynomials, the condition on $`s`$ is as severe as $`s>(d1)2^d`$, though for diagonal equations the methods of §2 relax this condition to $`s>(1+o(1))d\mathrm{log}d`$. However, there is a class of varieties with small dimension relative to degree, for which the circle method supplies non-trivial information concerning the density of rational points. The idea is to apply a descent process in order to interpret points on the original variety in terms of corresponding points on a new variety, with higher dimension relative to degree, more amenable to the circle method.

To illustrate this principle, consider a field extension $`K`$ of $``$ of degree $`n`$ with associated norm form $`N(𝕩)[x_1,\mathrm{},x_n]`$. Also, let $`l`$ and $`k`$ be natural numbers with $`(k,l)=1`$, and let $`\alpha `$ be a non-zero rational number. Then Heath-Brown and Skorobogatov  descend from the variety $`t^l(1t)^k=\alpha N(𝕩)`$ to the associated variety $`aN(𝕦)+bN(𝕧)=z^n`$, for suitable integers $`a`$ and $`b`$. The circle method establishes weak approximation for the latter variety, and thereby it is shown that the Brauer-Manin obstruction is the only possible obstruction to the Hasse principle and weak approximation on any smooth projective model of the former variety. One can artificially construct further examples amenable to the circle method. For example, if we take linearly independent linear forms $`L_i(𝕩)[x_1,\mathrm{},x_n]`$ $`(1in+r)`$, then one can establish non-trivial lower bounds for the density of rational points on the variety $`z^k=L_1(𝕩)\mathrm{}L_{n+r}(𝕩)`$ by descending to a variety that resembles a system of $`r`$ diagonal forms of degree $`k`$, with constrained varying coefficients. The investigation of such matters will likely provide an active area of research into the future. In this context we point to work of Peyre , which addresses the interaction between descent and the circle method in some generality.

## References


format: text

text:
# 1 Introduction

## 1 Introduction

With the advance of electronics, data acquisition (DAQ) and computing capabilities, it is more and more common to use Flash Analog to Digital Convertors (FADCs) as the general-purpose readout device for different detectors. Typically, the entire pulse shape from the individual detector module, as well as the relative timing between pulses from different detector channels, are recorded. Such readout scheme replaces Analog Digital Convertors and Time-to-Digital Convertors, which give information only on the amplitude or integrated area and the timing of the pulses, respectively.

With FADCs as readout, the information related to a single pulse is represented by many measurements. Frequently, multi-channel (order of hundreds and thousands is common) FADC systems are deployed. As a result, a large volume of data have to be extracted, stored and processed. The n-bit resolution per time-bin and the sampling rate of the FADC are important parameters which will have tremendous effect on the characterization and performance of the experiment such as detector resolution, DAQ dead time, data size, data processing time, cost and so on.

A motivated question to address is therefore “What is the dynamic range of an n-bit FADC, beyond the nominal $`2^\mathrm{n}`$ value?” That is, can the crucial information like the total area as well as the rise and fall times be reconstructed, even though part of the pulse $``$ the largest and usually the most important part $``$ is not recorded? These are the focuses of investigation of this article.

In Section 2, we describe the research program on which this project is based, and the hardware used to generate the data which allow such an investigation. Two approaches are surveyed and compared. One of the algorithms, the Charge Summation Correction Method, showed better potentials and were pursued in depth. The goals are to achieve both energy reconstruction and particle identification even for saturated events. The details and results are discussed in Section 3.

Though the algorithm was developed and tested with scintillating CsI(Tl) crystal detectors, the approach is generic and can be applicable to other detector schemes using FADC readout as well.

## 2 Scintillating CsI(Tl) Crystal Detector

The present investigation is part of the research program of the TEXONO Collaboration . A CsI(Tl) scintillating crystal  array with 2 kg modular mass and a total mass of 186 kg is taking data at the Kuo-Sheng (KS) Nuclear Power Plant to study low energy neutrino physics. Each module is read out by two photo-multipliers (PMTs) and the signals are recorded by a 20 MHz FADCs with 8-bit resolution . The FADC value at the $`\mathrm{t}_\mathrm{i}^{\mathrm{th}}`$ time-bin is denoted by $`\mathrm{y}_\mathrm{i}`$. Throughout this work, we use the normalized value $`\mathrm{Y}_\mathrm{i}`$ which is $`\mathrm{y}_\mathrm{i}`$ divided by the saturation level $`\mathrm{L}_{\mathrm{sat}}=2^\mathrm{n}1=255`$ : $`\mathrm{Y}_\mathrm{i}=\mathrm{y}_\mathrm{i}/\mathrm{L}_{\mathrm{sat}}`$ .

The physics requirements of the experiment are that events with energy from 10 keV to 10 MeV can be reconstructed, and that pulse shape discrimination (PSD) can be applied to these events to achieve event ($`\alpha `$/$`\gamma `$) identification. That is, a dynamic range of at least 1000 is necessary with an 8-bit resolution readout per time-bin.

Special data were taken to devise techniques and to develop software algorithms addressing this issue. A $`5\times 5\times 5\mathrm{cm}^3`$ CsI(Tl) crystal was used and signals were read out by a single PMT and a FADC channel at identical settings to the KS experiment. The trigger was set at the $`\mathrm{t}_0=100^{\mathrm{th}}`$ FADC time-bin, such that the pulse information of $`\delta \mathrm{t}_{}`$=5 $`\mu `$s (100 time-bins) and $`\delta \mathrm{t}_+`$=25.6$`\mu `$s (512 time-bins) before and after the trigger were recorded, respectively. Both $`\alpha `$ and $`\gamma `$ events were taken by activating the CsI(Tl) crystal with 241Am-Be and 137Cs sources, respectively.

The amplitude of the pulse was adjusted by selecting the appropriate PMT high voltage. Typical unsaturated and saturated pulses are displayed in Figures 1a and 1b, respectively. A large sample of unsaturated events are added together to produce a smooth reference pulse shape $`\overline{\mathrm{Y}}`$(t), where the peak amplitude is normalized to 1. The profiles of $`\overline{\mathrm{Y}}`$(t) due to $`\alpha `$\- and $`\gamma `$-events are depicted in Figure 2, and can be fitted to the function

$$\overline{\mathrm{Y}}(\mathrm{t})=\mathrm{N}[1\mathrm{exp}(\frac{\mathrm{t}}{\tau _0})][\frac{1}{\tau _1}\mathrm{exp}(\frac{\mathrm{t}}{\tau _1})+\frac{\mathrm{r}}{\tau _2}\mathrm{exp}(\frac{\mathrm{t}}{\tau _2})].$$
(1)
The rise time($`\tau _0`$) and fall times($`\tau _1,\tau _2`$) as well as the ratio between the slow and fast decay components(r) are listed in Table 1. The normalization constant N is chosen to make the peak amplitude equal to unity. The values of $`\tau _0`$ are dominated by the electronics shaping rise time of 250 ns , while ($`\tau _1,\tau _2`$) are due to the characteristic fluorescence time profile of CsI(Tl).

The integrated charge (Q) of an event is obtained by numerically integrating the pedestal-subtracted FADC pulses:

$$\mathrm{Q}=\underset{\mathrm{t}_\mathrm{i}>\mathrm{t}_0}{}[\mathrm{Y}_\mathrm{i}\mathrm{Y}_0]=\underset{\mathrm{t}_\mathrm{i}>\mathrm{t}_0}{}\frac{[\mathrm{y}_\mathrm{i}\mathrm{y}_0]}{\mathrm{L}_{\mathrm{sat}}},$$
(2)
where $`\mathrm{y}_0`$ is the average pedestal level determined from the pre-trigger part ($`\mathrm{t}_\mathrm{i}<\mathrm{t}_0`$) of the pulse, and $`\mathrm{Y}_0=\mathrm{y}_0/\mathrm{L}_{\mathrm{sat}}`$. The typical full-width-half-maximum (FWHM) resolution for unsaturated events at 662 keV is about 8%.

## 3 Charge Summation Method

### 3.1 Initial Survey

Saturated events are those with pulse amplitude larger than those corresponding to the maximum FADC value of $`\mathrm{L}_{\mathrm{sat}}=255`$. The focus of this section is to devise an algorithm with which energy reconstruction and event identification by PSD can be achieved for these class of events.

Two approaches were surveyed as initial studies:

1. Statistical Approach In principle, one can derive the values of N in Equation 1, using maximum likelihood analysis or minimum $`\chi ^2`$ fitting methods, based only on the ($`t_i`$,$`\mathrm{y}_\mathrm{i}`$) data where $`\mathrm{y}_\mathrm{i}<\mathrm{L}_{\mathrm{sat}}`$. The values of $`(\tau _0,\tau _1,\tau _2,\mathrm{r})`$ are adopted from those in Table 1. In practice, the fits are computer-time consuming and unstable. The efficiencies are bad, especially when the saturation is severe. The main reason is that the values of $`\tau _0`$ and $`\tau _1`$ are derived from small portions of the pulse shape. When the data from these portions are missing due to saturation, there will be insufficient constraints such that it is difficult for the analysis procedures to converge to the optimal values for N.
2. Pulse Analysis Approach The objective of this approach is to identify measure-able values which are related to the physics variables of interest, namely, the total energy and the decay time. The observables are the “saturated charge” and “mean time”. This method uses all the available information from the saturated pulses. Such an approach can be compared to $``$ and in fact inspired by $``$ the “Double Charge Method”  commonly-used to achieve $`\alpha `$/$`\gamma `$ identification in liquid  and crystal  scintillators.

Results of the initial surveys indicated that the pulse analysis approach would show better promises. The “Charge Summation Method” based on this approach was therefore pursued in details. We discuss the algorithm, its consistency checks and performance in energy reconstruction and event identification in the following sub-sections.

### 3.2 Reference Profile and Correction Function

From a standard reference pulse shape, a relationship is obtained between the “saturated charge ($`\mathrm{Q}_{\mathrm{sat}}`$)” defined by the integration of the saturated events and the “total charge ($`\mathrm{Q}_{\mathrm{tot}}`$)” from that of the raw fully-recorded pulses.

The saturation level for an n=8-bit FADC is $`\mathrm{L}_{\mathrm{sat}}=2^\mathrm{n}1=255`$. To investigate the saturation effects on a pulse shape with amplitude $`\mathrm{\Lambda }>\mathrm{L}_{\mathrm{sat}}`$, a profile $`\mathrm{Z}(\mathrm{\Lambda })`$ is constructed from the reference shape $`\overline{\mathrm{Y}}(\mathrm{t})`$ of Figure 2, scaled-up by a factor of of $`\mathrm{\Lambda }`$/$`\mathrm{L}_{\mathrm{sat}}`$. Therefore, the value at the $`\mathrm{t}_\mathrm{i}^{\mathrm{th}}`$ time-bin is given by

$$\mathrm{Z}_\mathrm{i}(\mathrm{\Lambda })=\overline{\mathrm{Y}}_\mathrm{i}[\frac{\mathrm{\Lambda }}{\mathrm{L}_{\mathrm{sat}}}].$$
(3)
The pedestal level $`\mathrm{Z}_0`$ is the average of $`\mathrm{Z}_\mathrm{i}`$ for $`\mathrm{t}_\mathrm{i}<\mathrm{t}_0`$. The dependence of the total and saturated charge as functions of $`\mathrm{\Lambda }`$ for $`\gamma `$-events can be defined as:

$$\mathrm{Q}_{\mathrm{tot}}(\mathrm{\Lambda })=\underset{\mathrm{t}_\mathrm{i}}{}[\mathrm{Z}_\mathrm{i}\mathrm{Z}_0]$$
(4)
and

$$\mathrm{Q}_{\mathrm{sat}}(\mathrm{\Lambda })=\underset{\mathrm{t}_\mathrm{i}}{}[\xi _\mathrm{i}\mathrm{Z}_0],\mathrm{where}\{\begin{array}{cc}\xi _\mathrm{i}=\mathrm{Z}_\mathrm{i}\mathrm{for}\mathrm{Z}_\mathrm{i}<1\hfill & \\ \xi _\mathrm{i}=1\mathrm{for}\mathrm{Z}_\mathrm{i}1\hfill & \end{array},$$
(5)
respectively. That is, given the profile $`\mathrm{Z}(\mathrm{\Lambda })`$ with peak amplitude $`\mathrm{\Lambda }`$/$`\mathrm{L}_{\mathrm{sat}}`$, $`\mathrm{Q}_{\mathrm{tot}}`$ and $`\mathrm{Q}_{\mathrm{sat}}`$ are the integrated areas of the complete pulse and the pulse truncated at a maximum amplitude of $`\mathrm{\Lambda }=\mathrm{L}_{\mathrm{sat}}`$, respectively. By varying the values of $`\mathrm{\Lambda }/\mathrm{L}_{\mathrm{sat}}`$, the dependence of $`\mathrm{Q}_{\mathrm{tot}}`$ and $`\mathrm{Q}_{\mathrm{sat}}`$ on $`\mathrm{\Lambda }/\mathrm{L}_{\mathrm{sat}}`$ for $`\gamma `$\- and $`\alpha `$-events can be derived. They are depicted in Figure 3a and 3b, respectively. As expected, the variations shows a linear relationship between $`\mathrm{Q}_{\mathrm{tot}}`$ and $`\mathrm{\Lambda }/\mathrm{L}_{\mathrm{sat}}`$, while $`\mathrm{Q}_{\mathrm{sat}}`$ would level off at large $`\mathrm{\Lambda }/\mathrm{L}_{\mathrm{sat}}`$ due to saturation.

A relationship between $`\mathrm{Q}_{\mathrm{tot}}`$ and $`\mathrm{Q}_{\mathrm{sat}}`$ for the saturated events (that is, when $`\mathrm{Q}_{\mathrm{tot}}>\mathrm{Q}_{\mathrm{sat}}`$) can then be obtained. It is described by a polynomial function

$$\mathrm{Q}_{\mathrm{tot}}=\underset{\mathrm{k}=0}{\overset{3}{}}\alpha _\mathrm{k}\mathrm{Q}_{\mathrm{sat}}^{}{}_{}{}^{\mathrm{k}}$$
(6)
for both the $`\gamma `$\- and $`\alpha `$-samples as displayed in Figure 4. The coefficients $`\alpha _\mathrm{i}`$ are obtained from minimum $`\chi ^2`$-fits.

This function is then used to compensate the missing area in saturated events. Integration of the saturated pulses provide the values of $`\mathrm{Q}_{\mathrm{sat}}`$, from which the correct values of the total charge $`\mathrm{Q}_{\mathrm{tot}}`$ can be derived in an event-by-event basis.

### 3.3 Consistency Tests

The correction algorithm discussed in the previous section does not depend on the exact saturation level $`\mathrm{L}_{\mathrm{sat}}`$. One can use this feature to perform consistency checks on data sample where the total charge is known. An unsaturated event is made saturated by software. The algorithm is applied to the simulated saturated event, and the corrected total charge is compared with the raw measurements.

A level of $`\mathrm{L}_{\mathrm{test}}`$=100 is chosen. A new set of $`\alpha _i`$ is obtained by going through Eqs. 3 to 6 again, with $`\mathrm{L}_{\mathrm{sat}}`$ replaced by $`\mathrm{L}_{\mathrm{test}}`$. A data set of $`\gamma `$-events was taken with a 137Cs source. The PMT high-voltage was adjusted such that the amplitude for the 662 keV line is below the actual saturation level of 255 but significantly above $`\mathrm{L}_{\mathrm{test}}`$. The charge $`\mathrm{Q}`$ is measured by integrating the normalized FADC counts $`\mathrm{Y}_\mathrm{i}`$, and the spectrum is displayed in Figure 5a. The same data is then made saturated by re-defining $`\mathrm{Y}_\mathrm{i}\mathrm{Y}_\mathrm{i}^{}`$ such that

$$\mathrm{Y}_\mathrm{i}^{}=\{\begin{array}{cc}\mathrm{Y}_\mathrm{i}\mathrm{for}\mathrm{Y}_\mathrm{i}\mathrm{L}_{\mathrm{test}}/\mathrm{L}_{\mathrm{sat}}\hfill & \\ \mathrm{L}_{\mathrm{test}}/\mathrm{L}_{\mathrm{sat}}\mathrm{for}\mathrm{Y}_\mathrm{i}>\mathrm{L}_{\mathrm{test}}/\mathrm{L}_{\mathrm{sat}}\hfill & \end{array}.$$
(7)
The saturated charge is derived by integrating the normalized $`\mathrm{Y}_\mathrm{i}^{}`$ and the spectra is shown in Figure 5b, while that for the corrected charge $`\mathrm{Q}^{}`$ after the correction procedures is shown in Figure 5c. There is excellent agreement on the 137Cs peaks in terms of their integrated charge, resolution and area between Figures 5a and 5c. The event-by-event deviations are parametrized by

$$\mathrm{\Delta }=\frac{\mathrm{Q}\mathrm{Q}^{}}{(\mathrm{Q}+\mathrm{Q}^{})/2}.$$
(8)
The distribution of $`\mathrm{\Delta }`$ is depicted in Figure 6. The results show that the mean bias is less than 0.5% which is small compared to the intrinsic FWHM resolution of 8.7%. Therefore, the correction procedures are valid and justified.

### 3.4 Results

#### 3.4.1 Energy Reconstruction

The algorithms are applied to actual saturated data. The data is generated by increasing the PMT high voltages such that the events corresponding to the 662 keV $`\gamma `$-line from 137Cs have amplitude larger than $`\mathrm{L}_{\mathrm{sat}}`$=255, while those from the 32 keV line remain unsaturated.

An example of the energy spectra of the saturated events, as well as the corrected spectra are illustrated in Figure 7a and Figure 7b, respectively. The ratio of the peak positions at 32 keV and 662 keV (denoted by $`\mathrm{Q}_{32}`$ and $`\mathrm{Q}_{662}`$, respectively) is measured and compared to the nominal values:

$$\mathrm{R}=\frac{\mathrm{Q}_{32}/\mathrm{Q}_{662}}{32/662}.$$
(9)
This can be used as a figure of merit to quantify the performance. The values of R for saturated and corrected events as functions of the expected peak amplitudes are displayed in Figure 8a. It can be seen that the charge summation correction methods extend the effective dynamic range of the FADC by at least a factor of 14, or equivalently, 4 more bits.

The dependence of performance of the algorithm on the FADC recording period $`\delta \mathrm{t}_+`$ is studied. The data are truncated by software to a shorter interval and the same correction procedures are repeated. The results are displayed in Figure 8b for the variation of R with peak amplitudes at $`\delta \mathrm{t}_+=25.6,12.8,6.4\mu \mathrm{s}`$ (that is, 512, 256 and 128 time-bins, respectively). It can be seen that the charge reconstruction is no longer valid for $`\delta \mathrm{t}_+=6.4\mu \mathrm{s}`$ at an amplitude of $`>`$8 times of $`\mathrm{L}_{\mathrm{sat}}`$. Alternatively, a plot of R versus $`\delta \mathrm{t}_{\mathrm{unsat}}`$/$`\delta \mathrm{t}_+`$ is displayed in Figure 9, where $`\delta \mathrm{t}_{\mathrm{unsat}}`$ is the unsaturated interval, where the FADC counts are less than $`\mathrm{L}_{\mathrm{sat}}`$. It can be seen the values of R deviate from 1 only for $`\delta \mathrm{t}_{\mathrm{unsat}}`$/$`\delta \mathrm{t}_+<0.02`$. Accordingly, the algorithm is valid so long as the the unsaturated interval of the event is more than 2% of the FADC recording period. In comparison, the fitting method as discussed in Section 3.1 would start to fail when the saturation interval is more than twice the value of $`\tau _1`$, or equivalently, the figure-of-merit R deviates from 1 when $`\delta \mathrm{t}_{\mathrm{unsat}}`$/$`\delta \mathrm{t}_+<0.74`$.

The dependence of energy resolution of the 662 keV line on amplitudes is shown in Figure 8c. There is no intrinsic bias in the correction, though the energy resolution gets worse from around 9% to 16% FWHM when the amplitude increases by a factor of 14 at $`\delta \mathrm{t}_+`$=25.6 $`\mu `$s. The degradation in energy resolution is larger for the shorter FADC recording time.

#### 3.4.2 Pulse Shape Discrimination

As depicted in Figure 2, the scintillation light in CsI(Tl) has different decay time constants for $`\gamma `$\- and $`\alpha `$-activated events. The $`\alpha `$-events decay faster in inorganic crystal scintillators $``$ opposite to the responses in organic liquid scintillators . Accordingly, the pulse shape information recorded by the FADC can provide event identification capabilities.

Identification of $`\gamma `$/$`\alpha `$ events is achieved by studying and comparing the different responses from two observables. One of them, the saturated charge ($`\mathrm{Q}_{\mathrm{sat}}`$), has been adopted in energy re-construction as discussed in Section 3. The other is the “mean time” ($`<\mathrm{t}>`$) which reflects the decay profiles, defined as:

$$<\mathrm{t}>=\frac{\underset{\mathrm{t}_\mathrm{i}>\mathrm{t}_0}{}[\mathrm{t}_\mathrm{i}\mathrm{t}_0][\mathrm{Y}_\mathrm{i}\mathrm{Y}_0]}{\underset{\mathrm{t}_\mathrm{i}>\mathrm{t}_0}{}[\mathrm{Y}_\mathrm{i}\mathrm{Y}_0]}.$$
(10)
The average variations of $`<\mathrm{t}>`$ with $`\mathrm{Q}_{\mathrm{sat}}`$ for both $`\gamma `$\- and $`\alpha `$-events can be derived from the reference pulse shape of Figure 2, and are displayed in Figure 10a. The scattered plot of $`<\mathrm{t}>`$ versus $`\mathrm{Q}_{\mathrm{sat}}`$ for individual events is shown in Figure 10b. It can be seen that $`\alpha `$/$`\gamma `$ events are well separated for $`\mathrm{Q}_{\mathrm{sat}}120`$. From Figure 3, it can be concluded that $`\alpha `$/$`\gamma `$ differentiation can be achieved for $`\alpha `$-events with amplitude up to 20 times of the saturation level.

## 4 Summary

An algorithm to process saturated events in Flash Analog-to-Digital Convertor was developed and studied. Charge reconstruction and event identification via pulse shape analysis are achieved, effectively extending an 8-bit FADC by at least 4 more bits. The algorithm is valid so long as the unsaturated interval is more than 2% of the FADC recording period.

Although this work originates from the goals of handling data from CsI(Tl) crystal scintillators, the approach should be generic and can be applicable to other detector systems where the entire pulse shape are recorded while part of which are saturated. These studies will facilitate the specification and design of large FADC-based data acquisition systems.

This work was supported by contracts 90-2112-M-001-037 and 91-2112-M-001-036 from the National Science Council, Taiwan, and 19975050 from the National Science Foundation, China.


text:
# Spin readout and initialization in a semiconductor quantum dot

## Abstract

Electron spin qubits in semiconductors are attractive from the viewpoint of long coherence times. However, single spin measurement is challenging. Several promising schemes incorporate ancillary tunnel couplings that may provide unwanted channels for decoherence. Here, we propose a novel spin-charge transduction scheme, converting spin information to orbital information within a single quantum dot by microwave excitation. The same quantum dot can be used for rapid initialization, gating, and readout. We present detailed modeling of such a device in silicon to confirm its feasibility.

Quantum dot qubits lie at the point of intersection between quantum computing and few-electron spintronics. The long coherence times and novel quantum mechanical properties of electron spins, the attraction of scalability, and the shared technology base with the digital electronics industry fuel a large research effort to develop these systems. Several proposals have been put forward for spin qubits in semiconductors loss98 ; kane98 ; vrijen ; levy ; friesen03 . In particular, exchange coupled gates in semiconductor dots appear to be workable friesen03 , and several components of qubit technology have been demonstrated recently elzerman . However, the combined challenge of preparing, storing and measuring spins is formidable.

Measurements of spin qubits pose a special challenge. On the one hand, qubits should be well isolated from their environment to avoid decoherence. On the other hand, it is necessary to individually couple the qubits to an external measurement device. Qubit initialization involves an additional dissipative coupling to the environment. For quantum computing, we must initiate such couplings selectively, and with sufficient strength to perform the operations quickly. Indeed, scalable quantum computing relies on fault-tolerant quantum error correction algorithms, involving frequent, parallel measurements, and a steady supply of initialized qubits Sho96 . Fortunately, rapid and sensitive quantum measurement techniques involving radio frequency single electron transistors (rf-SETs) have been developed scholkopf . rf-SETs have been used to detect the tunneling of individual electrons in semiconductor devices rimberg03 .

We have previously designed a quantum dot architecture specifically for the purpose of manipulating electron spins for fast and accurate two-qubit operations that serve as universal gates for quantum computations friesen03 . Recent experimental results have shown that decoherence does not pose a fundamental problem for such gate operations lyon03 . Using special qubit geometries friesen02 , it should be possible to perform reliable gate operations in silicon quantum dots at rates between about 1 MHz and 1 GHz note . The question is now whether similar speeds and reliable operation can be achieved for measurement and initialization operations. In this paper, we show how to extend our architecture to enable the readout and initialization of the state of a single spin. The general scheme does not assume a specific material system. However, in this Letter we focus upon silicon-based devices because of their desirable coherence properties.

Our method relies on the idea of converting spin information to charge information discussed by Loss and DiVincenzo loss98 , and using single electron transistors to read out the resulting spin state as proposed by Kane et al. kane98 , who posit spin-dependent charge motion onto impurities in Si. Martin et al. martin  have proposed a scheme for single spin readout that also converts spin information into charge information in an electron trap near a conducting channel. The resistance of the channel depends on the occupation of the trap, which in turn can be made to depend on the spin. Our proposed device has the additional feature in that it enables both readout and rapid initialization of the spin state. Rapid initialization (as opposed to initialization by thermalization) is just as essential as readout, and just as difficult to carry out. Our design is simple in that it obviates the need for spin-polarized leads or ancillary qubits. We go beyond conceptual design–the full three-dimensional system is simulated self-consistently and the operating parameters are thereby optimized. This enables us to demonstrate quantitatively that the device is in a working regime.

The basic Si-SiGe heterostructure is described in Ref. friesen03 . The main point for present purposes is that the active layer is pure strained Si, which minimizes decoherence from spin-phonon coupling tahan02 . The qubits are gated quantum dots which hold one electron, and the gate geometry is the key to our scheme. Specifically, the electrons are confined in asymmetric lateral wells, such that orbital excitation results in lateral center-of-charge movement. Figure 1 shows the first two orbital states of an electron confined to an asymmetric box, with center-of-charge positions varying by a distance $`\mathrm{\Delta }y`$. If an external radiation source of frequency $`E_{12}/h`$ drives the system between the ground and first excited states, $`n=1`$ and $`n=2`$, the center-of-charge of the electron will oscillate in time at the Rabi frequency $`\nu _R`$. This motion can be detected by a sensitive electrometer, such as a single electron transistor or a quantum point contact. With the system placed in a magnetic field, charge motion can be generated in a *spin-dependent* fashion. We define the logical 0 and 1 of the qubit as the spin up and spin down states of the electron in its orbital ground state, and use the orbital excited state only during initialization and readout. The system is now driven at the readout frequency $`\nu _{12}(B)=E_{12}/hg\mu _BB/h`$. This causes charge motion only if the qubit was in the down state at the time of measurement. Although this transition is forbidden at the electric dipole level because of the spin flip, spin-orbit coupling allows for a nonzero transition rate, as described below.

Crucially, this architecture also allows for rapid initialization of the qubit. Consider exposing a random ensemble of qubits to radiation of frequency $`\nu _{12}`$. An electron in the $`(1,)`$ state will be excited to the level $`(2,)`$ and experience a relatively fast relaxation to the ground state $`(1,)`$, compared to a spin-flip relaxation to the level $`(1,)`$. The net result is to rapidly polarize and thereby initialize the qubits. By varying the gate voltages (and thus $`\nu _{12}`$) on individual dots we ensure that only desired qubits are brought into resonance. Clearly, understanding the competition between the three time scales, $`1/\nu _R`$, $`1/\mathrm{\Gamma }_I`$, and $`T_1`$, as a function of material parameters and gate potentials is the key to utilizing this device for readout and initialization. Note that $`T_1`$ represents a thermal initialization time from $`(1,)`$ to $`(1,)`$. Since $`T_1`$ is a decoherence time, we require $`1/\nu _RT_1`$. The Rabi oscillation frequency $`\nu _R`$ depends on the incident intensity and is therefore controllable, within limits. Robust measurement requires that many Rabi oscillations occur before orbital decay: $`1/\nu _R1/\mathrm{\Gamma }_I`$.

In order to understand the competition of time scales, we introduce an rms interaction energy, time averaged (as indicated by {}) over an optical cycle, which expresses the strength of the interaction in the electric dipole approximation. This defines the Rabi oscillation frequency shore : $`|h\nu _R|^2=\{|V^{E1}|^2\}`$. Here, $`V^{E1}=(e\mathrm{}E_0/m^{}E_{12})\widehat{ϵ}𝐩`$ is the dipole term in the Hamiltonian, $`|E_0|^2`$ is twice the mean value of $`|𝐄(t)|^2`$ averaged in time, and $`\widehat{ϵ}`$ is the polarization unit vector. The electric field $`E_0`$ inside the semiconductor with dielectric constant $`\epsilon =\epsilon _r\epsilon _0`$ is related to the intensity of the external radiation $`I`$ by $`E_0=\sqrt{2I/c\epsilon _0\sqrt{\epsilon _r}}`$.

The dipole Hamiltonian does not flip the spin directly, but spin-orbit coupling causes each qubit state to be an additive mixture of up and down spin. In the 2D limit, the spin-orbit (so) Hamiltonian is dominated by the bulk \[Dresselhaus (D)\] and structural \[Rashba (R)\] inversion asymmetry terms, $`H_{\mathrm{so}}=H_D+H_R`$, where $`H_D=\beta (p_y\sigma _yp_x\sigma _x)`$ and $`H_R=\alpha (p_x\sigma _yp_y\sigma _x)`$. Including $`H_{so}`$ perturbatively gives a nonzero dipole matrix element, and for light polarized in the $`y`$ direction, the readout frequency is given by

$$\left|\nu _R\right|\frac{eE_0}{2\pi h\nu _{12}(B)}\sqrt{\alpha ^2+\beta ^2}|2|y_y|1|.$$
(1)
Note that the applied radiation need not be circularly polarized for this readout scheme. The Dresselhaus and Rashba parameters $`\alpha `$ and $`\beta `$ depend on intrinsic material properties, device design, and external electric field. For GaAs, both theoretical and experimental values vary widely silva : $`\alpha =11000`$ m/s and $`\beta =10003000`$ m/s. In a centrosymmetric crystal like silicon which has no bulk inversion asymmetry, $`\beta =0`$. The one known data point for a SiGe two-dimensional electron gas gives $`\alpha 8`$ m/s alphainsi-rossler , which we use in our estimates below.

The relaxation of our quantum dot to its ground state enables spin polarization, but limits or even inhibits readout if it occurs too quickly. This problem has been addressed by Khaetskii and Nazarov in GaAs quantum dots khaetskii . In silicon, where there is no piezoelectric interaction, we calculate the relaxation rate via the Golden Rule with the usual deformation potential electron-phonon Hamiltonian ridley . At sufficiently low temperatures $`(T<1K)`$, optical polar phonons and multiphonon processes do not contribute. By considering only longitudinal phonons, with dispersion $`\omega =v_lq`$, and using the long-wavelength approximation $`e^{i𝐤𝐫}1+i𝐤𝐫`$, we obtain the orbital decay rate due to electron-lattice coupling:

$$\mathrm{\Gamma }_I=\frac{(E_{12})^5}{6\pi \mathrm{}^6v_l^7\rho }\left(\mathrm{\Xi }_d+\mathrm{\Xi }_u/3\right)^2\underset{i}{}|1|x_i|2|^2,$$
(2)
where $`\rho `$ is the mass density, and $`\mathrm{\Xi }_d`$ and $`\mathrm{\Xi }_u`$ are the deformation constants. In strained systems, transverse phonons can also be important, as we shall discuss in a future publication friesen04 .

We perform a numerical analysis to obtain performance characteristics for the proposed measurement scheme. The numerical techniques used are an extension of those used in Refs. \[friesen03, ; friesen02, \]. The gate potentials, the electronic orbitals, and their corresponding image potentials (arising predominantly from the metallic gates) are computed self-consistently by a combination of three-dimensional finite element and diagonalization techniques. Specifically, we determine the readout oscillation frequency Eq. (1), the orbital decay rate Eq. (2), and the coupling sensitivity of the qubit electron to an integrated SET. We have modeled numerous devices and a full account will be published elsewhere friesen04 . One of the most promising is shown in Fig. 2. This is a 6 nm strained silicon quantum well sandwiched between layers of strain-relaxed silicon-germanium (Si85Ge15). The bottom barrier (30 nm) separates the quantum well from a grounded metallic back-gate. The top barrier (30 nm) separates the quantum well from lithographically patterned Schottky top-gates, whose voltages can be controlled independently. We consider negative gate potentials, which provide lateral confinement of the quantum dot through electrostatic repulsion. The main features of the quantum dot are that it is narrow, long, and slightly asymmetric. The narrow feature ensures that the excited orbitals are non-degenerate, so a microwave field with narrow line width will not induce unwanted transitions. The long feature ensures that the energy splitting $`E_{12}`$ (and thus $`\mathrm{\Gamma }_I`$) will be small. The asymmetry ensures that excited orbitals will provide charge motion. Figure 3 shows the confinement potential for the device of Fig. 2; the asymmetry of the dot potential is apparent.

Our most interesting results are obtained in the regime where the gate and image potentials are comparable in size. The effective confinement potential, including images, is rather complicated due to the inhomogeneous gate arrangement. Indeed, it is a very poor approximation to neglect screening in this system. The computed wavefunctions for our dot are shown in Fig. 3.

We have also incorporated a SET into our proposed readout device, as shown in Fig. 2. The island of the SET is sheathed by a thin, 2 nm layer of silicon dioxide, and tunnel-coupled to adjoining source and drain gates. A third, capacitively-coupled gate is placed nearby, to provide full control over both the charge and the potential of the island. When the device is operated in the Coulomb blockade regime, it becomes a sensitive electrometer schoelkopf98 . In the gate arrangement of Fig. 2, the SET plays two roles. First, voltage control of the island allows us to vary the size of the dot and thus the energy splitting $`E_{12}`$. Second, capacitive coupling to the quantum dot enables detection of its orbital state. The scheme works as as follows. Since the dot orbitals are spatially distinct, they induce different amounts of charge on the SET island. Consequently, transport currents through the SET will reflect the orbital states of the dot. The device exhibits optimal sensitivity if biased at the half-maximum of the conductance peak. The third SET gate has been introduced specifically to adjust this working point. As expected for our geometry, we find that the SET couples most strongly to the excited electronic orbital.

For the device shown in Fig. 2 we obtain an energy splitting of $`E_{12}=0.129\mathrm{meV}=31.2\mathrm{GHz}`$ between the two lowest orbital states, and the dominant matrix elements $`|1|y|2|=48`$ nm and $`|2|y_y|1|=3.6`$. From these results we obtain the readout oscillation frequency $`\nu _R=5.5\times 10^5\sqrt{I}\mathrm{Hz}`$ (for microwave intensity $`I`$ in units of W/m2), and the orbital decay rate for spontaneous phonon emission $`\mathrm{\Gamma }_I=12.7`$ MHz. The orbital decay rate for emission of a photon is very much less: $`\mathrm{\Gamma }_\gamma =e^2E_{12}^3|1|y|2|^2\sqrt{\epsilon }/3\pi \mathrm{}^4c^3\epsilon _0^{3/2}=6.5`$ Hz, and therefore is not a limiting factor in this scheme. All results are obtained at an ambient temperature of $`T=100`$ mK and a magnetic field of $`H=0.05\mathrm{T}=1.40\mathrm{GHz}`$. The dielectric constant of Si is $`\epsilon =11.9\epsilon _0`$. Note that the natural widths of the states are small compared to their separation. The charge on the SET island was computed by integrating $`𝐧𝐃`$ over the surface. The excess or induced charge of the excited quantum dot orbital, relative to its ground state, was found to be $`\mathrm{\Delta }Q=0.052e`$. The corresponding center-of-charge motion in the dot is $`\mathrm{\Delta }y=4.3`$ nm. However, the charge motion does not track closely with $`\mathrm{\Delta }Q`$; the latter is determined primarily by the capacitive coupling between the dot and the SET. Finally, we find that by changing the bias voltage on the SET by 20% (thus reducing the dot size), the excitation resonance frequency $`\nu _{12}`$ changes by 8 GHz.

The prospects for initialization and readout in this scheme are partially summarized in Fig. 4. Theoretical considerations of shot noise schoelkopf98 ; korotkov99  place an upper bound of about $`4\times 10^6e/\sqrt{\mathrm{Hz}}`$ on the detection sensitivity for charge induced on the island of an optimized rf-SET, as a function of the measurement bandwidth. Similarly, the decay rate of the excited electronic orbital places a lower bound on the readout oscillation frequency $`\nu _R`$. The latter is a function of microwave intensity. By directing $`70`$ pW of microwave power onto a dot of size $`0.1\times 0.1\mu \mathrm{m}^2`$, as consistent with low-temperature transport experiments, we would obtain the $`\nu _R`$ working point marked in Fig. 4, with $`1/\nu _R1/\mathrm{\Gamma }_I`$. Sample heating outside the dot can be minimized by using integrated on-chip antennas to focus the power blick95 . We must also satisfy the requirement $`T_s1/\mathrm{\Gamma }_I=78.7\mathrm{ns}`$ where $`T_s`$ is the spin decoherence time. Theoretical estimates for $`T_1`$ in Si quantum dots exceed 1 ms tahan02 . Experiments on dots have yet to be performed, but spins of electrons in impurity donor states should have similar lifetimes. Recent work on electron spins in 28Si:P lyon03  would imply that $`T_s>\text{60 ms}`$. Hence the crucial hierarchy of timescales $`1/\nu _R1/\mathrm{\Gamma }_I`$ $`T_s`$ should be achievable.

In summary, we have proposed, simulated, and analyzed a scheme for measuring single electron spins in a quantum dot via spin-charge transduction. Adequate charge detection sensitivity is provided by an integrated rf-SET. The scheme also enables fast initialization.

We thank Alex Rimberg and the quantum computing group at the University of Wisconsin-Madison for many beneficial discussions. Our work was supported by the U.S. Army Research Office through ARDA, and the National Science Foundation through the QuBIC program, DMR-0079983, and DMR-0081039.


text:
# Global Stationary Phase and the Sign Problem

## Abstract

We present a computational strategy for reducing the sign problem in the evaluation of high dimensional integrals with non-positive definite weights. The method involves stochastic sampling with a positive semidefinite weight that is adaptively and optimally determined during the course of a simulation. The optimal criterion, which follows from a variational principle for analytic actions $`S(z)`$, is a global stationary phase condition that the average gradient of the phase $`\mathrm{Im}S`$ along the sampling path vanishes. Numerical results are presented from simulations of a model adapted from statistical field theories of classical fluids.

A familiar problem that arises in the context of lattice gauge theorylgt , quantum chemistryquantum1 , correlated electron physicsloh , and equilibrium field theories of classical fluidsghfrev , is the evaluation of integrals of the form

$$Z=_{C_1}𝑑x\mathrm{exp}[S(x)]$$
(1)
where the path of integration $`C_1`$ is the real axis and the action (or effective Hamiltonian) $`S(x)`$ is complex. In the cases of primary interest $`xR^n`$ is a n-vector representing a discrete representation (lattice sites or spectral elements) of one or more classical or quantum fields. The dimension $`n`$ is typically large, of order $`10^310^6`$. Here we shall use one-dimensional notation, although the formalism is primarily intended for cases of $`n1`$.

For real $`S(x)`$, there are a variety of powerful methods available for evaluating $`Z`$, including Monte Carlo and (real) Langevin simulationsbinder . However, in the case of complex $`S=S_R+iS_I`$, the integrand is not positive semidefinite, so the Monte Carlo method is not immediately applicable. Simulations can be carried out using the positive semidefinite weight $`\mathrm{exp}(S_R)`$, but then an oscillatory phase factor of $`\mathrm{exp}(iS_I)`$ must be included in the computation of averageslin . The rapid oscillations in this factor (the “sign problem”), which become more pronounced for large $`n`$, can dramatically slow convergence in such simulations. Alternatively, a “complex Langevin” simulation technique has been devised in which the field variables $`x`$ are extended to the complex plane and a Langevin trajectory prescribed for the purpose of generating Markov chains of statesparisi . Unfortunately this method is not guaranteed to converge and pathological behavior has been noted for specific modelslee ; schoenmaker . In the present letter we describe a new simulation approach that is useful for reducing the sign problem in integrals of the form of Eq. (1), where $`S(z)`$ is an analytic function of the complex n-vector $`z=x+iy`$.

We begin by considering a displacement of the original integration path along the real $`x`$ axis, $`C_1`$ to a new parallel path $`C_y`$ defined by $`z=x+iy,x_j(\mathrm{},\mathrm{})`$, in which $`yR^n`$ is an arbitrary displacement of $`C_1`$ along the imaginary axis. Note that the displacement $`y_j`$ need not be uniform in $`j`$ for the $`n>1`$ case. Provided $`S(z)`$ is analytic in the rectangular strip bounded by $`C_1`$ and $`C_y`$ and $`|\mathrm{exp}[S(R+iy)]|0`$ for $`R\pm \mathrm{}`$, it follows that

$$Z=_{C_y}𝑑z\mathrm{exp}[S(z)]=_{C_y}𝑑x\mathrm{exp}[S(x+iy)]$$
(2)
and the resulting $`Z`$ is independent of the choice of $`y`$. Upon decomposing $`S`$ into real and imaginary parts $`S_R(x,y)+iS_I(x,y)`$, $`Z`$ can be rewritten as

$$Z=Z_y_{C_y}𝑑xP_y(x)\mathrm{exp}[iS_I(x,y)]$$
(3)
where $`Z_y_{C_y}𝑑x\mathrm{exp}[S_R(x,y)]`$ and $`P_y(x)`$ is a normalized, positive semidefinite, probability distribution for a random variable $`x`$ at the fixed value of $`y`$:

$$P_y(x)=\mathrm{exp}[S_R(x,y)]/Z_y$$
(4)
It follows that the average of an analytic observable $`f(x)`$ can be evaluated alternatively from the formulas

$`f(x)`$ $`=`$ $`Z^1{\displaystyle _{C_1}}𝑑x\mathrm{exp}[S(x)]f(x)`$ (5)
$`=`$ $`{\displaystyle \frac{\mathrm{exp}[iS_I(x,y)]f(x+iy)_y}{\mathrm{exp}[iS_I(x,y)]_y}}`$
where $`h(x)_y_{C_y}𝑑xP_y(x)h(x)`$ denotes an average with probability weight $`P_y(x)`$.

It is the second expression in Eq. (5) that is of interest in the present letter. A poor choice of $`y`$ will lead to significant oscillations in the phase factor $`\mathrm{exp}[iS_I(x,y)]`$ as $`x`$ is stochastically varied along the sampling path $`C_y`$ in a simulation. This would drive both numerator and denominator in Eq. (5) to zero and dramatically slow or prevent convergence of average quantities of interest. One approach to alleviate this difficulty would be to choose $`y=y^{}`$, where $`y^{}`$ is the imaginary component of a saddle point $`z^{}`$ defined by $`S^{}(z^{})=0`$. The deformed integration path $`C_y^{}`$ would then be a line passing through the saddle point parallel to the real axis. If this path happened to be a constant phase (steepest ascent) path locally around the saddle point, then the phase oscillations would be reduced on trajectories that remain close to the saddle pointbender . In general, however, path $`C_y^{}`$ will not be a a constant phase path, even in the close vicinity of $`z^{}`$. A local analysis about each saddle point, costing $`O(n^2)`$ in computational effort, can be used to identify proper constant phase paths. However, in typical problems where field fluctuations are strong, significant weight is given to trajectories that are not localized around saddle points.

The essence of our method is a global strategy for selecting an optimal displacement $`y`$, denoted $`\overline{y}`$. To this end, we introduce a “generating” function (functional)

$$G(y)=\mathrm{ln}_{C_y}𝑑x\mathrm{exp}[S_R(x,y)]$$
(6)
Invoking the Cauchy-Riemann (CR) equations, it is straightforward to show that the first derivative of $`G(y)`$ is given by

$$\frac{G(y)}{y_j}=\frac{}{x_j}S_I(x,y)_y$$
(7)
The second derivative follows from repeated application of the CR equations and an integration by parts

$`{\displaystyle \frac{^2G(y)}{y_jy_k}}`$ $`=`$ $`[{\displaystyle \frac{}{x_j}}S_I{\displaystyle \frac{}{x_j}}S_I_y][{\displaystyle \frac{}{x_k}}S_I{\displaystyle \frac{}{x_k}}S_I_y]_y`$ (8)
$`+`$ $`[{\displaystyle \frac{}{x_j}}S_R][{\displaystyle \frac{}{x_k}}S_R]_y,`$
which is the sum of two positive definite forms. It follows that $`G(y)`$ is manifestly a convex function for any $`y`$.

We now claim that the “optimal” choice $`y=\overline{y}`$ is such that

$$\frac{G(y)}{y_j}|_{\overline{y}}=\frac{}{x_j}S_I(x,\overline{y})_{\overline{y}}=0$$
(9)
Evidently such a point would be a local minimum of $`G(y)`$. Moreover, it implies that $`S_I`$ has vanishing gradients on average along the sampling path $`C_{\overline{y}}`$. This condition can be viewed as a global, rather than localbender , stationary phase criterion and would seem to be an excellent way to minimize the effect of phase fluctuations. Since $`G(y)`$ has a unique minimum, it follows that $`\overline{y}_j`$ is homogeneous in $`j`$ for bulk systems with translationally invariant actions. The method evidently produces nontrivial inhomogeneous $`\overline{y}`$ when applied to field theories in bounded geometries.

It remains to discuss how to incorporate this optimal choice of sampling path into a simulation algorithm. We propose the following “optimal path sampling” (OPS) algorithm:

1. Initialize vectors $`x`$ and $`y=y^k`$ with $`k=0`$.
2. Carry out a stochastic simulation in $`x`$ at fixed $`y^k`$ to generate a Markov chain of $`x`$ states of length $`M`$. $`P_{y^k}(x)`$ should be used as a statistical weight for importance sampling. The simulation method could be Metropolis Monte Carlo, its “smart” or “hybrid” variantskennedy , or a real Langevin technique.
3. Evaluate $`G(y^k)`$ and $`G(y^k)/y^k`$ by averaging over the $`x`$ configurations accumulated in the $`M`$-state simulation. Update $`y`$ to approach $`\overline{y}`$ by making a steepest descent step

$$y^{k+1}=y^k\lambda \frac{G(y^k)}{y^k}$$

where $`\lambda `$ is an adjustable relaxation parameter. Alternatively, the accumulated information on $`G(y)`$ could be used to carry out approximate line minimizations, which would permit conjugate gradient updates from $`y^k`$ to $`y^{k+1}`$.
4. Repeat steps 2 and 3 for $`k=1,2,\mathrm{}`$ until the sequence of $`y^k`$ converges to within some prescribed tolerance to $`\overline{y}`$. The simulation has now equilibrated.
5. Carry out a long stochastic simulation (“production run”) with statistical weight $`P_{\overline{y}}(x)`$.
6. Compute averages over the simulated states according to Eq. (5) with $`y=\overline{y}`$.

Evidently, the parameters $`M`$ and $`\lambda `$ can be adjusted to accelerate the “equilibration” period.

Our OPS method has some similarities to (and was inspired by) the complex Langevin (CL) simulation technique. In that approach, one generates a Markov chain of states in the complex plane by integrating the Langevin equationsparisi

$$\frac{x}{t}=\mathrm{Re}\frac{dS}{dz}+\eta (t)$$
(10)
$$\frac{y}{t}=\mathrm{Im}\frac{dS}{dz}$$
(11)
where $`\eta (t)`$ is a real Gaussian white noise with $`\eta (t)=0`$ and $`\eta _j(t)\eta _k(t^{})=2\delta (tt^{})\delta _{jk}`$. Ensemble averages $`f(x)`$ are computed as time averages of $`f(x+iy)`$ over the chain of states. Under conditions where the CL method converges, we have observed that $`y`$ drifts to a nearly constant value that is not associated with any saddle point $`y^{}`$. Eq. (11) reduces approximately in this case to $`\mathrm{Im}dS/dz_y=0`$, which is equivalent to the condition (9). The OPS technique is also distinct from so-called “stationary phase Monte Carlo” methods, which apply filtering and sparse sampling methods to suppress phase oscillationsquantum1 ; sabo . These methods are effective but apparently have no variational basis.

Before providing a numerical example of the OPS method, it is illustrative to see how our global stationary phase criterion works in a simple one dimensional example

$$\mathrm{Ai}(t)=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑x\mathrm{exp}[i(x^3/3+tx)],$$
(12)
which is a representation of the Airy function. In this case $`S(x,t)=i(x^3/3+tx)`$ and Eq. (9) leads to $`\overline{y}^2x^2_{\overline{y}}t=0`$. This equation has a single root, corresponding to the minimum of $`G(y)`$, that yields $`\overline{y}(t)`$. For example, $`\overline{y}(1)=1.19149`$. Of particular interest is the effect of the optimal displacement on phase oscillations. In Fig. 1 we plot $`\mathrm{Re}[\mathrm{exp}(iS_I(x,y,t))]`$ verses $`x`$ at $`t=1`$ for $`y=0`$ (no shift) and $`y=\overline{y}`$ (optimal). Clearly the optimal shift dramatically suppresses phase oscillations over the interval $`2x2`$. The global stationary phase criterion has no effect outside this interval, because $`P_{\overline{y}}(x)`$ decays supra-exponentially there as $`\mathrm{exp}(x^2\overline{y})`$ and so no statistical weight is given to $`|x|2`$.

As a numerical test of the OPS method, we have carried out simulations of the model

$$S(x)=\underset{j=1}{\overset{n}{}}[\alpha x_j^2+(x_{j+1}x_j)^2\chi \mathrm{exp}(ix_j)]$$
(13)
which can be viewed as a lattice field theory for the one-dimensional classical Yukawa fluid in the grand canonical ensemble ($`\alpha `$ is a measure of interaction strength and $`\chi `$ is the activity). For the case of $`n>1`$, periodic boundary conditions are applied. The model has a saddle point $`z_j^{}=iy_j^{}`$ that lies on the imaginary axis and is homogeneous in the index $`j`$ (as well as an inhomogeneous “1d crystal-like” saddle point). Its location is given by the solution of $`\chi \mathrm{exp}(y_j^{})+2\alpha y_j^{}=0`$. The optimal displaced path $`\overline{y}_j`$ is homogeneous in $`j`$ and is given by the solution of $`\chi \mathrm{exp}(\overline{y}_j)\mathrm{cos}x_j_{\overline{y}}+2\alpha \overline{y}_j=0`$. We see that $`y^{}`$ and $`\overline{y}`$ are coincident under conditions ($`\alpha 1`$) where the random variable $`x`$ fluctuates closely about the saddle point $`x^{}=0`$. In the strongly fluctuating regime ($`\alpha 1`$), $`\mathrm{cos}x_j_{\overline{y}}`$ will be dramatically reduced, resulting in a large shift of $`\overline{y}`$ away from $`y^{}`$. These expectations are borne out in numerical simulations of the model.

We have carried out conventional Metropolis Monte Carlo (MC)\[i.e. Eq. (5) with $`y=0`$\], OPS, and CL simulations of the model with action Eq. (13). The results were obtained from runs with a total of $`10^7`$ MC cycles or Langevin steps, a time step of 0.001 in the case of CL, and parameters $`M=1000`$, $`\lambda =0.05`$ for OPS. In Fig. 2 we compare the results obtained from OPS and CL simulations with $`n=1`$ and $`\chi =1`$. The top panel (a) shows $`z`$ as a function of $`\alpha `$, while the bottom panel (b) displays the real and imaginary parts of the “sign” $`\mathrm{exp}(iS_I)`$. In contrast to OPS, CL fails to converge, or converges very slowly, for $`\alpha 0.15`$. Conventional MC also converges, but the average sign is approximately 0.8, as opposed to $`1`$ shown by the OPS.

It is often observedloh  that the sign in conventional MC simulations decreases exponentially with $`n`$, causing a breakdown of the method. This is illustrated for the present model in Fig. 3 with parameters $`n=10`$ and $`\chi =1`$. The conventional MC method fails to converge for $`\alpha 0.1`$ in contrast to OPS. Moreover, the real part of the sign is strongly suppressed in the MC results, even at large values of $`\alpha `$. The sign problem is evidently strongly suppressed, if not eliminated entirely for this model in OPS.

The OPS method is applicable to any field theory with an action $`S(z)`$ that is analytic throughout a domain of $`z`$ relevant to numerical simulations. This includes the important cases of classical fluids in the grand canonical ensemble and path integral formulations of time-dependent quantum chemical problems. Other situations including fluids in the canonical ensemble, strongly correlated electrons, and lattice gauge theories are characterized by analytic $`\mathrm{exp}(S)`$, but with zeros along the real axis and hence logarithmic singularities in $`S`$. We believe that OPS will also be useful in such problems, however precautions should be taken to avoid crossing branch cuts in the steepest descent approach to the optimal displacement $`\overline{y}`$. Finally, we note that the displaced paths considered here were parallel to the real axis. Generalization of the method to optimize both the displacement and shape of the path could prove even more powerful.

In summary, we have identified a variational principle that permits a global stationary phase analysis of integrals of arbitrary dimension with analytic integrands. We expect that this technique will have important implications for analytical and numerical investigations of field theories in the complex plane.

###### Acknowledgements.

This work was supported in part by the NSF under the MRSEC program award No. DMR00-80034 and DMR98-70785. We are grateful to H. Metiu, C. Garcia-Cervera, R. Sugar, J. S. Langer, M. P. A. Fisher, and D. Scalapino for helpful discussions.


text:
# Response to “Comment on ‘time-dependent density-matrix renormalization group: a systematic method for the study of quantum many-body out-of-equilibrium systems’ ” by H. G. Luo, T. Xiang, and X. Q. Wang

The authors of the CommentLXW  propose an improvement to our TdDMRG algorithm as described in our LetterMiguel , providing few details about the proposed extension. Luo, Xiang, and Wang (LXW) then point out that the oscillations at late times ($`t>18`$) shown in two of the lines plotted in Fig. 2 of our manuscript are an artifact of the truncation of the Hilbert space.

LXW are correct that the oscillations we reported in the one case of $`V/w=1.1`$ are in fact spurious. We had systematically examined the errors induced by truncations of the Hilbert space in our Letter. However we failed to notice that, for the case of $`V/w=1.1`$ when the leads are insulating, the eigenvalues of the reduced density matrices become so small so quickly (due to the gap to excitations) that our code discarded the same portion of the Hilbert space for truncations M=256 and 512. Thus in this one case we were left with the illusion that we had achieved convergence, when in fact we had not.

Unfortunately the Comment does not provide enough details regarding the construction of the (reduced) density matrix. As Eq. 1 of the Comment shows, LXW’s density matrix requires knowledge of the wavefunction at different times; but to obtain the time-evolved wavefunction, one needs the density matrix in the first place. Thus it seems that LXW use an iterative procedure: LXW apparently carry out repeated integrations forward in time for each chain size starting with the shortest chain, of length 4 sites. Once the density matrix is constructed for a chain of a given length, the chain is then lengthened by the addition of two sites via the usual (infinite-size) DMRG algorithm, and the time evolution is repeated again from the beginning. The process is then repeated until a sufficiently long chain is built up, tailored to the particular choice of applied bias. Thus the proposed modification of our TdDMRG algorithm is computationally intensive and time consuming. We further point out that the whole procedure must repeated each time a change is made to the bias or other parameters, as the density matrix is tailored to each particular choice of applied bias or parameters. This is in contrast to the original TdDMRG procedure which runs more quickly as it does not require that the chain be built up again each time the bias or the parameters are changed. And while it seems that LXW have found one way to improve the long-time behavior of the TdDMRG algorithm, it is not at all clear that their procedure is optimal. It is unclear, for instance, how much improvement can be expected in the physical case of electrons with spin, as the Hilbert space is so much larger in that case. Other possible improvements to the original TdDMRG algorithm can be imagined.

Finally, we emphasize that our Letter is substantially correct as the Comment corrects only the case $`V/w=1.1`$ shown in Fig. 2 of our Letter. The original TdDMRG scheme itself permits the identification of spurious oscillations at late times. Convergence is tested by checking that observables do not change significantly as $`M`$ is increased. The results of LXW do not reveal new phenomena not already presented in our Letter, and it remains to be demonstrated whether or not the proposed extension is the optimal one.