You are currently browsing the category archive for the ‘math.NA’ category.

Peter Denton, Stephen Parke, Xining Zhang, and I have just uploaded to the arXiv a completely rewritten version of our previous paper, now titled “Eigenvectors from Eigenvalues: a survey of a basic identity in linear algebra“. This paper is now a survey of the various literature surrounding the following basic identity in linear algebra, which we propose to call the eigenvector-eigenvalue identity:

Theorem 1 (Eigenvector-eigenvalue identity) Let ${A}$ be an ${n \times n}$ Hermitian matrix, with eigenvalues ${\lambda_1(A),\dots,\lambda_n(A)}$. Let ${v_i}$ be a unit eigenvector corresponding to the eigenvalue ${\lambda_i(A)}$, and let ${v_{i,j}}$ be the ${j^{th}}$ component of ${v_i}$. Then

$\displaystyle |v_{i,j}|^2 \prod_{k=1; k \neq i}^n (\lambda_i(A) - \lambda_k(A)) = \prod_{k=1}^{n-1} (\lambda_i(A) - \lambda_k(M_j))$

where ${M_j}$ is the ${n-1 \times n-1}$ Hermitian matrix formed by deleting the ${j^{th}}$ row and column from ${A}$.

When we posted the first version of this paper, we were unaware of previous appearances of this identity in the literature; a related identity had been used by Erdos-Schlein-Yau and by myself and Van Vu for applications to random matrix theory, but to our knowledge this specific identity appeared to be new. Even two months after our preprint first appeared on the arXiv in August, we had only learned of one other place in the literature where the identity showed up (by Forrester and Zhang, who also cite an earlier paper of Baryshnikov).

The situation changed rather dramatically with the publication of a popular science article in Quanta on this identity in November, which gave this result significantly more exposure. Within a few weeks we became informed (through private communication, online discussion, and exploration of the citation tree around the references we were alerted to) of over three dozen places where the identity, or some other closely related identity, had previously appeared in the literature, in such areas as numerical linear algebra, various aspects of graph theory (graph reconstruction, chemical graph theory, and walks on graphs), inverse eigenvalue problems, random matrix theory, and neutrino physics. As a consequence, we have decided to completely rewrite our article in order to collate this crowdsourced information, and survey the history of this identity, all the known proofs (we collect seven distinct ways to prove the identity (or generalisations thereof)), and all the applications of it that we are currently aware of. The citation graph of the literature that this ad hoc crowdsourcing effort produced is only very weakly connected, which we found surprising:

The earliest explicit appearance of the eigenvector-eigenvalue identity we are now aware of is in a 1966 paper of Thompson, although this paper is only cited (directly or indirectly) by a fraction of the known literature, and also there is a precursor identity of Löwner from 1934 that can be shown to imply the identity as a limiting case. At the end of the paper we speculate on some possible reasons why this identity only achieved a modest amount of recognition and dissemination prior to the November 2019 Quanta article.

The Polymath15 paper “Effective approximation of heat flow evolution of the Riemann ${\xi}$ function, and a new upper bound for the de Bruijn-Newman constant“, submitted to Research in the Mathematical Sciences, has just been uploaded to the arXiv. This paper records the mix of theoretical and computational work needed to improve the upper bound on the de Bruijn-Newman constant ${\Lambda}$. This constant can be defined as follows. The function

$\displaystyle H_0(z) := \frac{1}{8} \xi\left(\frac{1}{2} + \frac{iz}{2}\right),$

where ${\xi}$ is the Riemann ${\xi}$ function

$\displaystyle \xi(s) := \frac{s(s-1)}{2} \pi^{-s/2} \Gamma\left(\frac{s}{2}\right) \zeta(s)$

has a Fourier representation

$\displaystyle H_0(z) = \int_0^\infty \Phi(u) \cos(zu)\ du$

where ${\Phi}$ is the super-exponentially decaying function

$\displaystyle \Phi(u) := \sum_{n=1}^\infty (2\pi^2 n^4 e^{9u} - 3\pi n^2 e^{5u} ) \exp(-\pi n^2 e^{4u} ).$

The Riemann hypothesis is equivalent to the claim that all the zeroes of ${H_0}$ are real. De Bruijn introduced (in different notation) the deformations

$\displaystyle H_t(z) := \int_0^\infty e^{tu^2} \Phi(u) \cos(zu)\ du$

of ${H_0}$; one can view this as the solution to the backwards heat equation ${\partial_t H_t = -\partial_{zz} H_t}$ starting at ${H_0}$. From the work of de Bruijn and of Newman, it is known that there exists a real number ${\Lambda}$ – the de Bruijn-Newman constant – such that ${H_t}$ has all zeroes real for ${t \geq \Lambda}$ and has at least one non-real zero for ${t < \Lambda}$. In particular, the Riemann hypothesis is equivalent to the assertion ${\Lambda \leq 0}$. Prior to this paper, the best known bounds for this constant were

$\displaystyle 0 \leq \Lambda < 1/2$

with the lower bound due to Rodgers and myself, and the upper bound due to Ki, Kim, and Lee. One of the main results of the paper is to improve the upper bound to

$\displaystyle \Lambda \leq 0.22. \ \ \ \ \ (1)$

At a purely numerical level this gets “closer” to proving the Riemann hypothesis, but the methods of proof take as input a finite numerical verification of the Riemann hypothesis up to some given height ${T}$ (in our paper we take ${T \sim 3 \times 10^{10}}$) and converts this (and some other numerical verification) to an upper bound on ${\Lambda}$ that is of order ${O(1/\log T)}$. As discussed in the final section of the paper, further improvement of the numerical verification of RH would thus lead to modest improvements in the upper bound on ${\Lambda}$, although it does not seem likely that our methods could for instance improve the bound to below ${0.1}$ without an infeasible amount of computation.

We now discuss the methods of proof. An existing result of de Bruijn shows that if all the zeroes of ${H_{t_0}(z)}$ lie in the strip ${\{ x+iy: |y| \leq y_0\}}$, then ${\Lambda \leq t_0 + \frac{1}{2} y_0^2}$; we will verify this hypothesis with ${t_0=y_0=0.2}$, thus giving (1). Using the symmetries and the known zero-free regions, it suffices to show that

$\displaystyle H_{0.2}(x+iy) \neq 0 \ \ \ \ \ (2)$

whenever ${x \geq 0}$ and ${0.2 \leq y \leq 1}$.

For large ${x}$ (specifically, ${x \geq 6 \times 10^{10}}$), we use effective numerical approximation to ${H_t(x+iy)}$ to establish (2), as discussed in a bit more detail below. For smaller values of ${x}$, the existing numerical verification of the Riemann hypothesis (we use the results of Platt) shows that

$\displaystyle H_0(x+iy) \neq 0$

for ${0 \leq x \leq 6 \times 10^{10}}$ and ${0.2 \leq y \leq 1}$. The problem though is that this result only controls ${H_t}$ at time ${t=0}$ rather than the desired time ${t = 0.2}$. To bridge the gap we need to erect a “barrier” that, roughly speaking, verifies that

$\displaystyle H_t(x+iy) \neq 0 \ \ \ \ \ (3)$

for ${0 \leq t \leq 0.2}$, ${x = 6 \times 10^{10} + O(1)}$, and ${0.2 \leq y \leq 1}$; with a little bit of work this barrier shows that zeroes cannot sneak in from the right of the barrier to the left in order to produce counterexamples to (2) for small ${x}$.

To enforce this barrier, and to verify (2) for large ${x}$, we need to approximate ${H_t(x+iy)}$ for positive ${t}$. Our starting point is the Riemann-Siegel formula, which roughly speaking is of the shape

$\displaystyle H_0(x+iy) \approx B_0(x+iy) ( \sum_{n=1}^N \frac{1}{n^{\frac{1+y-ix}{2}}} + \gamma_0(x+iy) \sum_{n=1}^N \frac{n^y}{n^{\frac{1+y+ix}{2}}} )$

where ${N := \sqrt{x/4\pi}}$, ${B_0(x+iy)}$ is an explicit “gamma factor” that decays exponentially in ${x}$, and ${\gamma_0(x+iy)}$ is a ratio of gamma functions that is roughly of size ${(x/4\pi)^{-y/2}}$. Deforming this by the heat flow gives rise to an approximation roughly of the form

$\displaystyle H_t(x+iy) \approx B_t(x+iy) ( \sum_{n=1}^N \frac{b_n^t}{n^{s_*}} + \gamma_t(x+iy) \sum_{n=1}^N \frac{n^y}{n^{\overline{s_*}}} ) \ \ \ \ \ (4)$

where ${B_t(x+iy)}$ and ${\gamma_t(x+iy)}$ are variants of ${B_0(x+iy)}$ and ${\gamma_0(x+iy)}$, ${b_n^t := \exp( \frac{t}{4} \log^2 n )}$, and ${s_*}$ is an exponent which is roughly ${\frac{1+y-ix}{2} + \frac{t}{4} \log \frac{x}{4\pi}}$. In particular, for positive values of ${t}$, ${s_*}$ increases (logarithmically) as ${x}$ increases, and the two sums in the Riemann-Siegel formula become increasingly convergent (even in the face of the slowly increasing coefficients ${b_n^t}$). For very large values of ${x}$ (in the range ${x \geq \exp(C/t)}$ for a large absolute constant ${C}$), the ${n=1}$ terms of both sums dominate, and ${H_t(x+iy)}$ begins to behave in a sinusoidal fashion, with the zeroes “freezing” into an approximate arithmetic progression on the real line much like the zeroes of the sine or cosine functions (we give some asymptotic theorems that formalise this “freezing” effect). This lets one verify (2) for extremely large values of ${x}$ (e.g., ${x \geq 10^{12}}$). For slightly less large values of ${x}$, we first multiply the Riemann-Siegel formula by an “Euler product mollifier” to reduce some of the oscillation in the sum and make the series converge better; we also use a technical variant of the triangle inequality to improve the bounds slightly. These are sufficient to establish (2) for moderately large ${x}$ (say ${x \geq 6 \times 10^{10}}$) with only a modest amount of computational effort (a few seconds after all the optimisations; on my own laptop with very crude code I was able to verify all the computations in a matter of minutes).

The most difficult computational task is the verification of the barrier (3), particularly when ${t}$ is close to zero where the series in (4) converge quite slowly. We first use an Euler product heuristic approximation to ${H_t(x+iy)}$ to decide where to place the barrier in order to make our numerical approximation to ${H_t(x+iy)}$ as large in magnitude as possible (so that we can afford to work with a sparser set of mesh points for the numerical verification). In order to efficiently evaluate the sums in (4) for many different values of ${x+iy}$, we perform a Taylor expansion of the coefficients to factor the sums as combinations of other sums that do not actually depend on ${x}$ and ${y}$ and so can be re-used for multiple choices of ${x+iy}$ after a one-time computation. At the scales we work in, this computation is still quite feasible (a handful of minutes after software and hardware optimisations); if one assumes larger numerical verifications of RH and lowers ${t_0}$ and ${y_0}$ to optimise the value of ${\Lambda}$ accordingly, one could get down to an upper bound of ${\Lambda \leq 0.1}$ assuming an enormous numerical verification of RH (up to height about ${4 \times 10^{21}}$) and a very large distributed computing project to perform the other numerical verifications.

This post can serve as the (presumably final) thread for the Polymath15 project (continuing this post), to handle any remaining discussion topics for that project.

This is the eleventh research thread of the Polymath15 project to upper bound the de Bruijn-Newman constant ${\Lambda}$, continuing this post. Discussion of the project of a non-research nature can continue for now in the existing proposal thread. Progress will be summarised at this Polymath wiki page.

There are currently two strands of activity.  One is writing up the paper describing the combination of theoretical and numerical results needed to obtain the new bound $\Lambda \leq 0.22$.  The latest version of the writeup may be found here, in this directory.  The theoretical side of things have mostly been written up; the main remaining tasks to do right now are

1. giving a more detailed description and illustration of the two major numerical verifications, namely the barrier verification that establishes a zero-free region for $H_t(x+iy)=0$ for $0 \leq t \leq 0.2, 0.2 \leq y \leq 1, |x - 6 \times 10^{10} - 83952| \leq 0.5$, and the Dirichlet series bound that establishes a zero-free region for $t = 0.2, 0.2 \leq y \leq 1, x \geq 6 \times 10^{10} + 83952$; and
2. giving more detail on the conditional results assuming more numerical verification of RH.

Meanwhile, several of us have been exploring the behaviour of the zeroes of $H_t$ for negative $t$; this does not directly lead to any new progress on bounding $\Lambda$ (though there is a good chance that it may simplify the proof of $\Lambda \geq 0$), but there have been some interesting numerical phenomena uncovered, as summarised in this set of slides.  One phenomenon is that for large negative $t$, many of the complex zeroes begin to organise themselves near the curves

$\displaystyle y = -\frac{t}{2} \log \frac{x}{4\pi n(n+1)} - 1.$

(An example of the agreement between the zeroes and these curves may be found here.)  We now have a (heuristic) theoretical explanation for this; we should have an approximation

$\displaystyle H_t(x+iy) \approx B_t(x+iy) \sum_{n=1}^\infty \frac{b_n^t}{n^{s_*}}$

in this region (where $B_t, b_n^t, n^{s_*}$ are defined in equations (11), (15), (17) of the writeup, and the above curves arise from (an approximation of) those locations where two adjacent terms $\frac{b_n^t}{n^{s_*}}$, $\frac{b_{n+1}^t}{(n+1)^{s_*}}$ in this series have equal magnitude (with the other terms being of lower order).

However, we only have a partial explanation at present of the interesting behaviour of the real zeroes at negative t, for instance the surviving zeroes at extremely negative values of $t$ appear to lie on the curve where the quantity $N$ is close to a half-integer, where

$\displaystyle \tilde x := x + \frac{\pi t}{4}$

$\displaystyle N := \sqrt{\frac{\tilde x}{4\pi}}$

The remaining zeroes exhibit a pattern in $(N,u)$ coordinates that is approximately 1-periodic in $N$, where

$\displaystyle u := \frac{4\pi |t|}{\tilde x}.$

A plot of the zeroes in these coordinates (somewhat truncated due to the numerical range) may be found here.

We do not yet have a total explanation of the phenomena seen in this picture.  It appears that we have an approximation

$\displaystyle H_t(x) \approx A_t(x) \sum_{n=1}^\infty \exp( -\frac{|t| \log^2(n/N)}{4(1-\frac{iu}{8\pi})} - \frac{1+i\tilde x}{2} \log(n/N) )$

where $A_t(x)$ is the non-zero multiplier

$\displaystyle A_t(x) := e^{\pi^2 t/64} M_0(\frac{1+i\tilde x}{2}) N^{-\frac{1+i\tilde x}{2}} \sqrt{\frac{\pi}{1-\frac{iu}{8\pi}}}$

and

$\displaystyle M_0(s) := \frac{1}{8}\frac{s(s-1)}{2}\pi^{-s/2} \sqrt{2\pi} \exp( (\frac{s}{2}-\frac{1}{2}) \log \frac{s}{2} - \frac{s}{2} )$

The derivation of this formula may be found in this wiki page.  However our initial attempts to simplify the above approximation further have proven to be somewhat inaccurate numerically (in particular giving an incorrect prediction for the location of zeroes, as seen in this picture).  We are in the process of using numerics to try to resolve the discrepancies (see this page for some code and discussion).

This is the tenth “research” thread of the Polymath15 project to upper bound the de Bruijn-Newman constant ${\Lambda}$, continuing this post. Discussion of the project of a non-research nature can continue for now in the existing proposal thread. Progress will be summarised at this Polymath wiki page.

Most of the progress since the last thread has been on the numerical side, in which the various techniques to numerically establish zero-free regions to the equation $H_t(x+iy)=0$ have been streamlined, made faster, and extended to larger heights than were previously possible.  The best bound for $\Lambda$ now depends on the height to which one is willing to assume the Riemann hypothesis.  Using the conservative verification up to height (slightly larger than) $3 \times 10^{10}$, which has been confirmed by independent work of Platt et al. and Gourdon-Demichel, the best bound remains at $\Lambda \leq 0.22$.  Using the verification up to height $2.5 \times 10^{12}$ claimed by Gourdon-Demichel, this improves slightly to $\Lambda \leq 0.19$, and if one assumes the Riemann hypothesis up to height $5 \times 10^{19}$ the bound improves to $\Lambda \leq 0.11$, contingent on a numerical computation that is still underway.   (See the table below the fold for more data of this form.)  This is broadly consistent with the expectation that the bound on $\Lambda$ should be inversely proportional to the logarithm of the height at which the Riemann hypothesis is verified.

As progress seems to have stabilised, it may be time to transition to the writing phase of the Polymath15 project.  (There are still some interesting research questions to pursue, such as numerically investigating the zeroes of $H_t$ for negative values of $t$, but the writeup does not necessarily have to contain every single direction pursued in the project. If enough additional interesting findings are unearthed then one could always consider writing a second paper, for instance.

Below the fold is the detailed progress report on the numerics by Rudolph Dwars and Kalpesh Muchhal.

This is the seventh “research” thread of the Polymath15 project to upper bound the de Bruijn-Newman constant ${\Lambda}$, continuing this post. Discussion of the project of a non-research nature can continue for now in the existing proposal thread. Progress will be summarised at this Polymath wiki page.

The most recent news is that we appear to have completed the verification that ${H_t(x+iy)}$ is free of zeroes when ${t=0.4}$ and ${y \geq 0.4}$, which implies that ${\Lambda \leq 0.48}$. For very large ${x}$ (for instance when the quantity ${N := \lfloor \sqrt{\frac{x}{4\pi} + \frac{t}{16}} \rfloor}$ is at least ${300}$) this can be done analytically; for medium values of ${x}$ (say when ${N}$ is between ${11}$ and ${300}$) this can be done by numerically evaluating a fast approximation ${A^{eff} + B^{eff}}$ to ${H_t}$ and using the argument principle in a rectangle; and most recently it appears that we can also handle small values of ${x}$, in part due to some new, and significantly faster, numerical ways to evaluate ${H_t}$ in this range.

One obvious thing to do now is to experiment with lowering the parameters ${t}$ and ${y}$ and see what happens. However there are two other potential ways to bound ${\Lambda}$ which may also be numerically feasible. One approach is based on trying to exclude zeroes of ${H_t(x+iy)=0}$ in a region of the form ${0 \leq t \leq t_0}$, ${X \leq x \leq X+1}$ and ${y \geq y_0}$ for some moderately large ${X}$ (this acts as a “barrier” to prevent zeroes from flowing into the region ${\{ 0 \leq x \leq X, y \geq y_0 \}}$ at time ${t_0}$, assuming that they were not already there at time ${0}$). This require significantly less numerical verification in the ${x}$ aspect, but more numerical verification in the ${t}$ aspect, so it is not yet clear whether this is a net win.

Another, rather different approach, is to study the evolution of statistics such as ${S(t) = \sum_{H_t(x+iy)=0: x,y>0} y e^{-x/X}}$ over time. One has fairly good control on such quantities at time zero, and their time derivative looks somewhat manageable, so one may be able to still have good control on this quantity at later times ${t_0>0}$. However for this approach to work, one needs an effective version of the Riemann-von Mangoldt formula for ${H_t}$, which at present is only available asymptotically (or at time ${t=0}$). This approach may be able to avoid almost all numerical computation, except for numerical verification of the Riemann hypothesis, for which we can appeal to existing literature.

Participants are also welcome to add any further summaries of the situation in the comments below.

This is the sixth “research” thread of the Polymath15 project to upper bound the de Bruijn-Newman constant ${\Lambda}$, continuing this post. Discussion of the project of a non-research nature can continue for now in the existing proposal thread. Progress will be summarised at this Polymath wiki page.

The last two threads have been focused primarily on the test problem of showing that ${H_t(x+iy) \neq 0}$ whenever ${t = y = 0.4}$. We have been able to prove this for most regimes of ${x}$, or equivalently for most regimes of the natural number parameter ${N := \lfloor \sqrt{\frac{x}{4\pi} + \frac{t}{16}} \rfloor}$. In many of these regimes, a certain explicit approximation ${A^{eff}+B^{eff}}$ to ${H_t}$ was used, together with a non-zero normalising factor ${B^{eff}_0}$; see the wiki for definitions. The explicit upper bound

$\displaystyle |H_t - A^{eff} - B^{eff}| \leq E_1 + E_2 + E_3$

has been proven for certain explicit expressions ${E_1, E_2, E_3}$ (see here) depending on ${x}$. In particular, if ${x}$ satisfies the inequality

$\displaystyle |\frac{A^{eff}+B^{eff}}{B^{eff}_0}| > \frac{E_1}{|B^{eff}_0|} + \frac{E_2}{|B^{eff}_0|} + \frac{E_3}{|B^{eff}_0|}$

then ${H_t(x+iy)}$ is non-vanishing thanks to the triangle inequality. (In principle we have an even more accurate approximation ${A^{eff}+B^{eff}-C^{eff}}$ available, but it is looking like we will not need it for this test problem at least.)

We have explicit upper bounds on ${\frac{E_1}{|B^{eff}_0|}}$, ${\frac{E_2}{|B^{eff}_0|}}$, ${\frac{E_3}{|B^{eff}_0|}}$; see this wiki page for details. They are tabulated in the range ${3 \leq N \leq 2000}$ here. For ${N \geq 2000}$, the upper bound ${\frac{E_3^*}{|B^{eff}_0|}}$ for ${\frac{E_3}{|B^{eff}_0|}}$ is monotone decreasing, and is in particular bounded by ${1.53 \times 10^{-5}}$, while ${\frac{E_2}{|B^{eff}_0|}}$ and ${\frac{E_1}{|B^{eff}_0|}}$ are known to be bounded by ${2.9 \times 10^{-7}}$ and ${2.8 \times 10^{-8}}$ respectively (see here).

Meanwhile, the quantity ${|\frac{A^{eff}+B^{eff}}{B^{eff}_0}|}$ can be lower bounded by

$\displaystyle |\sum_{n=1}^N \frac{b_n}{n^s}| - |\sum_{n=1}^N \frac{a_n}{n^s}|$

for certain explicit coefficients ${a_n,b_n}$ and an explicit complex number ${s = \sigma + i\tau}$. Using the triangle inequality to lower bound this by

$\displaystyle |b_1| - \sum_{n=2}^N \frac{|b_n|}{n^\sigma} - \sum_{n=1}^N \frac{|a_n|}{n^\sigma}$

we can obtain a lower bound of ${0.18}$ for ${N \geq 2000}$, which settles the test problem in this regime. One can get more efficient lower bounds by multiplying both Dirichlet series by a suitable Euler product mollifier; we have found ${\prod_{p \leq P} (1 - \frac{b_p}{p^s})}$ for ${P=2,3,5,7}$ to be good choices to get a variety of further lower bounds depending only on ${N}$, see this table and this wiki page. Comparing this against our tabulated upper bounds for the error terms we can handle the range ${300 \leq N \leq 2000}$.

In the range ${11 \leq N \leq 300}$, we have been able to obtain a suitable lower bound ${|\frac{A^{eff}+B^{eff}}{B^{eff}_0}| \geq c}$ (where ${c}$ exceeds the upper bound for ${\frac{E_1}{|B^{eff}_0|} + \frac{E_2}{|B^{eff}_0|} + \frac{E_3}{|B^{eff}_0|}}$) by numerically evaluating ${|\frac{A^{eff}+B^{eff}}{B^{eff}_0}|}$ at a mesh of points for each choice of ${N}$, with the mesh spacing being adaptive and determined by ${c}$ and an upper bound for the derivative of ${|\frac{A^{eff}+B^{eff}}{B^{eff}_0}|}$; the data is available here.

This leaves the final range ${N \leq 10}$ (roughly corresponding to ${x \leq 1600}$). Here we can numerically evaluate ${H_t(x+iy)}$ to high accuracy at a fine mesh (see the data here), but to fill in the mesh we need good upper bounds on ${H'_t(x+iy)}$. It seems that we can get reasonable estimates using some contour shifting from the original definition of ${H_t}$ (see here). We are close to finishing off this remaining region and thus solving the toy problem.

Beyond this, we need to figure out how to show that ${H_t(x+iy) \neq 0}$ for ${y > 0.4}$ as well. General theory lets one do this for ${y \geq \sqrt{1-2t} = 0.447\dots}$, leaving the region ${0.4 < y < 0.448}$. The analytic theory that handles ${N \geq 2000}$ and ${300 \leq N \leq 2000}$ should also handle this region; for ${N \leq 300}$ presumably the argument principle will become relevant.

The full argument also needs to be streamlined and organised; right now it sprawls over many wiki pages and github code files. (A very preliminary writeup attempt has begun here). We should also see if there is much hope of extending the methods to push much beyond the bound of ${\Lambda \leq 0.48}$ that we would get from the above calculations. This would also be a good time to start discussing whether to move to the writing phase of the project, or whether there are still fruitful research directions for the project to explore.

Participants are also welcome to add any further summaries of the situation in the comments below.

This is the fifth “research” thread of the Polymath15 project to upper bound the de Bruijn-Newman constant ${\Lambda}$, continuing this post. Discussion of the project of a non-research nature can continue for now in the existing proposal thread. Progress will be summarised at this Polymath wiki page.

We have almost finished off the test problem of showing that ${H_t(x+iy) \neq 0}$ whenever ${t = y = 0.4}$. We have two useful approximations for ${H_t}$, which we have denoted ${A^{eff}+B^{eff}}$ and ${A^{eff}+B^{eff}-C^{eff}}$, and a normalising quantity ${B^{eff}_0}$ that is asymptotically equal to the above expressions; see the wiki page for definitions. In practice, the ${A^{eff}+B^{eff}}$ approximation seems to be accurate within about one or two significant figures, whilst the ${A^{eff}+B^{eff}-C^{eff}}$ approximation is accurate to about three or four. We have an effective upper bound

$\displaystyle |H_t - A^{eff} - B^{eff}| \leq E_1 + E_2 + E_3^*$

where the expressions ${E_1,E_2,E_3^*}$ are quite small in practice (${E_3^*}$ is typically about two orders of magnitude smaller than the main term ${B^{eff}_0}$ once ${x}$ is moderately large, and the error terms ${E_1,E_2}$ are even smaller). See this page for details. In principle we could also obtain an effective upper bound for ${|H_t - (A^{eff} + B^{eff} - C^{eff})|}$ (the ${E_3^*}$ term would be replaced by something smaller).

The ratio ${\frac{A^{eff}+B^{eff}}{B^{eff}_0}}$ takes the form of a difference ${\sum_{n=1}^N \frac{b_n}{n^s} - e^{i\theta} \sum_{n=1}^N \frac{a_n}{n^s}}$ of two Dirichlet series, where ${e^{i\theta}}$ is a phase whose value is explicit but perhaps not terribly important, and the coefficients ${b_n, a_n}$ are explicit and relatively simple (${b_n}$ is ${\exp( \frac{t}{4} \log^2 n)}$, and ${a_n}$ is approximately ${(n/N)^y b_n}$). To bound this away from zero, we have found it advantageous to mollify this difference by multiplying by an Euler product ${\prod_{p \leq P} (1 - \frac{b_p}{p^s})}$ to cancel much of the initial oscillation; also one can take advantage of the fact that the ${b_n}$ are real and the ${a_n}$ are (approximately) real. See this page for details. The upshot is that we seem to be getting good lower bounds for the size of this difference of Dirichlet series starting from about ${x \geq 5 \times 10^5}$ or so. The error terms ${E_1,E_2,E_3^*}$ are already quite small by this stage, so we should soon be able to rigorously keep ${H_t}$ from vanishing at this point. We also have a scheme for lower bounding the difference of Dirichlet series below this range, though it is not clear at present how far we can continue this before the error terms ${E_1,E_2,E_3^*}$ become unmanageable. For very small ${x}$ we may have to explore some faster ways to compute the expression ${H_t}$, which is still difficult to compute directly with high accuracy. One will also need to bound the somewhat unwieldy expressions ${E_1,E_2}$ by something more manageable. For instance, right now these quantities depend on the continuous variable ${x}$; it would be preferable to have a quantity that depends only on the parameter ${N = \lfloor \sqrt{ \frac{x}{4\pi} + \frac{t}{16} }\rfloor}$, as this could be computed numerically for all ${x}$ in the remaining range of interest quite quickly.

As before, any other mathematical discussion related to the project is also welcome here, for instance any summaries of previous discussion that was not covered in this post.

This is the fourth “research” thread of the Polymath15 project to upper bound the de Bruijn-Newman constant ${\Lambda}$, continuing https://terrytao.wordpress.com/2018/01/24/polymath-proposal-upper-bounding-the-de-bruijn-newman-constant/. Progress will be summarised at this Polymath wiki page.

We are getting closer to finishing off the following test problem: can one show that ${H_t(x+iy) \neq 0}$ whenever ${t = y = 0.4}$, ${x \geq 0}$? This would morally show that ${\Lambda \leq 0.48}$. A wiki page for this problem has now been created here. We have obtained a number of approximations ${A+B, A'+B', A^{eff}+B^{eff}, A^{toy}+B^{toy}}$ to ${H_t}$ (see wiki page), though numeric evidence indicates that the approximations are all very close to each other. (Many of these approximations come with a correction term ${C}$, but thus far it seems that we may be able to avoid having to use this refinement to the approximations.) The effective approximation ${A^{eff} + B^{eff}}$ also comes with an effective error bound

$\displaystyle |H_t - A^{eff} - B^{eff}| \leq E_1 + E_2 + E_3$

for some explicit (but somewhat messy) error terms ${E_1,E_2,E_3}$: see this wiki page for details. The original approximations ${A+B, A'+B'}$ can be considered deprecated at this point in favour of the (slightly more complicated) approximation ${A^{eff}+B^{eff}}$; the approximation ${A^{toy}+B^{toy}}$ is a simplified version of ${A^{eff}+B^{eff}}$ which is not quite as accurate but might be useful for testing purposes.

It is convenient to normalise everything by an explicit non-zero factor ${B^{eff}_0}$. Asymptotically, ${(A^{eff} + B^{eff}) / B^{eff}_0}$ converges to 1; numerically, it appears that its magnitude (and also its real part) stays roughly between 0.4 and 3 in the range ${10^5 \leq x \leq 10^6}$, and we seem to be able to keep it (or at least the toy counterpart ${(A^{toy} + B^{toy}) / B^{toy}_0}$) away from zero starting from about ${x \geq 4 \times 10^6}$ (here it seems that there is a useful trick of multiplying by Euler-type factors like ${1 - \frac{1}{2^{1-s}}}$ to cancel off some of the oscillation). Also, the bounds on the error ${(H_t - A^{eff} - B^{eff}) / B^{eff}_0}$ seem to be of size about 0.1 or better in these ranges also. So we seem to be on track to be able to rigorously eliminate zeroes starting from about ${x \geq 10^5}$ or so. We have not discussed too much what to do with the small values of ${x}$; at some point our effective error bounds will become unusable, and we may have to find some more faster ways to compute ${H_t}$.

In addition to this main direction of inquiry, there have been additional discussions on the dynamics of zeroes, and some numerical investigations of the behaviour Lehmer pairs under heat flow. Contributors are welcome to summarise any findings from these discussions from previous threads (or on any other related topic, e.g. improvements in the code) in the comments below.

This is the third “research” thread of the Polymath15 project to upper bound the de Bruijn-Newman constant ${\Lambda}$, continuing this previous thread. Discussion of the project of a non-research nature can continue for now in the existing proposal thread. Progress will be summarised at this Polymath wiki page.

We are making progress on the following test problem: can one show that ${H_t(x+iy) \neq 0}$ whenever ${t = 0.4}$, ${x \geq 0}$, and ${y \geq 0.4}$? This would imply that

$\displaystyle \Lambda \leq 0.4 + \frac{1}{2} (0.4)^2 = 0.48$

which would be the first quantitative improvement over the de Bruijn bound of ${\Lambda \leq 1/2}$ (or the Ki-Kim-Lee refinement of ${\Lambda < 1/2}$). Of course we can try to lower the two parameters of ${0.4}$ later on in the project, but this seems as good a place to start as any. One could also potentially try to use finer analysis of dynamics of zeroes to improve the bound ${\Lambda \leq 0.48}$ further, but this seems to be a less urgent task.

Probably the hardest case is ${y=0.4}$, as there is a good chance that one can then recover the ${y>0.4}$ case by a suitable use of the argument principle. Here we appear to have a workable Riemann-Siegel type formula that gives a tractable approximation for ${H_t}$. To describe this formula, first note that in the ${t=0}$ case we have

$\displaystyle H_0(z) = \frac{1}{8} \xi( \frac{1+iz}{2} )$

and the Riemann-Siegel formula gives

$\displaystyle \xi(s) = \frac{s(s-1)}{2} \pi^{-s/2} \Gamma(s/2) \sum_{n=1}^N \frac{1}{n^s}$

$\displaystyle + \frac{s(s-1)}{2} \pi^{-(1-s)/2} \Gamma((1-s)/2) \sum_{m=1}^M \frac{1}{m^{1-s}}$

$\displaystyle + \frac{s(s-1)}{2} \pi^{-s/2} \Gamma(s/2) \frac{e^{-i\pi s} \Gamma(1-s)}{2\pi i} \int_{C_M} \frac{w^{s-1} e^{-Nw}}{e^w-1}\ dw$

for any natural numbers ${N,M}$, where ${C_M}$ is a contour from ${+\infty}$ to ${+\infty}$ that winds once anticlockwise around the zeroes ${e^{2\pi im}, |m| \leq M}$ of ${e^w-1}$ but does not wind around any other zeroes. A good choice of ${N,M}$ to use here is

$\displaystyle N=M=\lfloor \sqrt{\mathrm{Im}(s)/2\pi}\rfloor = \lfloor \sqrt{\mathrm{Re}(z)/4\pi} \rfloor. \ \ \ \ \ (1)$

In this case, a classical steepest descent computation (see wiki) yields the approximation

$\displaystyle \int_{C_M} \frac{w^{s-1} e^{-Nw}}{e^w-1}\ dw \approx - (2\pi i M)^{s-1} \Psi( \frac{s}{2\pi i M} - N )$

where

$\displaystyle \Psi(\alpha) := 2\pi \frac{\cos \pi(\frac{1}{2}\alpha^2 - \alpha - \pi/8)}{\cos(\pi \alpha)} \exp( \frac{i\pi}{2} \alpha^2 - \frac{5\pi i}{8} ).$

Thus we have

$\displaystyle H_0(z) \approx A^{(0)} + B^{(0)} - C^{(0)}$

where

$\displaystyle A^{(0)} := \frac{1}{8} \frac{s(s-1)}{2} \pi^{-s/2} \Gamma(s/2) \sum_{n=1}^N \frac{1}{n^s}$

$\displaystyle B^{(0)} := \frac{1}{8} \frac{s(s-1)}{2} \pi^{-(1-s)/2} \Gamma((1-s)/2) \sum_{m=1}^M \frac{1}{m^{1-s}}$

$\displaystyle C^{(0)} := \frac{s(s-1)}{2} \pi^{-s/2} \Gamma(s/2) \frac{e^{-i\pi s} \Gamma(1-s)}{2\pi i} (2\pi i M)^{s-1} \Psi( \frac{s}{2\pi i M} - N )$

with ${s := \frac{1+iz}{2}}$ and ${N,M}$ given by (1).

Heuristically, we have derived (see wiki) the more general approximation

$\displaystyle H_t(z) \approx A + B - C$

for ${t>0}$ (and in particular for ${t=0.4}$), where

$\displaystyle A := \frac{1}{8} \frac{s(s-1)}{2} \pi^{-s/2} \Gamma(s/2) \sum_{n=1}^N \frac{\exp(\frac{t}{16} \log^2 \frac{s+4}{2\pi n^2} )}{n^s}$

$\displaystyle B := \frac{1}{8} \frac{s(s-1)}{2} \pi^{-(1-s)/2} \Gamma((1-s)/2) \sum_{m=1}^M \frac{\exp(\frac{t}{16} \log^2 \frac{5-s}{2\pi m^2} )}{m^{1-s}}$

$\displaystyle C := \exp(-\frac{t \pi^2}{64}) C^{(0)}.$

In practice it seems that the ${C}$ term is negligible once the real part ${x}$ of ${z}$ is moderately large, so one also has the approximation

$\displaystyle H_t(z) \approx A + B.$

For large ${x}$, and for fixed ${t,y>0}$, e.g. ${t=y=0.4}$, the sums ${A,B}$ converge fairly quickly (in fact the situation seems to be significantly better here than the much more intensively studied ${t=0}$ case), and we expect the first term

$\displaystyle B_0 := \frac{1}{8} \frac{s(s-1)}{2} \pi^{-(1-s)/2} \Gamma((1-s)/2) \exp( \frac{t}{16} \log^2 \frac{5-s}{2\pi} )$

of the ${B}$ series to dominate. Indeed, analytically we know that ${\frac{A+B-C}{B_0} \rightarrow 1}$ (or ${\frac{A+B}{B_0} \rightarrow 1}$) as ${x \rightarrow \infty}$ (holding ${y}$ fixed), and it should also be provable that ${\frac{H_t}{B_0} \rightarrow 1}$ as well. Numerically with ${t=y=0.4}$, it seems in fact that ${\frac{A+B-C}{B_0}}$ (or ${\frac{A+B}{B_0}}$) stay within a distance of about ${1/2}$ of ${1}$ once ${x}$ is moderately large (e.g. ${x \geq 2 \times 10^5}$). This raises the hope that one can solve the toy problem of showing ${H_t(x+iy) \neq 0}$ for ${t=y=0.4}$ by numerically controlling ${H_t(x+iy) / B_0}$ for small ${x}$ (e.g. ${x \leq 2 \times 10^5}$), numerically controlling ${(A+B)/B_0}$ and analytically bounding the error ${(H_t - A - B)/B_0}$ for medium ${x}$ (e.g. ${2 \times 10^5 \leq x \leq 10^7}$), and analytically bounding both ${(A+B)/B_0}$ and ${(H_t-A-B)/B_0}$ for large ${x}$ (e.g. ${x \geq 10^7}$). (These numbers ${2 \times 10^5}$ and ${10^7}$ are arbitrarily chosen here and may end up being optimised to something else as the computations become clearer.)

Thus, we now have four largely independent tasks (for suitable ranges of “small”, “medium”, and “large” ${x}$):

1. Numerically computing ${H_t(x+iy) / B_0}$ for small ${x}$ (with enough accuracy to verify that there are no zeroes)
2. Numerically computing ${(A+B)/B_0}$ for medium ${x}$ (with enough accuracy to keep it away from zero)
3. Analytically bounding ${(A+B)/B_0}$ for large ${x}$ (with enough accuracy to keep it away from zero); and
4. Analytically bounding ${(H_t - A - B)/B_0}$ for medium and large ${x}$ (with a bound that is better than the bound away from zero in the previous two tasks).

Note that tasks 2 and 3 do not directly require any further understanding of the function ${H_t}$.

Below we will give a progress report on the numeric and analytic sides of these tasks.

— 1. Numerics report (contributed by Sujit Nair) —

There is some progress on the code side but not at the pace I was hoping. Here are a few things which happened (rather, mistakes which were taken care of).

1. We got rid of code which wasn’t being used. For example, @dhjpolymath computed ${H_t}$ based on an old version but only realized it after the fact.
2. We implemented tests to catch human/numerical bugs before a computation starts. Again, we lost some numerical cycles but moving forward these can be avoided.
3. David got set up on GitHub and he is able to compare his output (in C) with the Python code. That is helping a lot.

Two areas which were worked on were

1. Computing ${H_t}$ and zeroes for ${t}$ around ${0.4}$
2. Computing quantities like ${(A+B-C)/B_0}$, ${(A+B)/B_0}$, ${C/B_0}$, etc. with the goal of understanding the zero free regions.

Some observations for ${t=0.4}$, ${y=0.4}$, ${x \in ( 10^4, 10^7)}$ include:

• ${(A+B) / B_0}$ does seem to avoid the negative real axis
• ${|(A+B) / B0| > 0.4}$ (based on the oscillations and trends in the plots)
• ${|C/B_0|}$ seems to be settling around ${10^{-4}}$ range.

See the figure below. The top plot is on the complex plane and the bottom plot is the absolute value. The code to play with this is here.

— 2. Analysis report —

The Riemann-Siegel formula and some manipulations (see wiki) give ${H_0 = A^{(0)} + B^{(0)} - \tilde C^{(0)}}$, where

$\displaystyle A^{(0)} = \frac{2}{8} \sum_{n=1}^N \int_C \exp( \frac{s+4}{2} u - e^u - \frac{s}{2} \log(\pi n^2) )\ du$

$\displaystyle - \frac{3}{8} \sum_{n=1}^N \int_C \exp( \frac{s+2}{2} u - e^u - \frac{s}{2} \log(\pi n^2) )\ du$

$\displaystyle B^{(0)} = \frac{2}{8} \sum_{m=1}^M \int_{\overline{C}} \exp( \frac{5-s}{2} u - e^u - \frac{1-s}{2} \log(\pi m^2) )\ du$

$\displaystyle - \frac{3}{8} \sum_{m=1}^M \int_C \exp( \frac{3-s}{2} u - e^u - \frac{1-s}{2} \log(\pi m^2) )\ du$

$\displaystyle \tilde C^{(0)} := -\frac{2}{8} \sum_{n=0}^\infty \frac{e^{-i\pi s/2} e^{i\pi s n}}{2^s \pi^{1/2}} \int_{\overline{C}} \int_{C_M} \frac{w^{s-1} e^{-Nw}}{e^w-1} \exp( \frac{5-s}{2} u - e^u)\ du dw$

$\displaystyle +\frac{3}{8} \sum_{n=0}^\infty \frac{e^{-i\pi s/2} e^{i\pi s n}}{2^s \pi^{1/2}} \int_{\overline{C}} \int_{C_M} \frac{w^{s-1} e^{-Nw}}{e^w-1} \exp( \frac{3-s}{2} u - e^u)\ du dw$

where ${C}$ is a contour that goes from ${+i\infty}$ to ${+\infty}$ staying a bounded distance away from the upper imaginary and right real axes, and ${\overline{C}}$ is the complex conjugate of ${C}$. (In each of these sums, it is the first term that should dominate, with the second one being about ${O(1/x)}$ as large.) One can then evolve by the heat flow to obtain ${H_t = \tilde A + \tilde B - \tilde C}$, where

$\displaystyle \tilde A := \frac{2}{8} \sum_{n=1}^N \int_C \exp( \frac{s+4}{2} u - e^u - \frac{s}{2} \log(\pi n^2) + \frac{t}{16} (u - \log(\pi n^2))^2)\ du$

$\displaystyle - \frac{3}{8} \sum_{n=1}^N \int_C \exp( \frac{s+2}{2} u - e^u - \frac{s}{2} \log(\pi n^2) + \frac{t}{16} (u - \log(\pi n^2))^2)\ du$

$\displaystyle \tilde B := \frac{2}{8} \sum_{m=1}^M \int_{\overline{C}} \exp( \frac{5-s}{2} u - e^u - \frac{1-s}{2} \log(\pi m^2) + \frac{t}{16} (u - \log(\pi m^2))^2)\ du$

$\displaystyle - \frac{3}{8} \sum_{m=1}^M \int_C \exp( \frac{3-s}{2} u - e^u - \frac{1-s}{2} \log(\pi m^2) + \frac{t}{16} (u - \log(\pi m^2))^2)\ du$

$\displaystyle \tilde C := -\frac{2}{8} \sum_{n=0}^\infty \frac{e^{-i\pi s/2} e^{i\pi s n}}{2^s \pi^{1/2}} \int_{\overline{C}} \int_{C_M}$

$\displaystyle \frac{w^{s-1} e^{-Nw}}{e^w-1} \exp( \frac{5-s}{2} u - e^u + \frac{t}{4} (i \pi(n-1/2) + \log \frac{w}{2\sqrt{\pi}} - \frac{u}{2})^2) \ du dw$

$\displaystyle +\frac{3}{8} \sum_{n=0}^\infty \frac{e^{-i\pi s/2} e^{i\pi s n}}{2^s \pi^{1/2}} \int_{\overline{C}} \int_{C_M}$

$\displaystyle \frac{w^{s-1} e^{-Nw}}{e^w-1} \exp( \frac{3-s}{2} u - e^u + \frac{t}{4} (i \pi(n-1/2) + \log \frac{w}{2\sqrt{\pi}} - \frac{u}{2})^2)\ du dw.$

Steepest descent heuristics then predict that ${\tilde A \approx A}$, ${\tilde B \approx B}$, and ${\tilde C \approx C}$. For the purposes of this project, we will need effective error estimates here, with explicit error terms.

A start has been made towards this goal at this wiki page. Firstly there is a “effective Laplace method” lemma that gives effective bounds on integrals of the form ${\int_I e^{\phi(x)} \psi(x)\ dx}$ if the real part of ${\phi(x)}$ is either monotone with large derivative, or has a critical point and is decreasing on both sides of that critical point. In principle, all one has to do is manipulate expressions such as ${\tilde A - A}$, ${\tilde B - B}$, ${\tilde C - C}$ by change of variables, contour shifting and integration by parts until it is of the form to which the above lemma can be profitably applied. As one may imagine though the computations are messy, particularly for the ${\tilde C}$ term. As a warm up, I have begun by trying to estimate integrals of the form

$\displaystyle \int_C \exp( s (1+u-e^u) + \frac{t}{16} (u+b)^2 )\ du$

for smallish complex numbers ${b}$, as these sorts of integrals appear in the form of ${\tilde A, \tilde B, \tilde C}$. As of this time of writing, there are effective bounds for the ${b=0}$ case, and I am currently working on extending them to the ${b \neq 0}$ case, which should give enough control to approximate ${\tilde A - A}$ and ${\tilde B-B}$. The most complicated task will be that of upper bounding ${\tilde C}$, but it also looks eventually doable.

This is the first official “research” thread of the Polymath15 project to upper bound the de Bruijn-Newman constant ${\Lambda}$. Discussion of the project of a non-research nature can continue for now in the existing proposal thread. Progress will be summarised at this Polymath wiki page.

The proposal naturally splits into at least three separate (but loosely related) topics:

• Numerical computation of the entire functions ${H_t(z)}$, with the ultimate aim of establishing zero-free regions of the form ${\{ x+iy: 0 \leq x \leq T, y \geq \varepsilon \}}$ for various ${T, \varepsilon > 0}$.
• Improved understanding of the dynamics of the zeroes ${z_j(t)}$ of ${H_t}$.
• Establishing the zero-free nature of ${H_t(x+iy)}$ when ${y \geq \varepsilon > 0}$ and ${x}$ is sufficiently large depending on ${t}$ and ${\varepsilon}$.

Below the fold, I will present each of these topics in turn, to initiate further discussion in each of them. (I thought about splitting this post into three to have three separate discussions, but given the current volume of comments, I think we should be able to manage for now having all the comments in a single post. If this changes then of course we can split up some of the discussion later.)

To begin with, let me present some formulae for computing ${H_t}$ (inspired by similar computations in the Ki-Kim-Lee paper) which may be useful. The initial definition of ${H_t}$ is

$\displaystyle H_t(z) := \int_0^\infty e^{tu^2} \Phi(u) \cos(zu)\ du$

where

$\displaystyle \Phi(u) := \sum_{n=1}^\infty (2\pi^2 n^4 e^{9u} - 3 \pi n^2 e^{5u}) \exp(- \pi n^2 e^{4u} )$

is a variant of the Jacobi theta function. We observe that ${\Phi}$ in fact extends analytically to the strip

$\displaystyle \{ u \in {\bf C}: -\frac{\pi}{8} < \mathrm{Im} u < \frac{\pi}{8} \}, \ \ \ \ \ (1)$

as ${e^{4u}}$ has positive real part on this strip. One can use the Poisson summation formula to verify that ${\Phi}$ is even, ${\Phi(-u) = \Phi(u)}$ (see this previous post for details). This lets us obtain a number of other formulae for ${H_t}$. Most obviously, one can unfold the integral to obtain

$\displaystyle H_t(z) = \frac{1}{2} \int_{-\infty}^\infty e^{tu^2} \Phi(u) e^{izu}\ du.$

In my previous paper with Brad, we used this representation, combined with Fubini’s theorem to swap the sum and integral, to obtain a useful series representation for ${H_t}$ in the ${t<0}$ case. Unfortunately this is not possible in the ${t>0}$ case because expressions such as ${e^{tu^2} e^{9u} \exp( -\pi n^2 e^{4u} ) e^{izu}}$ diverge as ${u}$ approaches ${-\infty}$. Nevertheless we can still perform the following contour integration manipulation. Let ${0 \leq \theta < \frac{\pi}{8}}$ be fixed. The function ${\Phi}$ decays super-exponentially fast (much faster than ${e^{tu^2}}$, in particular) as ${\mathrm{Re} u \rightarrow +\infty}$ with ${-\infty \leq \mathrm{Im} u \leq \theta}$; as ${\Phi}$ is even, we also have this decay as ${\mathrm{Re} u \rightarrow -\infty}$ with ${-\infty \leq \mathrm{Im} u \leq \theta}$ (this is despite each of the summands in ${\Phi}$ having much slower decay in this direction – there is considerable cancellation!). Hence by the Cauchy integral formula we have

$\displaystyle H_t(z) = \frac{1}{2} \int_{i\theta-\infty}^{i\theta+\infty} e^{tu^2} \Phi(u) e^{izu}\ du.$

Splitting the horizontal line from ${i\theta-\infty}$ to ${i\theta+\infty}$ at ${i\theta}$ and using the even nature of ${\Phi(u)}$, we thus have

$\displaystyle H_t(z) = \frac{1}{2} ( \int_{i\theta}^{i\theta+\infty} e^{tu^2} \Phi(u) e^{izu}\ du + \int_{-i\theta}^{-i\theta+\infty} e^{tu^2} \Phi(u) e^{-izu}\ du.$

Using the functional equation ${\Phi(\overline{u}) = \overline{\Phi(u)}}$, we thus have the representation

$\displaystyle H_t(z) = \frac{1}{2} ( K_{t,\theta}(z) + \overline{K_{t,\theta}(\overline{z})} ) \ \ \ \ \ (2)$

where

$\displaystyle K_{t,\theta}(z) := \int_{i\theta}^{i \theta+\infty} e^{tu^2} \Phi(u) e^{izu}\ du$

$\displaystyle = \sum_{n=1}^\infty 2 \pi^2 n^4 I_{t, \theta}( z - 9i, \pi n^2 ) - 3 \pi n^2 I_{t,\theta}( z - 5i, \pi n^2 )$

where ${I_{t,\theta}(b,\beta)}$ is the oscillatory integral

$\displaystyle I_{t,\theta}(b,\beta) := \int_{i\theta}^{i\theta+\infty} \exp( tu^2 - \beta e^{4u} + i b u )\ du. \ \ \ \ \ (3)$

The formula (2) is valid for any ${0 \leq \theta < \frac{\pi}{8}}$. Naively one would think that it would be simplest to take ${\theta=0}$; however, when ${z=x+iy}$ and ${x}$ is large (with ${y}$ bounded), it seems asymptotically better to take ${\theta}$ closer to ${\pi/8}$, in particular something like ${\theta = \frac{\pi}{8} - \frac{1}{4x}}$ seems to be a reasonably good choice. This is because the integrand in (3) becomes significantly less oscillatory and also much lower in amplitude; the ${\exp(ibu)}$ term in (3) now generates a factor roughly comparable to ${\exp( - \pi x/8 )}$ (which, as we will see below, is the main term in the decay asymptotics for ${H_t(x+iy)}$), while the ${\exp( - \beta e^{4u} )}$ term still exhibits a reasonable amount of decay as ${u \rightarrow \infty}$. We will use the representation (2) in the asymptotic analysis of ${H_t}$ below, but it may also be a useful representation to use for numerical purposes.