Computing the Tracy–Widom Distribution for Arbitrary $\beta>0$

Computing the Tracy–Widom Distribution
for Arbitrary $\boldsymbol{\beta>0}$††This paper is a contribution to the
Special Issue on Evolution Equations, Exactly Solvable Models and Random
Matrices in honor of Alexander Its’ 70th birthday. The full collection is
available at https://www.emis.de/journals/SIGMA/Its.html

Thomas TROGDON and Yiting ZHANG

T. Trogdon and Y. Zhang

Department of Applied Mathematics, University of Washington, Seattle,
Washington, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>
Received April 19, 2023, in final form January 03, 2024; Published online
January 13, 2024

We compute the Tracy–Widom distribution describing the asymptotic distribution
of the largest eigenvalue of a large random matrix by solving a boundary-value
problem posed by Bloemendal in his Ph.D. Thesis (2011). The distribution is
computed in two ways. The first method is a second-order finite-difference
method and the second is a highly accurate Fourier spectral method. Since
$\beta$ is simply a parameter in the boundary-value problem, any $\beta>0$ can
be used, in principle. The limiting distribution of the $n$th largest
eigenvalue can also be computed. Our methods are available in the Julia
package TracyWidomBeta.jl.

numerical differential equation; Tracy–Widom distribution; Fourier
transformation

65M06; 60B20; 60H25

##  1 Introduction

Tracy and Widom [23, 24, 25] introduced the Tracy–Widom distribution that
gives the limiting distribution of the rescaled largest eigenvalue of a random
matrix taken from an appropriate symmetry class. More precisely, the largest
eigenvalue $\lambda_{\max}$ satisfies the following fundamental limit

$\lim_{n\to\infty}\mathbb{P}\big{(}n^{1/6}(\lambda_{\max}(A_{n})-2\sqrt{n})\leq
x\big{)}=F_{\beta}(x)\qquad\text{for all}\ x\in\mathbb{R},$

where $F_{\beta}$ is the Tracy–Widom distribution and $\beta=1$, $2$, $4$ if
$A_{n}\sim$ Gaussian orthogonal ensemble, Gaussian unitary ensemble, Gaussian
symplectic ensemble, respectively.

For an $n\times n$ Gaussian ensemble $A_{n}$ with ordered eigenvalues
$\\{\lambda_{j}\\}_{j=1}^{n}$, the joint probability density function (jpdf)
for its eigenvalues is given by

$\rho\left(\lambda_{1},\lambda_{2},\dots,\lambda_{n}\right)=Z^{-1}_{n,\beta}\prod_{i}{\rm
e}^{-\beta\lambda^{2}_{i}/4}\prod_{i<j}\left|\lambda_{j}-\lambda_{i}\right|^{\beta},\qquad\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n},$
(1.1)

where $Z_{n,\beta}$ is the partition function. For $\beta=1,2,4$, (1.1) is
solvable: all correlation functions in finite dimensions can be explicitly
calculated using Hermite polynomials. This provides definitive local limit
theorems and establishes clear limits for the random points [15]. Importantly,
(1.1) also describes a one-dimensional Coulomb gas at inverse temperature
$\beta$ for any $\beta>0$. If $\beta\neq 1,2,4,6$, there are no known explicit
formulae which appear amenable to asymptotic analysis (see, for example, [10,
13, 20] for $\beta=6$). One, in general, must resort to numerical
computations, see [3] for $\beta=1,2,4$.

Dumitriu and Edelman [7] established that a family of symmetric tridiagonal
(Jacobi) matrix models have (1.1) as their eigenvalue density for any
$\beta>0$. More specifically, the matrix given by

$H_{n}^{\beta}=\frac{1}{\sqrt{\beta}}\begin{bmatrix}g&\chi_{(n-1)\beta}&&&\\\
\chi_{(n-1)\beta}&g&\chi_{(n-2)\beta}&&\\\ &\chi_{(n-2)\beta}&g&\ddots&\\\
&&\ddots&\ddots&\chi_{\beta}\\\ &&&\chi_{\beta}&g\end{bmatrix},$ (1.2)

has (1.1) as the jpdf of its eigenvalues. Here the entries are independent
random variables, up to symmetry, $g\sim N(0,2)$ and $\chi_{k}\sim$ Chi($k$).
We call $H^{\beta}_{n}$ the $\beta$-Hermite ensemble.

Sutton and Edelman [8, 21] then presented an argument111A version of this
argument was first presented by Edelman at the SIAM Conference on Applied
Linear Algebra held at the College of William $\&$ Mary in 2003. describing
how the spectrum of the rescaled operator

$\tilde{H}_{n}^{\beta}=n^{1/6}\big{(}2\sqrt{n}I-H_{n}^{\beta}\big{)},$

where $I$ is the $n\times n$ identity matrix, should be described by the
spectrum of the stochastic Airy operator

$\mathcal{H}_{\beta}=-\frac{{\rm d}^{2}}{{\rm
d}x^{2}}+x+\frac{2}{\sqrt{\beta}}b^{\prime}_{x},$

as $n\to\infty$, where $b^{\prime}_{x}$ is standard Gaussian white noise.

As a result, the following eigenvalue problem was considered in [19]

$\mathcal{H}_{\beta}f=\Lambda f\qquad\text{on}\ L^{2}(\mathbb{R}_{+})$

with a Dirichlet boundary condition $f(0)=0$. Ramírez, Rider, and Virág [19]
proved the following theorem.

######  Theorem 1.1.

With probability one, for each $k\geq 0$, the set of eigenvalues of
$\mathcal{H}_{\beta}$ has a well-defined $(k+1)$st lowest element
$\Lambda_{k}$. Moreover, let $\lambda_{1}\geq\lambda_{2}\geq\cdots$ denote the
eigenvalues of ${H}^{\beta}_{n}$. Then the vector

$\displaystyle\big{(}n^{1/6}(2\sqrt{n}-\lambda_{l})\big{)}_{l=1,\dots,k}$

converges in distribution to $(\Lambda_{0},\Lambda_{1},\dots,\Lambda_{k-1})$
as $n\to\infty$.

Ramírez, Rider, and Virág showed that the distribution of $-\Lambda_{0}$ is a
consistent definition of Tracy–Widom($\cdot$) for general $\beta>0$. Based on
[19, 26], Bloemendal and Viraǵ [1, 2] considered a generalized eigenvalue
problem

$\mathcal{H}_{\beta}f=\Lambda f\qquad\text{on}\ L^{2}(\mathbb{R}_{+})$

with boundary condition $f^{\prime}(0)=\omega f(0)$, where
$\omega\in\mathbb{R}$ represents a scaling parameter. Let
$\mathcal{H}_{\beta,\omega}$ denote $\mathcal{H}_{\beta}$ together with this
boundary condition. According to [19], the distribution $F_{\beta,\omega}$ of
$-\Lambda_{0}$ in the case $\omega=+\infty$ thus coincides with Tracy–Widom
for general $\beta>0$.

Bloemendal and Virág [2] further showed that $F=F_{\beta,\omega}$ can be
characterized using the solution of a boundary value problem, as we now
discuss in the following section. We also point out that Bloemendal [1]
provided Mathematica code to approximate $F$. The goal of the current work is
to expand and improve upon this scheme.

This paper is laid out as follows. In Section 2, we outline our two algorithms
to compute the Tracy–Widom distribution. In Section 3, we validate and compare
our methods. In Section 4, we present a number of additional numerical
results. We also include two appendices to discuss some nuances in the
numerical algorithms (Appendix A) and to discuss the large $\beta$ limit
(Appendix B). Code to produce all figures in this paper can be found here
[27].

##  2 Algorithm description

The Tracy–Widom distribution function $F_{\beta}$ can be characterized as
follows [1]. Consider

$\frac{\partial F}{\partial
x}+\frac{2}{\beta}\frac{\partial^{2}F}{\partial\omega^{2}}+\big{(}x-\omega^{2}\big{)}\frac{\partial
F}{\partial\omega}=0\qquad\text{for}\ (x,\omega)\in\mathbb{R}^{2},$ (2.1)

with boundary conditions given by

$\displaystyle F(x,\omega)\to 1\qquad\text{as}\ x,\omega\to\infty\
\text{together,}$ $\displaystyle F(x,\omega)\to 0\qquad\text{as}\
\omega\to-\infty\ \text{with}x\ \text{bounded above.}$

Then [2, Theorem 1.7]

$F_{\beta}(x)=\lim_{\omega\to\infty}F(x,\omega).$ (2.2)

Using $\omega=-\cot\theta$, we rewrite (2.1) as

$\frac{\partial H}{\partial
x}+\left(\frac{2}{\beta}\sin^{4}\theta\right)\frac{\partial^{2}H}{\partial\theta^{2}}+\left[\left(x+\frac{2}{\beta}\sin
2\theta\right)\sin^{2}\theta-\cos^{2}\theta\right]\frac{\partial
H}{\partial\theta}=0,$ (2.3)

with boundary condition

$\displaystyle H(x,0)=0.$

To address the boundary condition as $x\to\infty$ and $\theta\to\pi$, we
truncate the domain at a finite value, denoted as $x=x_{0}>0$. We then employ
the Gaussian asymptotics established by Bloemendal [1, Theorem 4.1.1] to
obtain the following approximate asymptotic initial condition [1, p. 103]:

$\displaystyle
H(x_{0},\theta)=\begin{cases}\Phi\left(\dfrac{x_{0}-\cot^{2}\theta}{\sqrt{(4/\beta)\cot\theta}}\right),&0\leq\theta\leq\pi/2,\\\
1,&\pi/2\leq\theta.\end{cases}$ (2.4)

Here $\Phi$ denotes the standard normal distribution function. Note that
$H(x_{0},\theta)$ is continuous in both $x_{0}>0$ and $\theta\geq 0$. One then
approximates the undeformed Tracy–Widom distribution $F_{\beta}(x)$ at
$\theta=\pi$ by $\omega=-\cot\theta$ and (2.2), i.e., $F_{\beta}(x)\approx
H(x,\pi)$.

###  2.1 Finite-difference discretization

One way to solve this boundary-value problem is by discretizing it using
finite differences, see, for example, [12]. For $0\leq n\leq N$, $0\leq m\leq
M$, define

$x_{n}=x_{0}+n\Delta x,\qquad\theta_{m}=mh.$

Given that $\beta>0$, one integrates equation (2.3) backward in “time” with
respect to the time-like variable $x$ to guarantee its well-posedness: $\Delta
x<0$. Denote the approximation of $H(x,\theta_{m})$ by $H_{m}(x)$. We then
replace partial derivatives of $H$ with respect to $\theta$ by centered
differences, with grid spacing $h$. The method of lines formulation reads

$\frac{\partial\boldsymbol{H}_{M}}{\partial
x}=-\bigg{(}\frac{2\sin^{4}\theta}{\beta
h^{2}}\bigg{)}\check{T}\boldsymbol{H}_{M}+\frac{\big{[}\big{(}x+\frac{2}{\beta}\sin
2\theta\big{)}\sin^{2}\theta-\cos^{2}\theta\big{]}}{2h}\check{U}\boldsymbol{H}_{M},$
(2.5)

where $\boldsymbol{H}_{M}(x)=[H_{1}(x),H_{2}(x),\dots,H_{M}(x)]^{\textsf{T}}$,
$Mh=\theta_{M}$, and

$\check{T}=\begin{bmatrix}-2&\hphantom{-}1\\\
\hphantom{-}1&-2&\hphantom{-}1\\\ &\hphantom{-}1&-2&1\\\
&&\ddots&\ddots&\ddots\\\ &&&1&-2&\hphantom{-}1\\\
&&&&\hphantom{-}1&-2\end{bmatrix},\qquad\check{U}=\begin{bmatrix}0&-1\\\
1&\hphantom{-}0&-1\\\ &\hphantom{-}1&\hphantom{-}0&-1\\\
&&\ddots&\ddots&\ddots\\\ &&&\hphantom{-}1&0&-1\\\ &&&-1&4&-3\end{bmatrix}.$
(2.6)

$\boldsymbol{H}_{M}$ excludes $H_{0}$ since $H_{0}(x)=0$. The coefficients in
the last row of $\check{U}$ come from the parameters of the two-step backward
difference formula (BDF2). Note that there is no need to replace the last row
of $\check{T}$ using a backward difference formula since $\theta_{M}=\pi$ and,
mathematically, there is no need to set an additional boundary condition since
(2.5) has vanishing diffusivity for $\theta=\pi$.

Let

$T(\beta,\boldsymbol{\theta}_{M},h):=-\bigg{(}\frac{2\sin^{4}\boldsymbol{\theta}_{M}}{\beta
h^{2}}\bigg{)}\check{T}=\begin{bmatrix}\dfrac{-2\sin^{4}\theta_{1}}{\beta
h^{2}}&&\\\ &\ddots&\\\ &&\dfrac{-2\sin^{4}\theta_{M}}{\beta
h^{2}}\end{bmatrix}\check{T},$ (2.7)

and

$\displaystyle
U(\beta,x,\boldsymbol{\theta}_{M},h):=\frac{\big{[}\big{(}x+\frac{2}{\beta}\sin
2\boldsymbol{\theta}_{M}\big{)}\sin^{2}\boldsymbol{\theta}_{M}-\cos^{2}\boldsymbol{\theta}_{M}\big{]}}{2h}\check{U}$
$\displaystyle=\begin{bmatrix}\dfrac{(x+\frac{2}{\beta}\sin
2\theta_{1})\sin^{2}\theta_{1}-\cos^{2}\theta_{1}}{2h}&&\\\ &\ddots&\\\
&&\dfrac{(x+\frac{2}{\beta}\sin
2\theta_{M})\sin^{2}\theta_{M}-\cos^{2}\theta_{M}}{2h}\end{bmatrix}\check{U},$
(2.8)

where
$\boldsymbol{\theta}_{M}=[\theta_{1},\theta_{2},\dots,\theta_{M}]^{\textsf{T}}$.
Then

$\frac{\partial\boldsymbol{H}_{M}}{\partial
x}=\left[T(\beta,\boldsymbol{\theta}_{M},h)+U(\beta,x,\boldsymbol{\theta}_{M},h)\right]\boldsymbol{H}_{M}.$
(2.9)

Noting that $x$ is the time-like variable, we apply the trapezoidal rule with
time step $\Delta x<0$ to (2.9) yielding

$\displaystyle\frac{\boldsymbol{H}^{n+1}_{M}-\boldsymbol{H}^{n}_{M}}{\Delta
x}=\frac{1}{2}\big{[}T(\beta,\boldsymbol{\theta}_{M},h)\boldsymbol{H}^{n}_{M}+U(\beta,x_{n},\boldsymbol{\theta}_{M},h)\boldsymbol{H}^{n}_{M}\big{]}$
$\displaystyle\hphantom{\frac{\boldsymbol{H}^{n+1}_{M}-\boldsymbol{H}^{n}_{M}}{\Delta
x}=}{}+\frac{1}{2}\big{[}T(\beta,\boldsymbol{\theta}_{M},h)\boldsymbol{H}^{n+1}_{M}+U(\beta,x_{n+1},\boldsymbol{\theta}_{M},h)\boldsymbol{H}^{n+1}_{M}\big{]},$

where $\boldsymbol{H}^{n}_{M}\approx\boldsymbol{H}_{M}(x_{n})$. Upon
rearranging, this gives

$\displaystyle\bigg{[}I-\frac{\Delta
x}{2}T(\beta,\boldsymbol{\theta}_{M},h)-\frac{\Delta
x}{2}U(\beta,x_{n+1},\boldsymbol{\theta}_{M},h)\bigg{]}\boldsymbol{H}^{n+1}_{M}$
$\displaystyle\quad\qquad=\bigg{[}I+\frac{\Delta
x}{2}T(\beta,\boldsymbol{\theta}_{M},h)+\frac{\Delta
x}{2}U(\beta,x_{n},\boldsymbol{\theta}_{M},h)\bigg{]}\boldsymbol{H}^{n}_{M}.$
(2.10)

Finally, we obtain $F_{\beta}(x)\approx\tilde{F}_{\beta}(x):=H_{M}(x)\approx
H^{n}_{M}$, $n=0,1,\dots,N$. Below, we also explore other time integration
methods, in addition to the trapezoidal method. We use the trapezoidal method
as our default due to its A-stability, but it can be outperformed by BDF
methods in this context.

###  2.2 Spectral discretization

To obtain better accuracy, we apply a Fourier spectral method. As before, for
$0\leq n\leq N$, define $x_{n}=x_{0}+n\Delta x$. Suppose
$H(x,\theta)=\int_{0}^{\theta}\rho(x,\theta^{\prime})\,{\rm
d}\theta^{\prime}$, then we rewrite (2.3) as

$\frac{\partial H}{\partial
x}+\bigg{(}\frac{2}{\beta}\sin^{4}\theta\bigg{)}\frac{\partial\rho}{\partial\theta}+\bigg{[}\bigg{(}x+\frac{2}{\beta}\sin
2\theta\bigg{)}\sin^{2}\theta-\cos^{2}\theta\bigg{]}\rho=0.$

Upon taking the derivative with respect to $\theta$ on both sides, we arrive
at a partial differential equation for $\rho(x,\theta)$,

$\displaystyle\frac{\partial\rho}{\partial x}$
$\displaystyle+\bigg{(}\frac{8}{\beta}\sin^{3}\theta\cos\theta\bigg{)}\frac{\partial\rho}{\partial\theta}+\bigg{(}\frac{2}{\beta}\sin^{4}\theta\bigg{)}\frac{\partial^{2}\rho}{\partial\theta^{2}}$
$\displaystyle+\bigg{[}(2x+2)\sin\theta\cos\theta+\frac{2}{\beta}\big{(}2\sin^{2}\theta\cos{2\theta}+2\sin\theta\cos\theta\sin{2\theta}\big{)}\bigg{]}\rho$
$\displaystyle+\bigg{[}\bigg{(}x+\frac{2}{\beta}\sin
2\theta\bigg{)}\sin^{2}\theta-\cos^{2}\theta\bigg{]}\frac{\partial\rho}{\partial\theta}=0.$
(2.11)

Now, suppose

$\rho(x,\theta)\approx\sum_{m=-M}^{M}a_{m}(x){\rm e}^{2{\rm
i}m\theta/l},\qquad\theta\in[0,l\pi),\quad\theta_{M}=l\pi.$ (2.12)

Substituting (2.12) in (2.11) gives a system of ordinary differential
equations for $a_{m}(x)$, $m=-M,\dots,M$, after truncation

$\frac{{\rm d}\boldsymbol{a}_{M}(x)}{{\rm
d}x}=(A+xB)\boldsymbol{a}_{M}(x),\qquad\boldsymbol{a}_{M}(x)=\begin{bmatrix}a_{-M}(x)\\\
\vdots\\\ a_{M}(x)\end{bmatrix},$ (2.13)

where

$\displaystyle A={}$
$\displaystyle-\frac{2}{\beta}\bigg{(}\frac{3}{8}I-\frac{1}{4}S_{-l}-\frac{1}{4}S_{+l}+\frac{1}{16}S_{+2l}+\frac{1}{16}S_{-2l}\bigg{)}D_{2}$
$\displaystyle{}+\bigg{[}\frac{1}{2}I+\frac{1}{4}S_{-l}+\frac{1}{4}S_{+l}-\frac{2}{\beta}\bigg{(}\frac{1}{2{\rm
i}}S_{-l}-\frac{1}{2{\rm
i}}S_{+l}\bigg{)}\bigg{(}\frac{1}{2}I-\frac{1}{4}S_{+l}-\frac{1}{4}S_{-l}\bigg{)}\bigg{]}D_{1}$
$\displaystyle{}-\bigg{[}\frac{8}{\beta}\bigg{(}\frac{1}{-8{\rm
i}}S_{+3l/2}+\frac{1}{8{\rm i}}S_{-3l/2}+\frac{3}{8{\rm
i}}S_{-l/2}-\frac{3}{8{\rm
i}}S_{+l/2}\bigg{)}\bigg{(}\frac{1}{2}S_{-l/2}+\frac{1}{2}S_{+l/2}\bigg{)}\bigg{]}D_{1}$
$\displaystyle{}-\frac{4}{\beta}\bigg{(}\frac{1}{2}I-\frac{1}{4}S_{+l}-\frac{1}{4}S_{-l}\bigg{)}\bigg{(}\frac{1}{2}S_{+l}+\frac{1}{2}S_{-l}\bigg{)}$
$\displaystyle{}-\frac{4}{\beta}\bigg{(}\frac{1}{2{\rm
i}}S_{-l/2}-\frac{1}{2{\rm
i}}S_{+l/2}\bigg{)}\bigg{(}\frac{1}{2}S_{+l/2}+\frac{1}{2}S_{-l/2}\bigg{)}\bigg{(}\frac{1}{2{\rm
i}}S_{-l}-\frac{1}{2{\rm i}}S_{+l}\bigg{)}$
$\displaystyle{}-2\bigg{(}\frac{1}{2{\rm i}}S_{-l/2}-\frac{1}{2{\rm
i}}S_{+l/2}\bigg{)}\bigg{(}\frac{1}{2}S_{+l/2}+\frac{1}{2}S_{-l/2}\bigg{)},$
(2.14)

and

$B=\bigg{(}-\frac{1}{2}I+\frac{1}{4}S_{-l}+\frac{1}{4}S_{+l}\bigg{)}D_{1}-2\bigg{(}\frac{1}{2{\rm
i}}S_{-l/2}-\frac{1}{2{\rm
i}}S_{+l/2}\bigg{)}\bigg{(}\frac{1}{2}S_{+l/2}+\frac{1}{2}S_{-l/2}\bigg{)}.$
(2.15)

Here $S_{\pm k}$ represent the (Fourier modes) shift matrices, $S_{\pm
k}=S^{k}_{\pm 1}$, where

$\displaystyle S_{+1}=\begin{bmatrix}0&1&&&&0\\\ &0&1\\\ &&0&1\\\
&&&\ddots&\ddots\\\ &&&&0&1\\\ 1&&&&&0\end{bmatrix},\qquad
S_{-1}=\begin{bmatrix}0&&&&&1\\\ 1&0&&&\\\ 0&1&0&&\\\ &&\ddots&\ddots&\\\
&&&1&0&0\\\ 0&&&&1&0\end{bmatrix}.$

Also $D_{1}$ and $D_{2}$ represent the first and second-order differentiation
matrices in the Fourier space, i.e., $D_{2}=D^{2}_{1}$, where

$\displaystyle D_{1}=\begin{bmatrix}-2{\rm i}M/l&&&&\\\ &-2{\rm
i}(M-1)/l&&&\\\ &&\ddots&&\\\ &&&2{\rm i}(M-1)/l&\\\ &&&&&2{\rm
i}M/l\end{bmatrix}.$

The Fourier coefficients of the initial condition are obtained via

$a_{m}(x_{0})=\frac{1}{l\pi}\int_{0}^{l\pi}\rho(x_{0},\theta){\rm e}^{-2{\rm
i}m\theta/l}\,{\rm d}\theta.$ (2.16)

Instead of applying the second-order accurate trapezoidal rule to integrate
the system of ODEs for $\rho(x,\theta)$, we suggest the use of a five-step
backward differentiation formula method (BDF5), see [12, Section 8.4]. To use
BDF5, we need four more starting conditions, i.e.,
$\boldsymbol{a}_{M}(x_{0}-i\Delta x)$, $i=1,2,3,4$, which can be obtained by
(2.4) and (2.16). We then proceed with BDF5 to solve for
$\boldsymbol{a}^{n}_{M}\approx\boldsymbol{a}_{M}(x_{n})$

$\displaystyle\left[137I-60\Delta
x\left(A+x_{n}B\right)\right]\boldsymbol{a}^{n}_{M}=300\boldsymbol{a}^{n-1}_{M}$
$\displaystyle-300\boldsymbol{a}^{n-2}_{M}+200\boldsymbol{a}^{n-3}_{M}-75\boldsymbol{a}^{n-4}_{M}+12\boldsymbol{a}^{n-5}_{M}.$
(2.17)

Other time integration methods can be used, and this will be discussed further
below.

Finally, the value of $H(x,\theta)$ can be recovered from

$H(x,\theta)\approx\int_{0}^{\theta}\sum_{m=-M}^{M}a_{m}(x){\rm e}^{2{\rm
i}m\theta^{\prime}/l}\,{\rm
d}\theta^{\prime}=\sum_{m=-M}^{M}a_{m}(x)\int_{0}^{\theta}{\rm e}^{2{\rm
i}m\theta^{\prime}/l}\,{\rm d}\theta^{\prime}.$

The Tracy–Widom distribution can then be approximated by setting $\theta=\pi$,

$F_{\beta}(x)\approx\tilde{F}_{\beta}(x)=\sum_{m=-M}^{M}a_{m}(x)\int_{0}^{\pi}{\rm
e}^{2{\rm i}m\theta^{\prime}/l}\,{\rm
d}\theta^{\prime}=\sum_{m=-M}^{M}a_{m}(x)\frac{l}{2{\rm i}m}\big{(}{\rm
e}^{2{\rm i}m\pi/l}-1\big{)},$

from which we find

$F_{\beta}(x_{n})\approx\tilde{F}_{\beta}(x_{n})\approx\sum_{m=-M}^{M}a^{n}_{m}\frac{l}{2{\rm
i}m}\big{(}{\rm e}^{2{\rm i}m\pi/l}-1\big{)},\qquad n=0,1,\dots,N.$

##  3 Algorithm validation and comparison

At this point, we have approximated the values of $F_{\beta}(x)$ on equally-
spaced grid points $x_{0},x_{1},\allowbreak\ldots,x_{N}$. To obtain a
continuous function, we perform Fourier interpolation. Since Fourier
interpolation has high accuracy when the function being interpolated is
periodic, we define

$\varphi(x)=\frac{\mathrm{erf}(x)+1}{2},\qquad x\in\mathbb{R},$

where $\mathrm{erf}$ denotes the error function [18]. Then consider

$G_{\beta}(x):=\tilde{F}_{\beta}(x)-\varphi(x),$

which is nearly a periodic function. Note that other functions can also be
used instead of the error function $\mathrm{erf}(x)$ in constructing a
periodic function. We perform Fourier interpolation to $G_{\beta}(x)$ on the
grid $x_{1},\dots,x_{N}$. The resulting interpolant is then evaluated on a
shifted and scaled Chebyshev grid on $[x_{N},x_{0}]$. By adding back in
$\varphi(x)$ evaluated on the same Chebyshev grid and then interpolating with
Chebyshev polynomials, we obtain a useful high-accuracy approximation of
$F_{\beta}(x)$. The number of Chebyshev coefficients is user-decided (we use
$10^{3}$ by default).

We point out that using either (2.5) or (2.13), we can also obtain the
approximation of $F^{\prime}_{\beta}(x)$ on the grid
$x_{0},x_{1},\ldots,x_{N}$ at nearly no extra cost. We then perform the same
procedure of interpolation to get an approximation of $F^{\prime}_{\beta}(x)$,
without using $\varphi(x)$ since the function being interpolated is already
nearly periodic.

Algorithm 1 shows the pseudocode for the Fourier interpolation, where T is the
grid points $x_{N},x_{N-1},\dots,x_{0}$ with grid spacing $\Delta x$, $K$ is
an integer that gives the number of Chebyshev coefficients, F is the value of
$\tilde{F}_{\beta}(x_{j})$, $j=1,\dots,N$, and f is the value of
$\tilde{F}^{\prime}_{\beta}(x_{j})$. We think of F, f as functions on T and
need to extend them to all of $[x_{N},x_{0}]$.

Algorithm 1 Fourier$\underline{~{}}$interp(T, $K$, F, f)

1:Set $\texttt{D=}[\min\texttt{T},\max\texttt{T}]$

2:Set
$\texttt{G}(\texttt{T})\texttt{=F}(\texttt{T})\texttt{-}\varphi(\texttt{T})$

3:Perform Fourier interpolation to $\texttt{G}(\texttt{T})$ and f(T) to obtain
the interpolant $\texttt{G}(x)$ and $\texttt{f}(x)$, $x\in\texttt{D}$

4:Set $\tilde{\texttt{F}}\texttt{=G}(\bar{x})\texttt{+}\varphi(\bar{x})$ and
$\tilde{\texttt{f}}(\bar{x})\texttt{=f}(\bar{x})$, where $\bar{x}$ is a
Chebyshev grid on D with $K$ points

5:Perform Chebyshev interpolation to $\tilde{\texttt{F}}$ and
$\tilde{\texttt{f}}$ and return the interpolants

For the finite-difference discretization using trapezoidal method to compute
either the cumulative distribution function (cdf) or the probability density
function (pdf) of $F_{\beta}$ using TW$(\beta)$, as provided in
TracyWidomBeta.jl, the following default values for the parameters are used:

$x_{0}=\bigg{\lfloor}\frac{13}{\sqrt{\beta}}\bigg{\rfloor},\quad
x_{N}=-10,\quad\Delta x=-10^{-3},\quad\theta_{M}=\pi,\quad K=10^{3},\quad
M=10^{3},$ (3.1)

where $\lfloor\cdot\rfloor$ denotes the floor function.

######  Notation.

Throughout we use TW($\beta$; params) to refer to our implementation with
different choices of parameters. For example,

$\text{{TW}$(\beta;\texttt{method="finite"},\texttt{step="trapz"},\texttt{pdf=true})$},$

refers to using the finite-difference discretization, time-stepped with the
trapezoidal method, outputting an approximation of $F^{\prime}_{\beta}$ with
the default parameters (3.1).

The default values of $x_{0}$ and $x_{N}$ are chosen to be $\lfloor
13/\sqrt{\beta}\rfloor$ and $-10$ respectively so that
$F_{\beta}(x_{0})\approx 1$ and $F_{\beta}(x_{N})\approx 0$ for $\beta\geq 1$.
Though not optimal, larger values of $x_{0}$ and smaller values of $x_{N}$ can
also be used. See Section 3.3.3 for a discussion on the selection of the
default value for $x_{0}$. For $0<\beta<1$, a larger domain for $x$ should be
used. The values of $\Delta x$ and $M$ are chosen so that $M=\lfloor{-1/\Delta
x}\rfloor$. In this way, the local truncation error of the trapezoidal method
is of the optimal order. Algorithm 2 shows the pseudocode for
TW$(\beta;\texttt{pdf=true})$. Step 6 may be replaced by solving an alternate
discretization of (2.9).

Algorithm 2 TW($\beta$; pdf=true)

1:Set the values of $x_{0}$, $x_{N}$, $\Delta x$, $\theta_{\text{M}}$, K, and
M as in (3.1)

2:Set $x\texttt{=}x_{0}:\Delta x:x_{N}$

3:Set up the initial condition as in (2.4)

4:Set up the matrices as in (2.6), (2.7), and (2.8)

5:for $\texttt{m=}1,\ldots,\texttt{length}(x)\texttt{-}1$ do

6: Solve (2.10) and denote the result by H

7: Use (2.5) and denote the result by h

8:end for

9:Take the values of H and h for $\theta\texttt{=}\pi$ and denote the results
by F and f

10:Perform Algorithm 1 with F, f, and $\texttt{T=}x$

11:Return the interpolant for $\tilde{\texttt{f}}$

Similarly, for the spectral discretization using BDF5 method to compute either
the cdf or the pdf using

$\displaystyle\begin{split}&\text{{TW}$(\beta;\texttt{method="spectral"},\texttt{step="bdf5"})$}\quad\text{or}\\\
&\text{{TW}$(\beta;\texttt{method="spectral"},\texttt{step="bdf5"},\texttt{pdf=true})$},\end{split}$

the following values for the parameters are used:

$\displaystyle x_{0}=\bigg{\lfloor}\frac{13}{\sqrt{\beta}}\bigg{\rfloor},\quad
x_{N}=-10,\quad\Delta x=-10^{-3},\quad\theta_{M}=20\pi,\quad K=10^{3},\quad
M=8\times 10^{3}.$ (3.2)

One thing to note here is that $\theta_{M}$ is set to be $20\pi$, which is
much larger than $\pi$ as used for finite-difference discretization. This
setting is necessary since a periodic boundary condition (2.12) is imposed on
$\rho(x,\theta)$, or equivalently, on $H(x,\theta)$. If the length of the
domain, $\theta_{M}$, is not large enough, due to periodicity, the
approximation to $H(x,\theta)$ will propagate to the end of the domain and
reappear at the lower boundary $\theta=0$. In other words, $\theta_{M}$
depends on the speed of propagation of the approximation to $H(x,\theta)$, and
it turns out that setting $\theta_{M}=20\pi$ is sufficient. $M$ is set to be
$8\times 10^{3}$ regardless of the value of $\Delta x$ to make sure we have
large enough number of Fourier modes to represent the initial condition. See
Section 3.4 for details on which value of $M$ to use for each method in terms
of both accuracy and computation time. Algorithm 3 shows the pseudocode for
TW$(\beta;\texttt{method="spectral"},\texttt{step="bdf5"},\texttt{pdf=true})$.
As with the finite-difference method, step 11 can be replaced with solving an
alternate discretization of (2.13).

Algorithm 3 TW($\beta$; method="spectral", step="bdf5", pdf=true)

1:Set the values of $x_{0}$, $x_{N}$, $\Delta x$, $\theta_{\text{M}}$, K, and
M as in (3.2)

2:Set $x\texttt{=}x_{0}:\Delta x:x_{N}$

3:Set D=$0:h:\theta_{\text{M}}-h$ with $h\texttt{=}\theta_{\text{M}}/\text{M}$

4:Set up the initial conditions using (2.4) and (2.16).

5:Construct matrices A and B as in (2.14) and (2.15)

6:Set
$\texttt{N=-}\texttt{floor}(\text{M}/2):1:\texttt{floor}(({\text{M}-}1)/2)$

7:for $\texttt{m}\in\texttt{N}$ do

8: Compute the vector $\texttt{integ=}\frac{l}{2{\rm i}m}\big{(}{\rm e}^{2{\rm
i}m\pi/l}\texttt{-}1\big{)}$ with $l\texttt{=}\theta_{\text{M}}/\pi$

9:end for

10:for $\texttt{k=}2:\texttt{length}(x)$ do

11: Solve (2.17) to find $\textbf{a}^{n}_{\text{M}}$

12: Use (2.13) and denote the result by $\textbf{b}^{n}_{\text{M}}$

13:end for

14:Compute
F=$\big{\langle}\textbf{a}^{n}_{\text{M}},\texttt{integ}\big{\rangle}$ and
f=$\big{\langle}\textbf{b}^{n}_{\text{M}},\texttt{integ}\big{\rangle}$

15:Perform Algorithm 1 with F, f, and $\texttt{T=}x$

16:Return the interpolant for $\tilde{\texttt{f}}$

######  Remark 3.1.

To optimize the use of the spectral method, one should consider $\theta$ from
the interval $[k_{1}(x)\pi,k_{2}(x)\pi]$, where $k_{1}$, $k_{2}$ are
$x$-dependent integers that adapt to where the solution is non-zero, see
Figure 10.

Note that when $\beta$ is very large, numerical instabilities can occur. This
instability arises primarily because the initial condition closely resembles a
step function. To get an accurate result in this case, first of all, one needs
to use finer grid for $x$ and larger value for $M$, i.e., refinement in both
time and space. Also, one needs to use more Chebyshev coefficients. With
current values for the parameters, TW$(\beta)$ and
TW$(\beta;\texttt{method="spectral"},\texttt{step="bdf5"})$ exhibits stability
roughly for $1\leq\beta\leq 30$.

###  3.1 Eigenvalues of $\boldsymbol{\Delta
x[T(\beta,\theta_{M},h)+U(\beta,x,\theta_{M},h)]}$
for finite-difference discretization

(a) Zoomed-out view

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

(b) Zoomed-in view

Figure 1: The absolute stability region of trapezoidal method applied to the
finite-difference discretization (the maroon-shaded region along with the
orange boundary curve) and the eigenvalues of $\Delta x(T+U)$ (the red dots
with white boundary) for $\beta=2$, $x=x_{0},x_{0}-1,\dots,x_{N}$, $\Delta
x=-10^{-3}$, and $M=10^{3}$.

(a) Zoomed-out view

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

(b) Zoomed-in view

Figure 2: The absolute stability region of BDF3 applied to the finite-
difference discretization (the maroon-shaded region along with the orange
boundary curve) and the eigenvalues of $\Delta x(T+U)$ (the red dots with
white boundary) for $\beta=2$, $x=x_{0},x_{0}-1,\dots,x_{N}$, $\Delta
x=-10^{-3}$, and $M=10^{3}$.

The default time-stepping routine for the finite-difference discretization is
the trapezoidal method but BDF3, BDF4, BDF5, and BDF6 can also be used on
(2.5) to solve for $\boldsymbol{H}_{M}$. It turns out that with values for the
parameters as in (3.1), convergence occurs for all these methods. Since with
the finite-difference discretization, the trapezoidal method and BDF3 are
usually used, Figures 1 and 2 show the absolute stability regions of
trapezoidal method and BDF3 along with eigenvalues of $\Delta x(T+U)$ for
$\beta=2$, $x=x_{0},x_{0}-1,\dots,x_{N}$, $\Delta x=-10^{-3}$, and $M=10^{3}$.
The eigenvalues in the left-half plane are firmly within the region of
absolute stability in each case.

###  3.2 Eigenvalues of $\boldsymbol{\Delta x(A+xB)}$ for the spectral
discretization

(a) Zoomed-out view

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

(b) Zoomed-in view

Figure 3: The absolute stability region of BDF5 applied to the spectral
discretization (the maroon-shaded region along with the orange boundary curve)
and the eigenvalues of $\Delta x(A+xB)$ (the red dots with white boundary) for
$\beta=2$, $x=x_{0},x_{0}-1,\dots,x_{N}$, $\Delta x=-10^{-3}$, and $M=8\times
10^{3}$.

(a) Zoomed-out view

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

(b) Zoomed-in view

Figure 4: The absolute stability region of BDF6 applied to the spectral
discretization (the maroon-shaded region along with the orange boundary curve)
and the eigenvalues of $\Delta x(A+xB)$ (the red dots with white boundary) for
$\beta=2$, $x=x_{0},x_{0}-1,\dots,x_{N}$, $\Delta x=-10^{-3}$, and $M=8\times
10^{3}$.

Figures 3 and 4 show the absolute stability regions of BDF5 and BDF6 along
with eigenvalues of $\Delta x(A+xB)$ for $\beta=2$,
$x=x_{0},x_{0}-1,\dots,x_{N}$, $\Delta x=-10^{-3}$, and $M=8\times 10^{3}$. It
is clear that, with the values of the parameters given as in (3.2), the
eigenvalues in the left-half plane are firmly within the region of absolute
stability for these methods. We choose BDF5 over BDF6 because with the same
values of these parameters, the performance is roughly the same yet BDF5 has a
larger absolute stability region. Selecting $\Delta x$ is more nuanced than
one might think, see Appendix A for more details.

###  3.3 Error analysis

We now compare the accuracy of these two algorithms when computing the cdf for
${\beta=1,2,4}$. For reference solutions, we implement the Fredholm
determinant representations for $\beta=1,2,4$ [3, 4, 9] by porting the code in
the Julia package RandomMatrices to Mathematica and implementing it in high-
precision arithmetic. This ensures that our reference solutions are accurate
to beyond the $\approx 10^{-16}$ machine precision for standard double
precision. Observe that for $F_{4}(x)$, $x$ should be divided by a factor of
$2^{1/6}$. This difference in variance convention has also been underscored in
[16, p. 47], [14, 17], and [1, Remark 5.1.4].

####  3.3.1 Error across the domain

(a) FD (Trapz)

Error

$x$

Error

$x$

(b) SD (BDF5)

Figure 5: Absolute errors from the two algorithms are presented over the
domain $x\in[-10,13]$ for $\beta=1,2,4$ with $x_{0}=\lfloor
13/\sqrt{\beta}\rfloor$ and $\Delta x=-10^{-3}$. For the finite-difference
discretization, $M=10^{3}$ is used, while for the spectral discretization,
$M=8\times 10^{3}$ is employed.

With values of the parameters given as in (3.1) and (3.2), Figure 5 shows the
absolute errors of the two algorithms across the domain. Tables 1 and 2 show
the absolute errors of some selected $x$-values with same values for the
parameters.

For finite-difference discretization, from either Figure 5 (a) or Table 1, the
error is roughly on the order of $10^{-7}$ around the peak of the
distribution, and it improves when $x$ approaches either end of the domain.
This is due to the fact that exact values, $0$ and $1$, for the initial
condition are imposed on the endpoints of the domain, and the Dirichlet
boundary condition ensures that the solution tends to zero.

For spectral discretization, from either Figure 5 (b) or Table 2, the error is
roughly on the order of $10^{-13}$ throughout the entire domain. Regardless of
the order of error, the reason for this discrepancy from finite-difference
discretization near the edges of the domain is that an error is introduced for
the initial condition when we represent $\rho(x,\theta)$ in terms of a
truncated Fourier series.

$x$ $\beta$ | $1$ | $2$ | $4$
---|---|---|---
$-8$ | $1.832(-15)$ | $5.024(-17)$ | $1.950(-17)$
$-6$ | $7.278(-10)$ | $9.602(-12)$ | $4.766(-16)$
$-4$ | $3.715(-7)$ | $3.600(-7)$ | $1.807(-7)$
$-2$ | $6.553(-7)$ | $1.027(-7)$ | $2.020(-6)$
$0$ | $1.295(-6)$ | $4.014(-7)$ | $3.319(-8)$
$2$ | $4.116(-7)$ | $1.137(-8)$ | $6.185(-12)$
$4$ | $3.352(-8)$ | $3.286(-11)$ | $1.565(-14)$
$6$ | $9.575(-10)$ | $3.186(-14)$ | $5.995(-15)$
Table 1: Absolute errors of finite-difference discretization using trapezoidal method with values of the parameters as in (3.1). $x$ $\beta$ | $1$ | $2$ | $4$
---|---|---|---
$-8$ | $5.912(-13)$ | $6.870(-13)$ | $3.509(-13)$
$-6$ | $3.278(-13)$ | $3.470(-14)$ | $2.021(-13)$
$-4$ | $9.057(-13)$ | $4.580(-13)$ | $4.525(-13)$
$-2$ | $2.663(-12)$ | $2.740(-12)$ | $4.809(-12)$
$0$ | $2.162(-12)$ | $1.090(-13)$ | $1.135(-12)$
$2$ | $7.924(-13)$ | $3.514(-13)$ | $9.293(-14)$
$4$ | $5.719(-13)$ | $1.423(-13)$ | $7.527(-14)$
$6$ | $1.819(-13)$ | $7.094(-14)$ | $6.439(-15)$
Table 2: Absolute errors of spectral discretization using BDF5 with values of
the parameters as in (3.2).

####  3.3.2 Order of error

(a) FD discretization

Error

$|\Delta x|$

Error

$|\Delta x|$

(b) Spectral discretization

Figure 6: The order of accuracy for the finite-difference (FD) discretization
and the spectral discretization. (a) The finite-difference discretization. The
errors for trapezoidal (Trapz) method are computed by treating
TW($\beta\texttt{=}2$; $x_{0}=\lfloor 13/\sqrt{\beta}\rfloor$, step="trapz",
interp=false, $\Delta x\texttt{=}-10^{-3}$, $M\texttt{=}10^{3}$) as the
reference and compared with $\Delta x=-0.005,-0.01,-0.02$ at $x=-2$. For BDF
methods, the same reference is used with step changed accordingly. For BDF3,
$\Delta x=-0.01,-0.02,-0.04$ are used. For BDF4, $\Delta x=-0.02,-0.04,-0.05$
are used. For BDF5, $\Delta x=-0.1,-0.2,-0.25$ are used. For BDF6, $\Delta
x=-0.25,-0.3,-0.4$ are used. (b) The spectral discretization. The errors for
trapezoidal (Trapz) method are computed by treating TW($\beta\texttt{=}2$;
$x_{0}=\lfloor 13/\sqrt{\beta}\rfloor$, method="spectral", step="trapz",
interp=false, $\Delta x\texttt{=}-10^{-3}$, $M\texttt{=}8000$) as the
reference and compared with $\Delta x=-0.025,-0.04,-0.05$ at $x=-2$. For BDF
methods, the same reference is used with step changed accordingly. For BDF3
and BDF4, $\Delta x=-0.025,-0.04,-0.05$ are used. For BDF5, $\Delta
x=-0.1,-0.2,-0.25$ are used. For BDF6, $\Delta x=-0.22,-0.25,-0.44$ are used.

(a) FD (Trapz)

Error

$x$

Error

$x$

(b) SD (BDF5)

Figure 7: Change of error of finite-difference discretization using
trapezoidal method with $\Delta x=-0.004,-0.002,-0.001$ and spectral
discretization using BDF5 with $\Delta x=-0.2,-0.1,-0.05$. Values of the other
parameters are used as in (3.1) and (3.2). When applying BDF5 to the spectral
discretization, errors are more pronounced for smaller time steps at the left
edge due to the onset of instability discussed in the appendix.

For $\beta=2$ evaluated at $x=-2$, Figure 6 shows the order of error plots of
finite-difference discretization and spectral discretization. We choose $x=-2$
since it is near the peak of the distribution. With $\Delta
x=-0.004,-0.002,-0.001$ for finite-difference discretization and $\Delta
x=-0.2,-0.1,-0.05$ for spectral discretization, Figure 7 shows the change of
error as the value of $|\Delta x|$ decreases.

For both the finite-difference discretization and the spectral discretization,
from Figure 6, the error has the expected order with respect to $|\Delta x|$.
However, if we decrease the value of $|\Delta x|$ further, each BDF method has
a corresponding range for $\Delta x$ that causes instability. Moreover, once
the value of $\Delta x$ exits this range, the corresponding BDF method becomes
accurate again. See Appendix A for more details.

######  Remark 3.2.

In addition to the trapezoidal and BDF methods, one could consider $A$-stable
and $L$-stable diagonally implicit Runge–Kutta methods as alternative options
for time-stepping. Furthermore, one could even consider adaptive time-
stepping. While these alternative methods can effectively address the
stability concerns arising from the spectrum of $\Delta x(A+xB)$ and $\Delta
x(T+U)$, they do tend to increase the computation time. These methods likely
require more than one linear solve at each time step (possibly one matrix
factorization and multiple uses of it). And our target time step is $|\Delta
x|=10^{-3}$ and, with this time step, BDF5 is stable with the default
parameters, requiring only one (sparse) linear solve per time step. Then the
comparison of methods becomes a complicated competition between more
expensive, higher-order methods that allow a larger time step and less
expensive, but still fairly high-order, methods requiring a smaller time step.
It is possible that savings can be achieved here, but we did not see it in our
experiments.

####  3.3.3 Error with respect to $\boldsymbol{x_{0}}$

(a) FD (Trapz)

Error

$x_{0}$

Error

$x_{0}$

(b) SD (BDF5)

Figure 8: Change of error of finite-difference discretization using
trapezoidal method and spectral discretization using BDF5 with respect to the
value of $x_{0}$. $\beta=1,2,4$ and $x_{0}=13,12,\dots,2$ are used. Values of
the other parameters are the same as in (3.1) and (3.2). For each value of
$\beta$ and $x_{0}$, the error is taken to be the maximum over
$x=-8,-7,\dots,x_{0}$.

Figure 8 shows how the maximum error over $x=-8,-7,\dots,x_{0}$ changes with
respect to the value of $x_{0}$. For $\beta=4$ with spectral discretization
using BDF5, the error for $x_{0}=12,13$ is roughly
$\mathcal{O}\big{(}10^{-11}\big{)}$, which can be improved to
$\mathcal{O}\big{(}10^{-12}\big{)}$ using more Fourier modes. Except for these
two cases, for both methods, as the value of $\beta$ increases, the value of
$\tilde{x}_{0}$, the minimum value of $x_{0}$ that can be used without
affecting the accuracy, decreases. This is exactly what we expect since the
larger value of $\beta$ is, the more concentrated the distribution becomes.
Based on [5, 6], one expects $\tilde{x}_{0}\approx C/\beta^{\frac{2}{3}}$ for
some constant $C$. This suggests a way to choose the optimal value of $x_{0}$
in terms of accuracy and computation time. For our algorithms, we choose
$C=13$ since for $\beta=1$, $F_{\beta}(13)\approx 1$. In practice, instead of
using $\beta^{\frac{2}{3}}$ for the denominator, we are more conservative and
use $\beta^{\frac{1}{2}}$ since, based on Figure 8 (b), when $\beta=4$, using
$x_{0}=13/\beta^{\frac{2}{3}}\approx 5.16$ brings in an error of
$\mathcal{O}\big{(}10^{-9}\big{)}$.

###  3.4 Computation time

The values in Table 3 are obtained by running on a computer with processor:
$11$th Gen Intel(R) Core(TM) i7-11800H $2.30$GHz with $16.0$GB of RAM.

Discretization | Method | Time (default) | Error (default) | Time ($10^{-6}$)
---|---|---|---|---
Finite difference | Trapz |  $12.684$s | $1.027(-7)$ |  $0.203$s
BDF3 |  $10.731$s | $5.461(-8)$ |  $0.138$s
BDF4 |  $10.458$s | $5.457(-8)$ |  $0.142$s
BDF5 |  $10.947$s | $5.538(-8)$ |  $0.141$s
BDF6 |  $11.067$s | $5.539(-8)$ |  $0.140$s
Spectral | Trapz |  $222.135$s | $4.790(-8)$ |  $10.906$s
BDF3 |  $205.888$s | $4.298(-11)$ |  $1.949$s
BDF4 |  $209.175$s | $2.346(-12)$ |  $1.388$s
BDF5 |  $201.819$s | $2.663(-12)$ |  $1.008$s
BDF6 |  $204.295$s | $2.677(-12)$ |  $204.295$s
Table 3: For each discretization in the first column and time-stepping method
in the second column, computation time to get the interpolated cdf for the
Tracy–Widom distribution and the corresponding error at $x=-2$ for $\beta=2$
are recorded in the next two columns. The computation times in the third
column are generated using the default parameters as in (3.1) and (3.2), with
the corresponding errors in the fourth column. If we aim for an error of
$\mathcal{O}\big{(}10^{-6}\big{)}$, the last column shows the corresponding
minimum observed computation times for each discretization and time-stepping
method to achieve this error.

Table 3 shows the computation time with default values of the parameters as in
(3.1) and (3.2) for $\beta=2$ along with the corresponding error at $x=-2$.
The last column provides the computation time if we aim for an error of
$\mathcal{O}\big{(}10^{-6}\big{)}$ at $x=-2$. For finite-difference
discretization, to have the error of $\mathcal{O}\big{(}10^{-6}\big{)}$ at
$x=-2$, $\Delta x=-0.01$ can be used instead of $\Delta x=-0.001$. For
spectral discretization, to have the error of
$\mathcal{O}\big{(}10^{-6}\big{)}$ at $x=-2$, $M=4\times 10^{3}$ can be used
instead of $M=8\times 10^{3}$ with $\Delta x=-0.01$ for trapezoidal method,
$\Delta x=-0.05$ for BDF3, $\Delta x=-0.07$ for BDF4, and $\Delta x=-0.1$ for
BDF5. It turns out that the error using BDF6 will jump from
$\mathcal{O}\big{(}10^{-5}\big{)}$ to $\mathcal{O}\big{(}10^{-12}\big{)}$.
From Table 3, it takes about $204.295$s using BDF6 to have the error of
$\mathcal{O}(10^{-12})$ with parameters as in (3.2). To have the error of
$\mathcal{O}\big{(}10^{-5}\big{)}$, it takes about $0.639$s with $M=5\times
10^{3}$ and $\Delta x=-0.2$. Work to improve the speed of the spectral
discretization is ongoing.

Using the same machine, Bloemendal’s code (see [1, Section 6]), which uses
Mathematica’s NDSolve, takes approximately $3$ seconds to compute the
approximation to the function $H$. We point out that Mathematica 13.2,
generates warnings from NDSolve regarding the errors in the approximate
solution. The accuracy of this approach, while simple, appears to be limited
to a maximum of 7 digits. Yet, the fact that this is indeed so simple, and
works, is an important reminder of the conceptual simplicity of this
representation of the Tracy–Widom distribution function.

Our methods expand upon this, providing both the cdfs and pdfs as output and
producing high-accuracy interpolants while leaving the trade-off between
accuracy and speed to the user’s discretion. And since we consider the
general-$\beta$ Tracy–Widom distribution functions as important nonlinear
special functions, developing methods, with the highest possible accuracy, is
critical. See Section 5 for some thoughts on further improvements.

##  4 Additional numerical results

In this section, we present additional plots to demonstrate the power and
flexibility of the code.

###  4.1 Comparison with large random matrices

We verify numerically that the pdf generated by our algorithm agrees with the
model presented by Dumitriu and Edelman [7]. Recall that $H_{n}^{\beta}$ in
(1.2) has (1.1) as the jpdf for its eigenvalues, and the distribution of its
largest eigenvalue, after rescaling, converges to $F_{\beta}$ for any
$\beta>0$ as $n\to\infty$.

Histograms in Figure 9 are the normalized histograms for
$n^{1/6}\big{(}\lambda_{\text{max}}\big{(}H^{\beta}_{n}\big{)}-2\sqrt{n}\big{)}$,
where $\lambda_{\text{max}}\big{(}H^{\beta}_{n}\big{)}$ denotes the largest
eigenvalue of the $\beta$-Hermite ensemble $H^{\beta}_{n}$.

(a) $\beta=3$

$F_{3}^{\prime}(x)$

$x$

$F_{5}^{\prime}(x)$

$x$

$F_{6}^{\prime}(x)$

$x$

$F_{7}^{\prime}(x)$

$x$

(b) $\beta=5$ (c) $\beta=6$ (d) $\beta=7$

Figure 9: The histograms for
$n^{1/6}\big{(}\lambda_{\text{max}}\big{(}H^{\beta}_{n}\big{)}-2\sqrt{n}\big{)}$
using $10^{6}$ samples with $n=10^{4}$. The density $F^{\prime}_{\beta}(x)$ is
generated by TW$(\beta;\texttt{pdf=true})$ for $\beta=3,5,6,7$. The small
positive bias in each histogram can be improved by using larger values of $n$.

(a) $\beta=3$

$\theta$$x$$\theta$$x$$\theta$$x$$\theta$$x$

(b) $\beta=5$ (c) $\beta=6$ (d) $\beta=7$

Figure 10: The evolution of $\rho(x,\theta):=\frac{\partial
H}{\partial\theta}$ as $x$ decreases from $x=0$ to $-10$ with step size $=-2$
for $\beta=3,5,6,7$. Finite-difference discretization with trapezoidal method
is used on (2.11) with $\Delta x=-10^{-3}$ and $M=10^{4}$ for
$\theta\in[0,10\pi]$.

###  4.2 Evolution of the density $\boldsymbol{\rho(x,\theta):=\partial
H(x,\theta)/\partial\theta}$

Figure 10 shows the waterfall plots of the approximation of
$\rho(x,\theta):=\partial{H}(x,\theta)/\partial{\theta}$ for
$\theta\in[0,10\pi]$ and $x=0,-2,\dots,-10$ using finite-difference
discretization with trapezoidal method on (2.11) with $\theta_{M}=10\pi$ and
$M=10^{4}$. Values of the other parameters are the same as in (3.1). As the
value of $x$ decreases, the density of the solution propagates to the right.

Figure 11 shows the contour plots of the approximation of $H(x,\theta)$ using
finite-difference discretization with trapezoidal method with values of the
parameters given as in (3.1). For each contour plot, the initial condition
$H(x_{0},\theta)$ is found along the top of the plot and the approximate
Tracy–Widom distribution, $F_{\beta}(x)\approx H(x,\theta_{M})$, is obtained
on the right-hand side of the plot.

###  4.3 Density of the Tracy–Widom distribution

Figure 12 (a) shows the plot of the approximation of $F_{\beta}^{\prime}$ for
$\beta=1$ to $10$ and $x\in[-6,4]$ using finite-difference discretization with
trapezoidal method with values of the other parameters as in (3.1). As we can
see from the plot, as the value of $\beta$ increases, $F_{\beta}^{\prime}$
becomes more concentrated, and its peak moves leftwards (See Appendix B.9 for
a discussion of the exact limiting behavior as $\beta\to\infty$).

Figure 12 (b) is a two-dimensional version of Figure 12 (a), which is also
generated using the finite-difference discretization with trapezoidal method.
It provides a closer view of $F_{\beta}^{\prime}$ for $\beta=1$ to $4$ with
step size $=0.2$. The red curves from right to left correspond to
$\beta=1,2,4$ respectively. The black curves show $F_{\beta}^{\prime}$ for the
other values of $\beta$.

(a) $\beta=3$

$\theta$

$x$

$\theta$

$x$

$\theta$

$x$

$\theta$

$x$

(b) $\beta=5$ (c) $\beta=6$ (d) $\beta=7$

Figure 11: The evolution of $H(x,\theta)$ as $x$ decreases from $\lfloor
13/\sqrt{\beta}\rfloor$ to $-10$ for $\beta=3,5,6,7$. Finite-difference
discretization with trapezoidal method is used with $\Delta x=-10^{-3}$ and
$M=10^{3}$.

(a) $\beta=1$ to $10$

$x$$\beta$

$F_{\beta}^{\prime}(x)$

$x$

(b) $\beta=1$ to $4$

Figure 12: Approximation of the density of the Tracy–Widom distribution for
different values of $\beta$ using the finite-difference discretization with
trapezoidal method with values for the other parameters as in (3.1).

###  4.4 Limiting densities of other eigenvalues

(a) $\beta=3$

$x$$x$$x$$x$

(b) $\beta=5$ (c) $\beta=6$ (d) $\beta=7$

Figure 13: Densities of $-\Lambda_{0}$, $-\Lambda_{1}$, and $-\Lambda_{2}$ of
$\mathcal{H}_{\beta}$ for $\beta=3,5,6,7$. Finite-difference discretization
with trapezoidal method is used with values of the parameters as in (3.1).

By [1, Theorem 2.4.3], one finds the limiting density of the $k$th largest
eigenvalue, after rescaling, of the $\beta$-Hermite ensemble $H^{\beta}_{n}$
at $\theta=k\pi$.

Using finite-difference discretization with trapezoidal method with values for
the parameters as in (3.1), Figure 13 shows the limiting densities of the
largest three eigenvalues of $H^{\beta}_{n}$, namely, $F^{\prime}(\cdot,\pi)$,
$F^{\prime}(\cdot,2\pi)$, and $F^{\prime}(\cdot,3\pi)$ from right to left
respectively. They can also be interpreted as the densities of $-\Lambda_{0}$,
$-\Lambda_{1}$, and $-\Lambda_{2}$ of the stochastic Airy operator
$\mathcal{H}_{\beta}$. As the value of $k$ increases, the limiting
distribution becomes more concentrated.

##  5 Outlook and open questions

The error analysis in this paper was empirical, as a proof of convergence is
elusive. It is likely that eigenvalue perturbation theory could be used to
help show that the spectrum of the $x$-dependent families of matrices we
consider remains close to, or inside of, the stability regions for the time-
steppers we have chosen — at least for sufficiently small time steps. A
challenge here is to obtain quantitative bounds, making “sufficiently small”
precise, and giving useful conclusions. Furthermore, even if the spectrum is
understood, it is currently not known how to estimate the eigenvector
condition numbers and pseudospectra of these families of matrices.

We also believe it to be possible to achieve better accuracy by adapting the
global spectral method of Olver and Townsend [22] and future work will be in
this direction. A high-precision implementation of this idea could help
corroborate conjectures about the tail behavior of the Tracy–Widom
distribution.

##  Appendix A The range of $\boldsymbol{\Delta x}$ that causes instability

(a) Error at $x=-2$

$|\Delta x|$

Error

$|\Delta x|$

$\langle\delta\left(x,\Delta x\right)\rangle$

$|\Delta x|$

$\langle\mu(x,\Delta x)\rangle$

(b) Mean $\langle\delta\left(x,\Delta x\right)\rangle$ (c) Mean
$\langle\mu(x,\Delta x)\rangle$

Figure 14: The non-monotonic convergence for BDF3 along with the mean fraction
$\langle\delta\left(x,\Delta x\right)\rangle$ of the unstable eigenvalues and
the mean value $\langle\mu(x,\Delta x)\rangle$. (a) TW($\beta\texttt{=}2$;
$x_{0}=\lfloor 13/\sqrt{\beta}\rfloor$, method="spectral", step="bdf3",
interp=false, $\Delta x$, $M\texttt{=}8000$) is implemented with $\Delta
x=-1,\allowbreak-0.5,-0.2,-0.1,-0.05,-0.04,-0.02,-0.01,-0.005,-0.004,-0.002,-0.001,-0.0005,-0.0002,\allowbreak-0.0001$.
Measured errors are capped at one for readability. (b) The blue (red) dots
correspond to the mean fraction $\langle\delta_{l}(x,\Delta
x)\rangle(\langle\delta_{r}(x,\Delta x)\rangle)$ as defined in (A.3) ((A.4)).
(c) The blue (red) dots correspond to the mean value $\langle\mu_{l}(x,\Delta
x)\rangle(\langle\mu_{r}(x,\Delta x)\rangle)$ as defined in (A.1)((A.2)).

(a) Zoomed-out view

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

(b) Zoomed-in view

Figure 15: The absolute stability region of BDF3 applying on spectral
discretization (the maroon-shaded region along with the orange boundary curve)
and the eigenvalues $z$ of $\Delta x(A+xB)$ (the red dots with white boundary
correspond to $\operatorname{Re}(z)>0$, and the pink dots with white boundary
correspond to $\operatorname{Re}(z)<0$) for $\beta=2$,
$x=x_{N},-9,\dots,x_{0}$, $\Delta x=-0.004,-0.002$, and $M=8\times 10^{3}$.

(a) Error at $x=-2$

$|\Delta x|$

Error

$|\Delta x|$

$\langle\delta\left(x,\Delta x\right)\rangle$

$|\Delta x|$

$\langle\mu(x,\Delta x)\rangle$

(b) Mean $\langle\delta\left(x,\Delta x\right)\rangle$ (c) Mean
$\langle\mu(x,\Delta x)\rangle$

Figure 16: The non-monotonic convergence for BDF4 along with the mean fraction
$\langle\delta\left(x,\Delta x\right)\rangle$ of the unstable eigenvalues and
the mean value $\langle\mu(x,\Delta x)\rangle$. The detailed computing process
for each plot is the same as for BDF3, see the caption of Figure 14.

(a) Zoomed-out view

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

(b) Zoomed-in view

Figure 17: The absolute stability region of BDF4 (the maroon-shaded region
along with the orange boundary curve) and the eigenvalues $z$ of $\Delta
x(A+xB)$ (the red dots with white boundary correspond to
$\operatorname{Re}(z)>0$, and the pink dots with white boundary correspond to
$\operatorname{Re}(z)<0$) for $\beta=2$, $x=x_{N},-9,\dots,x_{0}$, $\Delta
x=-0.01,-0.005,-0.004,-0.002$, and $M=8\times 10^{3}$.

(a) Error at $x=-2$

$|\Delta x|$

Error

$|\Delta x|$

$\langle\delta\left(x,\Delta x\right)\rangle$

$|\Delta x|$

$\langle\mu(x,\Delta x)\rangle$

(b) Mean $\langle\delta\left(x,\Delta x\right)\rangle$ (c) Mean
$\langle\mu(x,\Delta x)\rangle$

Figure 18: The non-monotonic convergence for BDF5 along with the mean fraction
$\langle\delta\left(x,\Delta x\right)\rangle$ of the unstable eigenvalues and
the mean value $\langle\mu(x,\Delta x)\rangle$. The detailed computing process
for each plot is the same as for BDF3, see the caption of Figure 14.

(a) Zoomed-out view

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

(b) Zoomed-in view

Figure 19: The absolute stability region of BDF5 (the maroon-shaded region
along with the orange boundary curve) and the eigenvalues $z$ of $\Delta
x(A+xB)$ (the red dots with white boundary correspond to
$\operatorname{Re}(z)>0$, and the pink dots with white boundary correspond to
$\operatorname{Re}(z)<0$) for $\beta=2$, $x=x_{N},-9,\dots,x_{0}$, $\Delta
x=-0.02,-0.01,-0.005,-0.004$, and $M=8\times 10^{3}$.

(a) Error at $x=-2$

$|\Delta x|$

Error

$|\Delta x|$

$\langle\delta\left(x,\Delta x\right)\rangle$

$|\Delta x|$

$\langle\mu(x,\Delta x)\rangle$

(b) Mean $\langle\delta\left(x,\Delta x\right)\rangle$ (c) Mean
$\langle\mu(x,\Delta x)\rangle$

Figure 20: The non-monotonic convergence for BDF6 along with the mean fraction
$\langle\delta\left(x,\Delta x\right)\rangle$ of the unstable eigenvalues and
the mean value $\langle\mu(x,\Delta x)\rangle$. The detailed computing process
for each plot is the same as for BDF3, see the caption of Figure 14.

(a) Zoomed-out view

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

$\operatorname{Im}(z)$

$\operatorname{Re}(z)$

(b) Zoomed-in view

Figure 21: The absolute stability region of BDF6 (the maroon-shaded region
along with the orange boundary curve) and the eigenvalues $z$ of $\Delta
x(A+xB)$ (the red dots with white boundary correspond to
$\operatorname{Re}(z)>0$, and the pink dots with white boundary correspond to
$\operatorname{Re}(z)<0$) for $\beta=2$, ${x=x_{N},-9,\dots,x_{0}}$, $\Delta
x=-0.1,-0.05,-0.04,-0.02,-0.01,-0.005,-0.004$, and $M=8\times 10^{3}$.

When applying BDF methods to the spectral discretization, we find that each
BDF method has a corresponding range for $\Delta x$ that causes instability.
Moreover, once the value of $\Delta x$ exits this range, convergence appears
to occur at the expected rate. As a result, when an error plot is generated
with respect to the value of $\left|\Delta x\right|$, we will observe a non-
monotonic pattern of convergence. Applying a $r$-step linear method with the
form [12, Section 7.3]

$\sum_{j=0}^{r}\alpha_{j}U^{n+j}=k\sum_{j=0}^{r}\beta_{j}f\big{(}U^{n+j},t_{n+j}\big{)}$

to $u^{\prime}=\lambda u$ gives

$\sum_{j=0}^{r}\alpha_{j}U^{n+j}=k\sum_{j=0}^{r}\beta_{j}\lambda U^{n+j}\
\Longrightarrow\ \sum_{j=0}^{r}(\alpha_{j}-z\beta_{j})U^{n+j}=0,$

where $z\equiv k\lambda$. We call
$\sum_{j=0}^{r}(\alpha_{j}-z\beta_{j})\xi^{j}$ the stability polynomial of
this method and denote it by $\pi\left(\xi;z\right)$ [12, Section 7.3]. The
absolute stability region of a general $r$-step method consists of the values
of $z$ such that the roots $\xi_{i}(z)$, $i=1,2,\dots,r$ of $\pi(\xi;z)$
satisfy the following two conditions:

  1. (1)

$|\xi_{i}|\leq 1$ for $i=1,2,\dots,r$,

  2. (2)

If $\xi_{i}$ is a repeated root, then $|\xi_{i}|<1$.

See, for example, [12, Definition 6.2] for the definition. To analyze the
roots of $\pi(\xi;z)$ for each BDF method as a function of $x$ and $\Delta x$,
suppose the eigenvalues of $(A+xB)$ are given by
$\lambda_{1},\lambda_{2},\dots,\lambda_{2M+1}$ and define

$\displaystyle\mu_{l}(x,\Delta x):=\max_{j}\left|\frac{\max_{i}|\xi_{i}(\Delta
x\lambda_{j})|-1}{\Delta x}\right|\qquad\text{for}\ \operatorname{Re}(\Delta
x\lambda_{j})<0,$ $\displaystyle\mu_{r}(x,\Delta
x):=\max_{j}\left|\frac{\max_{i}|\xi_{i}(\Delta x\lambda_{j})|-1}{\Delta
x}\right|\qquad\text{for}\ \operatorname{Re}(\Delta x\lambda_{j})>0.$

Similarly, define

$\displaystyle\delta_{l}(x,\Delta
x):=\frac{1}{2M+1}\\#\big{\\{}j\colon\max_{i}|\xi_{i}(\Delta
x\lambda_{j})|>1,\,\operatorname{Re}(\Delta x\lambda_{j})<0\big{\\}},$
$\displaystyle\delta_{r}(x,\Delta
x):=\frac{1}{2M+1}\\#\big{\\{}j\colon\max_{i}|\xi_{i}(\Delta
x\lambda_{j})|>1,\,\operatorname{Re}(\Delta x\lambda_{j})>0\big{\\}},$

where $A+xB$ is the coefficient matrix of the spectral discretization and
$\\#$ gives the cardinality of the set. The function $\mu_{l}(x,\Delta
x)\left(\mu_{r}(x,\Delta x)\right)$ helps capture the maximum growth rate
instabilities in the numerical method caused by the eigenvalues of $\Delta
x(A+xB)$ in the left(right)-half plane and $\delta_{l}(x,\Delta
x)\left(\delta_{r}(x,\Delta x)\right)$ gives the fraction of eigenvalues of
$\Delta x(A+xB)$ in the left(right)-half plane that give rise to a positive
growth rate.

We then compute the mean value $\langle\mu_{l}(x,\Delta x)\rangle$,
$\langle\mu_{r}(x,\Delta x)\rangle$, $\langle\delta_{l}(x,\Delta x)\rangle$,
and $\langle\delta_{r}(x,\Delta x)\rangle$ respectively over
$x=-10,-9,\dots,\lfloor 13/\sqrt{\beta}\rfloor$:

$\displaystyle\langle\mu_{l}(x,\Delta x)\rangle:=\frac{1}{10+\lfloor
13/\sqrt{\beta}\rfloor}\sum_{x=-10}^{\lfloor
13/\sqrt{\beta}\rfloor}\mu_{l}(x,\Delta x),$ (A.1)
$\displaystyle\langle\mu_{r}(x,\Delta x)\rangle:=\frac{1}{10+\lfloor
13/\sqrt{\beta}\rfloor}\sum_{x=-10}^{\lfloor
13/\sqrt{\beta}\rfloor}\mu_{r}(x,\Delta x),$ (A.2)
$\displaystyle\langle\delta_{l}(x,\Delta x)\rangle:=\frac{1}{10+\lfloor
13/\sqrt{\beta}\rfloor}\sum_{x=-10}^{\lfloor
13/\sqrt{\beta}\rfloor}\delta_{l}(x,\Delta x),$ (A.3)
$\displaystyle\langle\delta_{r}\left(x,\Delta
x\right)\rangle:=\frac{1}{10+\lfloor
13/\sqrt{\beta}\rfloor}\sum_{x=-10}^{\lfloor
13/\sqrt{\beta}\rfloor}\delta_{r}(x,\Delta x).$ (A.4)

Figures 14–21 show the non-monotonic convergence for each BDF method along
with the mean fraction $\langle r\left(x,\Delta x\right)\rangle$ of the
unstable eigenvalues of $\Delta x(A+xB)$, the value of $\langle\mu(x,\Delta
x)\rangle$, and the eigenvalues of $\Delta x(A+xB)$ within the unstable range
for $\Delta x$. On one hand, it is not the unstable eigenvalues in the right-
half plane causing the non-monotonic convergence. On the other hand, the range
for $\Delta x$ that causes instability coincides with the range for $\Delta x$
whose values of $\langle\mu(x,\Delta x)\rangle$ in the left-half plane are
larger than $\exp(1)$ and have the largest magnitudes. Moreover, it also
coincides with the range for $\Delta x$ that contains the largest mean
fraction $\langle r(x,\Delta x)\rangle$ of the unstable eigenvalues from the
left-half plane. Note that for BDF6, though the error at $x=-0.1$ does not
blow up, it still suggests that $x=-0.1$ causes instability. For BDF3, the
range for $\Delta x$ that causes instability is approximately
$[-0.004,-0.002]$. For BDF4, the range is approximately $[-0.01,-0.002]$. For
BDF5, the range is approximately $[-0.02,-0.004]$. For BDF6, the range is
approximately $[-0.1,-0.004]$.

##  Appendix B The case when $\boldsymbol{\beta=\infty}$

We consider the limiting behavior of $F_{\beta}$ as $\beta\to\infty$. If we
let $\beta=\infty$, the original boundary value problem becomes

$\frac{\partial F}{\partial x}+\big{(}x-\omega^{2}\big{)}\frac{\partial
F}{\partial\omega}=0\qquad\text{for}\ (x,\omega)\in\mathbb{R}^{2},$ (B.1)

with the following boundary conditions:

$\displaystyle F(x,\omega)\to 1\qquad\text{as}\ x,\omega\to\infty\
\text{together,}$ (B.2) $\displaystyle F(x,\omega)\to 0\qquad\text{as}\
\omega\to-\infty\ \text{with}\ x\ \text{bounded above.}$ (B.3)

Using the method of characteristics [11], rewrite (B.1) as

$\frac{{\rm d}F}{{\rm d}x}=0,$

where

$\frac{{\rm d}F}{{\rm d}x}=\frac{\partial F}{\partial x}+\frac{{\rm
d}\omega}{{\rm d}x}\frac{\partial F}{\partial\omega}$

along curves for which

$\frac{{\rm d}\omega}{{\rm d}x}=x-\omega^{2}.$ (B.4)

Substituting $\omega=u^{\prime}/u$ in (B.4) yields the Airy equation

$u^{\prime\prime}-xu=0,$

whose solution is given by

$u=c_{1}\operatorname{Ai}(x)+c_{2}\operatorname{Bi}(x),$

a linear combination of $\operatorname{Ai}(x)$, the Airy function of the first
kind, and $\operatorname{Bi}(x)$, the Airy function of the second kind [18,
Chapter 9]. Thus, the solution of (B.4) is given by

$\omega=\frac{u^{\prime}}{u}=\frac{c_{1}\operatorname{Ai}^{\prime}(x)+c_{2}\operatorname{Bi}^{\prime}(x)}{c_{1}\operatorname{Ai}(x)+c_{2}\operatorname{Bi}(x)},$
(B.5)

which are the characteristic curves of $F(x,\omega)$.

######  Claim B.1.

There exists a unique characteristic curve for $F(x,\omega)$ along which
$\omega\to-\infty$ as ${x\to\infty}$.

###### Proof.

For existence, consider the characteristic curves

$\tilde{\omega}(x)=\frac{\operatorname{Ai}^{\prime}(x)}{\operatorname{Ai}(x)},$
(B.6)

which is obtained from (B.5) by setting $c_{2}=0$ and $c_{1}\neq 0$. By [18,
Section 9.7 (ii)],

$\lim_{x\to\infty}\tilde{\omega}(x)=\lim_{x\to\infty}\frac{\operatorname{Ai}^{\prime}(x)}{\operatorname{Ai}(x)}=-\infty.$
(B.7)

Therefore, among the characteristic curves represented by (B.6), there exists
one along which $\omega\to-\infty$ as $x\to\infty$.

For uniqueness, note that for fixed values of $c_{1}$ and $c_{2}\neq 0$, as a
well-defined function of $x$, $\omega\to\infty$ as $x\to\infty$. Thus,
combined with (B.7), we conclude that only one such characteristic curve
exists. ∎

After the change of variable $\omega=-\cot\theta$, we obtain

$\cot\theta=-\frac{c_{1}\operatorname{Ai}^{\prime}(x)+c_{2}\operatorname{Bi}^{\prime}(x)}{c_{1}\operatorname{Ai}(x)+c_{2}\operatorname{Bi}(x)},$
(B.8)

and the boundary conditions become

$\displaystyle F(x,\theta)\to 1\qquad\text{as}\ x\to\infty\ \text{and}\
\theta\to\pi,$ $\displaystyle F(x,\theta)\to 0\qquad\text{as}\ \theta\to 0\
\text{with}\ x\ \text{bounded above.}$

(a)

$\theta$

$x$

Figure 22: Contour plot of (B.8) with respect to $x$ and $\theta$ for
$x\in[-10,10]$ and $\theta\in[0,\pi]$. The red curves correspond to the level
zero. The values of $x$ where the red curves cross $\theta=\pi$ coincide with
the zeros of $\mathrm{Ai}(x)$.

Figure 22 shows the contour plot of (B.8) with respect to $x$ and $\theta$ for
$x\in[-10,10]$ and ${\theta\in[0,\pi]}$. The red curves correspond to the
level zero, which occurs when $c_{1}\neq 0$ and $c_{2}=0.$ We can see that
there is one and only one red curve along which $\theta\to 0$ as $x\to\infty$.
Denote the first real zero (closest to $x=0$) of Ai$(x)$ by
$x_{1}\approx-2.33811$. It can be verified from (B.8) by setting $c_{2}=0$ and
$c_{1}\neq 0$ that $x\to x^{+}_{1}$ as $\theta\to\pi^{-}$. Set $F(x,\omega)$
in the following way:

$F(x,\omega)=\begin{dcases}0,&\omega<\dfrac{\operatorname{Ai}^{\prime}(x)}{\operatorname{Ai}(x)},\\\
1,&\omega>\dfrac{\operatorname{Ai}^{\prime}(x)}{\operatorname{Ai}(x)}.\end{dcases}$
(B.9)

Then (B.9) is the solution of (B.1) satisfying both boundary conditions (B.2)
and (B.3).

### Acknowledgements

This work is partially supported by NSFDMS-1945652. The authors would like to
thank the anonymous referees for their helpful comments and suggestions, which
have significantly contributed to the clarity of this paper.

## References

  * [1] Bloemendal A., Finite rank perturbations of random matrices and their continuum limits, Ph.D. Thesis, University of Toronto, Canada, 2011.
  * [2] Bloemendal A., Virág B., Limits of spiked random matrices I, Probab. Theory Related Fields 156 (2013), 795–825, arXiv:1011.1877.
  * [3] Bornemann F., On the numerical evaluation of distributions in random matrix theory: a review, Markov Process. Related Fields 16 (2010), 803–866, arXiv:0904.1581.
  * [4] Bornemann F., On the numerical evaluation of Fredholm determinants, Math. Comp. 79 (2010), 871–915, arXiv:0804.2543.
  * [5] Borot G., Nadal C., Right tail asymptotic expansion of Tracy–Widom beta laws, Random Matrices Theory Appl. 1 (2012), 1250006, 23 pages.
  * [6] Dumaz L., Virág B., The right tail exponent of the Tracy–Widom $\beta$ distribution, Ann. Inst. H. Poincaré Probab. Statist. 49 (2013), 915–933, arXiv:1102.4818.
  * [7] Dumitriu I., Edelman A., Matrix models for beta ensembles, J. Math. Phys. 43 (2002), 5830–5847, arXiv:math-ph/0206043.
  * [8] Edelman A., Sutton B.D., From random matrices to stochastic operators, J. Stat. Phys. 127 (2007), 1121–1165, arXiv:math-ph/0607038.
  * [9] Ferrari P.L., Spohn H., A determinantal formula for the GOE Tracy–Widom distribution, J. Phys. A 38 (2005), L557–L561, arXiv:math-ph/0505012.
  * [10] Grava T., Its A., Kapaev A., Mezzadri F., On the Tracy–Widomβ distribution for $\beta=6$, SIGMA 12 (2016), 105, 26 pages, arXiv:1607.01351.
  * [11] Kevorkian J., Partial differential equations. Analytical solution techniques, Wadsworth & Brooks/Cole Math. Ser., Wadsworth & Brooks/Cole Advanced Books & Software, Pacific Grove, CA, 1990.
  * [12] LeVeque R.J., Finite difference methods for ordinary and partial differential equations. Steady-state and time-dependent problems, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2007.
  * [13] Li Y., On the open question of the Tracy–Widom distribution of $\beta$-ensemble with $\beta=6$, arXiv:1812.00522.
  * [14] Mays A., Ponsaing A., Schehr G., Tracy–Widom distributions for the Gaussian orthogonal and symplectic ensembles revisited: a skew-orthogonal polynomials approach, J. Stat. Phys. 182 (2021), 28, 55 pages, arXiv:2007.14597.
  * [15] Mehta M.L., Random matrices, Pure Appl. Math. (Amsterdam), Vol. 142, 3rd ed., Elsevier/Academic Press, Amsterdam, 2004.
  * [16] Nadal C., Matrices aléatoires et leurs applications à la physique statistique et quantique, Ph.D. Thesis, Université Paris Sud - Paris XI, 2011, available at https://theses.hal.science/tel-00633266.
  * [17] Nadal C., Majumdar S.N., A simple derivation of the Tracy–Widom distribution of the maximal eigenvalue of a Gaussian unitary random matrix, J. Stat. Mech. Theory Exp. 2011 (2011), P04001, 29 pages, arXiv:1102.0738.
  * [18] Olver F.W.J., Olde Daalhuis A.B., Lozier D.W., Schneider B.I., Boisvert R.F., Clark C.W., Miller B.R., Saunders B.V., Cohl H.S., McClain M.A., NIST digital library of mathematical functions, Release 1.1.8 of 2022-12-15, aviable at http://dlmf.nist.gov/.
  * [19] Ramírez J.A., Rider B., Virág B., Beta ensembles, stochastic Airy spectrum, and a diffusion, J. Amer. Math. Soc. 24 (2011), 919–944, arXiv:math.PR/0607331.
  * [20] Rumanov I., Painlevé representation of Tracy–Widomβ distribution for $\beta=6$, Comm. Math. Phys. 342 (2016), 843–868, arXiv:1408.3779.
  * [21] Sutton B.D., The stochastic operator approach to random matrix theory, Ph.D. Thesis, Massachusetts Institute of Technology, 2005.
  * [22] Townsend A., Olver S., The automatic solution of partial differential equations using a global spectral method, J. Comput. Phys. 299 (2015), 106–123, arXiv:1409.2789.
  * [23] Tracy C.A., Widom H., Level-spacing distributions and the Airy kernel, Phys. Lett. B 305 (1993), 115–118, arXiv:hep-th/9210074.
  * [24] Tracy C.A., Widom H., Level-spacing distributions and the Airy kernel, Comm. Math. Phys. 159 (1994), 151–174, arXiv:hep-th/9211141.
  * [25] Tracy C.A., Widom H., On orthogonal and symplectic matrix ensembles, Comm. Math. Phys. 177 (1996), 727–754, arXiv:solv-int/9509007.
  * [26] Valkó B., Virág B., Continuum limits of random matrices and the Brownian carousel, Invent. Math. 177 (2009), 463–508, arXiv:0712.2000.
  * [27] Zhang Y., Trogdon T., TracyWidomBeta, 2023, available at https://github.com/Yiting687691/TracyWidomBeta.jl.