# The Glass Transition of Quantum Hard Spheres in High Dimensions

Michael Winer  Joint Quantum Institute, Department of Physics, University of
Maryland, College Park, Maryland 20742, USA  Christopher L. Baldwin
Department of Physics and Astronomy, Michigan State University, East Lansing,
Michigan 48824, USA  Richard Barney  Joint Quantum Institute, Department of
Physics, University of Maryland, College Park, Maryland 20742, USA  Victor
Galitski  Joint Quantum Institute, Department of Physics, University of
Maryland, College Park, Maryland 20742, USA  Brian Swingle  Department of
Physics, Brandeis University, Waltham, Massachusetts 02453, USA

###### Abstract

We study the equilibrium thermodynamics of quantum hard spheres in the
infinite-dimensional limit, determining the boundary between liquid and glass
phases in the temperature-density plane by means of the Franz-Parisi
potential. We find that as the temperature decreases from high values, the
effective radius of the spheres is enhanced by a multiple of the thermal de
Broglie wavelength, thus increasing the effective filling fraction and
decreasing the critical density for the glass phase. Numerical calculations
show that the critical density continues to decrease monotonically as the
temperature decreases further, suggesting that the system will form a glass at
sufficiently low temperatures for any density.

##  1 Introduction

Much of solid-state physics is built around the idea of periodic or almost-
periodic crystalline matter. This theory has produced enormous triumphs
including the notion of band structure and the Dulong-Petit theory of heat
capacity. However, a large amount of solid matter in both the natural and
human worlds is instead in a glassy phase — amorphous matter with no spatial
symmetry but with a fixed (or so slowly varying as to effectively be fixed)
structure. This phase has a well-studied phenomenology [1, 2, 3, 4], but
little can be derived from first principles. The most successful approaches
have been mean-field theory [5, 6, 7, 8, 9, 10] and mode-coupling theory [11,
12, 13, 14, 15].

An important example of the former is the infinite-dimensional glass, which is
an exactly solvable model that still displays many of the most important
behaviors of glasses, such as arrested dynamics and jamming. The study of the
infinite-dimensional glass has seen much progress in recent years, with
success in both dynamical [16, 17, 18, 19, 20] and thermodynamical [21, 22,
23, 24] methods.

Thus far, most of the work done has been for classical glasses. There has been
significant progress studying quantum fluctuations in related problems — hard-
core bosons on the Bethe lattice [25], the spherical perceptron [26, 27],
“canyon” landscapes [28], and a larger body of work on quantum spin glasses
(see Ref. [29] and references therein) — but quantum effects in the infinite-
dimensional structural glass itself have not yet been considered. The purpose
of the present paper is to do so, particularly for the prototypical model of
hard spheres.

Figure 1: Phase diagram of quantum hard spheres in high dimensions as a
function of inverse temperature $\beta$ and rescaled filling fraction
$\hat{\phi}$ (see Sec. 1.1 for the precise definitions of $\hat{\phi}$ and
$\hat{M}$). The relevant qualitative physics is sketched above — as $\beta$
increases (left to right), quantum fluctuations increase the effective size of
the spheres (see also Ref. [30]) and thus decrease the number of spheres per
unit volume needed for glassiness.

Of particular relevance for us is Ref. [31], which analyzes the mode-coupling
theory (MCT) of quantum hard spheres. MCT is an uncontrolled approximation,
but it allows one to directly study the physically relevant three-dimensional
system. Ref. [31] finds that weak quantum fluctuations increase the effective
radius of the spheres (see the sketch in Fig. 1), and thus the critical
density for glassiness decreases as one lowers from high temperature.
Interestingly, however, they observe that the liquid phase is re-entrant: as
one lowers the temperature, the system will transition from the liquid into
the glass phase but then back into the liquid phase. By contrast, in the
infinite-dimensional system considered here, we find that the critical density
decreases monotonically as a function of temperature. We analytically
determine the leading-order correction to the classical critical density at
high temperature, finding that it is indeed negative (consistent with what is
observed in Ref. [31]), and then numerically determine the phase diagram at
much lower temperatures. The result is shown in Fig. 1. While we cannot reach
arbitrarily low temperature in this manner, we do ultimately consider quite
small values (namely such that the properly scaled de Broglie wavelength is
much larger than the size of the spheres), and still observe a monotonically
decreasing critical density. We do not view these results as incompatible with
those of Ref. [31], however, since the infinite-dimensional limit and MCT
simply entail different approximations. This is rather an opportunity for
further investigation.

###  1.1 Definitions and setup

We are studying a system of (distinguishable) quantum hard spheres in $d$
dimensions, specifically in the large-$d$ limit. Some background discussion is
given in App. A. Our approach is based on the classical work found in, for
example, Refs. [32, 33, 34, 35].

The setup is as follows. We have $N$ spheres labeled by
$i\in\\{1,\cdots,N\\}$, each with radius $\ell/2$ and mass $M$, with momentum
operator $\textbf{P}_{i}$ and position operator $\textbf{X}_{i}$ (and
$\textbf{X}_{ij}\equiv\textbf{X}_{i}-\textbf{X}_{j}$). Bold letters denote
$d$-dimensional vectors. The system is described by a Hamiltonian

$H=\sum_{i}\frac{\textbf{P}_{i}^{2}}{2M}+\sum_{i<j}V(\textbf{X}_{ij}),$ (1)

with inverse temperature $\beta$ and particle density $\rho$. We take the
interactions $V$ to be rotationally symmetric but otherwise completely general
for as long as possible. Eventually, however, we will specialize to an
infinite hard-sphere potential. In this situation, the closest two spheres can
get to each other is $\ell$, with $V(\textbf{x}_{ij})=\infty$ for
$|\textbf{x}_{ij}|<\ell$ and $V(\textbf{x}_{ij})=0$ otherwise.

The filling fraction $\phi$ is defined as the product of $\rho$ with the
volume of an individual sphere. It will be convenient to define
$\hat{\phi}\equiv 2^{d}\phi/d$ and also $\hat{M}\equiv M\ell^{2}/d^{2}$. We
will find a non-trivial phase diagram when $\hat{\phi}$ and $\hat{M}$ are held
constant as $d$ goes to infinity. This means that in high dimensions, a glass
of hard spheres is mostly empty! The interaction potential $V$ must also scale
a certain way with $d$, namely
$V(\textbf{x}_{ij})\equiv\overline{v}(d(\frac{|\textbf{x}_{ij}|}{\ell}-1))$
with $\overline{v}(h)$ a $d$-independent function. Note that a hard-sphere
potential corresponds to $\overline{v}(h)=\infty$ for $h<0$ and
$\overline{v}(h)=0$ for $h>0$.

Our analysis will be entirely in terms of a path-integral description of the
system. Thus we will no longer consider operators $\textbf{P}_{i}$ and
$\textbf{X}_{i}$ but rather trajectories $\textbf{x}_{i}(\tau)$ with imaginary
time $\tau\in[0,\beta)$. The (Euclidean) Lagrangian is given by

$\mathcal{L}=\frac{\hat{M}d^{2}}{2\ell^{2}}\sum_{i}\left|\frac{\textrm{d}\textbf{x}_{i}(\tau)}{\textrm{d}\tau}\right|^{2}+\sum_{i<j}\overline{v}\left(d\Big{(}\frac{|\textbf{x}_{ij}(\tau)|}{\ell}-1\Big{)}\right).$
(2)

###  1.2 Franz-Parisi potential

We study this system using the Franz-Parisi potential (FPP), a method of
extracting dynamical information from purely thermodynamical calculations.
Here we give only a brief overview, and refer to Refs. [33, 36, 37] for
further details. The FPP is defined as the average free energy of the system
when it’s constrained to be at a given mean-squared distance $\Delta$ from a
reference configuration, with the reference configuration itself drawn from
the equilibrium Boltzmann distribution. In equation form, the FPP $\Phi$ is
defined as

$\beta\Phi\equiv-\int\prod_{i}\mathcal{D}\textbf{x}_{i}\frac{e^{-\int\textrm{d}\tau\mathcal{L}[\textbf{x}]}}{Z}\log\int\prod_{i}\mathcal{D}\textbf{y}_{i}e^{-\int\textrm{d}\tau\mathcal{L}[\textbf{y}]}\delta\Big{(}\frac{1}{N}\sum_{i}\big{|}\textbf{x}_{i}(0)-\textbf{y}_{i}(0)\big{|}^{2}-\Delta\Big{)},$
(3)

where $Z$ is the partition function and $\prod_{i}\mathcal{D}\textbf{x}_{i}$
denotes the path-integral measure. We will find that $\Phi$ is unbounded from
below as $\Delta\rightarrow\infty$, but that it can have a local minimum at
finite $\Delta$ for sufficiently large filling fractions. This is the arrested
(i.e., glass) phase — we refer to the above references for discussion of why
this is so. Thus our aim is specifically to determine whether $\Phi$ defined
by Eq. (3) has a stationary point at finite $\Delta$.

Following Ref. [33], we evaluate the FPP using the replica trick and a virial
expansion (see Refs. [6, 9, 37, 38] for details of the former and Refs. [39,
40, 41, 42] for details of the latter). The replica trick allows us to write
Eq. (3) as

$\displaystyle\beta\Phi$
$\displaystyle=-\int\prod_{i}\mathcal{D}\textbf{x}_{i}\frac{e^{-\int\textrm{d}\tau\mathcal{L}[\textbf{x}]}}{Z}\lim_{m\rightarrow
1}\frac{1}{m-1}\left[\left(\int\prod_{i}\mathcal{D}\textbf{y}_{i}e^{-\int\textrm{d}\tau\mathcal{L}[\textbf{y}]}\delta\Big{(}\frac{1}{N}\sum_{i}\big{|}\textbf{x}_{i}(0)-\textbf{y}_{i}(0)\big{|}^{2}-\Delta\Big{)}\right)^{m-1}-1\right]$
(4) $\displaystyle=-\lim_{m\rightarrow
1}\frac{1}{m-1}\left[\frac{1}{Z}\int\prod_{ia}\mathcal{D}\textbf{x}_{ia}e^{-\sum_{a}\int\textrm{d}\tau\mathcal{L}[\textbf{x}_{a}]}\prod_{a>1}\delta\Big{(}\frac{1}{N}\sum_{i}\big{|}\textbf{x}_{i1}(0)-\textbf{x}_{ia}(0)\big{|}^{2}-\Delta\Big{)}-1\right],$

where to obtain the bottom line, $\textbf{x}_{i}$ is redefined as
$\textbf{x}_{i1}$ and the $m-1$ different factors of $\textbf{y}_{i}$ are
redefined as $\textbf{x}_{ia}$ for $a\in\\{2,\cdots,m\\}$ (thus in total $a$
ranges from 1 to $m$). We thus need only to calculate the “replicated”
partition function

$Z_{\textrm{rep}}\equiv\int\prod_{ia}\mathcal{D}\textbf{x}_{ia}e^{-\sum_{a}\int\textrm{d}\tau\mathcal{L}[\textbf{x}_{a}]}\prod_{a>1}\delta\Big{(}\frac{1}{N}\sum_{i}\big{|}\textbf{x}_{i1}(0)-\textbf{x}_{ia}(0)\big{|}^{2}-\Delta\Big{)}.$
(5)

Eq. (4) then states that $Z_{\textrm{rep}}\sim Z[1-(m-1)\beta\Phi]$ for $m$
close to 1, i.e., $-\beta\Phi$ is the first-order coefficient in the expansion
of $\log{Z_{\textrm{rep}}}$ around $m=1$.

We calculate $\log{Z_{\textrm{rep}}}$ via a virial expansion [39, 40, 41, 42].
This expansion is in terms of the single-particle density $\rho(\textbf{x})$,
and thus the resulting free energy includes $\rho(\textbf{x})$ as a
variational function to be minimized over. Since we are considering a
replicated quantum system, $\rho(\textbf{x})$ is still a function of the
coordinates of all replicas at all imaginary times, i.e., the argument x in
$\rho(\textbf{x})$ has replica and imaginary-time indices —
$\textbf{x}\equiv\\{\textbf{x}_{a}(\tau)\\}_{a,\tau}$. Yet importantly, it no
longer has a particle index $i$. Furthermore, taking the large-$d$ limit
allows us to truncate the expansion at second order [33]. The replicated free
energy then amounts to (up to an overall constant concerning normalization of
the path integral)

$\displaystyle\beta f\equiv-\frac{\log{Z_{\textrm{rep}}}}{Nd}$
$\displaystyle=\min_{\rho}\left\\{\frac{1}{Nd}\int\prod_{a}\mathcal{D}\textbf{x}_{a}\rho(\textbf{x})\left[\log{\rho(\textbf{x})}+\frac{M}{2}\sum_{a}\int_{0}^{\beta}\textrm{d}\tau\left|\frac{\textrm{d}\textbf{x}_{a}(\tau)}{\textrm{d}\tau}\right|^{2}\right]\right.$
(6)
$\displaystyle\qquad\left.-\frac{1}{2Nd}\int\prod_{a}\mathcal{D}\textbf{x}_{a}\mathcal{D}\textbf{y}_{a}\rho(\textbf{x})\rho(\textbf{y})\left(\exp\left[-\sum_{a}\int_{0}^{\beta}\textrm{d}\tau
V\big{(}\textbf{x}_{a}(\tau)-\textbf{y}_{a}(\tau)\big{)}\right]-1\right)\right\\}.$

Keep in mind that we still have the constraint
$N^{-1}\sum_{i}|\textbf{x}_{i1}(0)-\textbf{x}_{ia}(0)|^{2}=\Delta$. This is
easy to account for in Eq. (6) — only consider densities $\rho(\textbf{x})$
such that $|\textbf{x}_{1}(0)-\textbf{x}_{a}(0)|^{2}=\Delta$ in the
minimization.

An additional simplification that occurs at large $d$ is that all terms in Eq.
(6) depend only on the “second moment” of $\rho(\textbf{x})$ (interpreting
$\rho(\textbf{x})/N$ as a probability distribution on x) [33]. Thus rather
than minimize over all functions $\rho(\textbf{x})$, we need only minimize
over its covariance matrix.

Even this can be further simplified by making use of the symmetries in the
system. Here it is helpful to think of the $m$ replicas of each sphere as
forming a “molecule”. Let us first consider the corresponding classical system
— each molecule is characterized by $md$ numbers (one for each component of
each $\textbf{x}_{a}$), and in principle $\rho(\textbf{x})$ could depend on
all of them. Yet since the system is translation- and rotation-invariant,
$\rho(\textbf{x})$ in fact depends only on the $m(m+1)/2$ numbers in the
matrix
$(\textbf{x}_{a}-\overline{\textbf{x}})\cdot(\textbf{x}_{b}-\overline{\textbf{x}})$,
where $\overline{\textbf{x}}\equiv m^{-1}\sum_{c}\textbf{x}_{c}$ is the
“center of mass” of the molecule. Furthermore, as noted above, the replicated
free energy depends only on the expectation value of this quantity, denoted
$\alpha_{ab}$:

$\alpha_{ab}\equiv\frac{d}{\ell^{2}}\big{<}(\textbf{x}_{a}-\overline{\textbf{x}})\cdot(\textbf{x}_{b}-\overline{\textbf{x}})\big{>},$
(7)

where the brackets denote an average with respect to $\rho(\textbf{x})$. The
explicit factors of $d$ and $\ell$ are chosen so that $\alpha_{ab}$ is $O(1)$
in the large-$d$ limit. Note that even the entries of $\alpha$ are not
entirely independent, since by definition the all-1 vector has eigenvalue 0:
$\sum_{b}(\textbf{x}_{b}-\overline{\textbf{x}})=0$ and thus
$\sum_{b}\alpha_{ab}=0$. Also note that the constraint amounts to requiring
$\alpha_{11}+\alpha_{aa}-2\alpha_{a1}=d\Delta/\ell^{2}$ for all $a>1$.

Recall that once we calculate the replicated free energy $f$ and then the FPP
$\Phi$, our goal is to identify points where
$\textrm{d}\Phi/\textrm{d}\Delta=0$. At such points, if they exist, we have a
further permutation symmetry between replicas (see Ref. [33]). Thus there are
only two distinct elements of $\alpha$ — the diagonal and off-diagonal entries
— and the fact that $\sum_{b}\alpha_{ab}=0$ fixes one in terms of the other.
We can therefore parametrize $\alpha$ by a single number $A$ (with this
specific parametrization chosen so as to match the expressions in Ref. [33]):

$\alpha=\frac{1}{2m}\begin{pmatrix}(m-1)A&-A&\cdots&-A\\\
-A&(m-1)A&\cdots&-A\\\ \vdots&\vdots&\ddots&\vdots\\\
-A&-A&\cdots&(m-1)A\end{pmatrix}.$ (8)

When we promote from classical to quantum thermodynamics, the $md$ coordinates
of each molecule become $md$ imaginary-time trajectories. The matrix $\alpha$
then acquires time dependence as well:

$\alpha_{ab}(\tau_{1},\tau_{2})\equiv\frac{d}{\ell^{2}}\big{<}\big{(}\textbf{x}_{a}(\tau_{1})-\overline{\textbf{x}}\big{)}\cdot\big{(}\textbf{x}_{b}(\tau_{2})-\overline{\textbf{x}}\big{)}\big{>},$
(9)

where, analogous to above,
$\overline{\textbf{x}}\equiv(m\beta)^{-1}\int\textrm{d}\tau\sum_{c}\textbf{x}_{c}(\tau)$.
However, each replica has a separate time-translation symmetry, meaning that
$\langle\textbf{x}_{a}(\tau_{1})\cdot\textbf{x}_{a}(\tau_{2})\rangle$ depends
only on $\tau_{1}-\tau_{2}\equiv\tau_{12}$ and
$\langle\textbf{x}_{a}(\tau_{1})\cdot\textbf{x}_{b}(\tau_{2})\rangle$ is
independent of both $\tau_{1}$ and $\tau_{2}$. Accounting for replica symmetry
as well, this means that $\alpha$ takes the form

$\alpha(\tau_{1},\tau_{2})=\frac{1}{2m}\begin{pmatrix}(m-1)A+mG(\tau_{12})&-A&\cdots&-A\\\
-A&(m-1)A+mG(\tau_{12})&\cdots&-A\\\ \vdots&\vdots&\ddots&\vdots\\\
-A&-A&\cdots&(m-1)A+mG(\tau_{12})\end{pmatrix}.$ (10)

As before, $\int\textrm{d}\tau_{2}\sum_{b}\alpha_{ab}(\tau_{1},\tau_{2})=0$,
which translates to $\int\textrm{d}\tau G(\tau)=0$. In other words, we are
fixing the zero-frequency component of $G_{\omega}\equiv\int\textrm{d}\tau
e^{i\omega\tau}G(\tau)$ to be identically zero.

###  1.3 Plan for the paper

To summarize, our strategy for determining the phase diagram of the quantum
hard-sphere liquid is as follows:

  * •

Express the replicated free energy given by Eq. (6) in terms of the dynamical
correlation function $G(\tau)$ (equivalently $G_{\omega}$) and the inter-
replica overlap $A$.

  * •

Obtain saddle-point equations for $G_{\omega}$ and $A$ in the limit
$m\rightarrow 1$.

  * •

We will find that the equations for $G_{\omega}$ do not depend on $A$ (and in
fact require only the zeroth-order term of $f$). Thus solve for $G_{\omega}$
first.

  * •

Insert the solution for $G_{\omega}$ into the equation for $A$ (which does
require the first-order term of $f$, i.e., the FPP $\Phi$). We will find that
$A=\infty$ is always a solution. If this is the only solution, then the system
is in the liquid phase. If there is also a finite-$A$ solution, on the other
hand, then the system is in the glass phase.

One might notice that in this procedure, we never actually calculate the FPP
as a function of the mean-squared distance $\Delta$. This is because, as
stated above, determining the phase diagram only requires knowledge of the
stationary points of the FPP. We can perform the extremization over $\Delta$
at the same time as we extremize over the other quantities needed to determine
the free energy, and within the replica-symmetric ansatz described above, this
amounts to extremizing over the parameter $A$.

In the remainder of the paper, we first evaluate the replicated free energy
and determine its saddle-point equations in Sec. 2. We analytically solve the
equations at high temperature in Sec. 3, finding that quantum fluctuations
“renormalize” the radius of the spheres to be slightly larger (and thus the
glass phase appears at a slightly lower density). We then numerically solve
the equations at lower temperatures in Sec. 4. We do not observe any sign of a
re-entrant liquid phase — the critical density instead decreases monotonically
as the temperature lowers. Lastly, we discuss future directions in Sec. 5.

##  2 The free energy and its saddle points

Here we evaluate the replicated free energy given in Eq. (6). It is helpful to
do so in stages, first determining the classical free energy (which will be
entirely a review of known results) and then incorporating quantum
fluctuations.

###  2.1 Classical free energy

Even for the classical free energy, we proceed in stages: first calculate the
unreplicated non-interacting free energy, then the unreplicated interacting
free energy, and lastly the replicated interacting free energy.

In the absence of interactions and without any replicas, the free energy is
trivial to evaluate — in any number of dimensions, the partition function is
$Z=V^{N}/N!$, where $V$ denotes the volume. This gives a free energy
$f_{\textrm{position}}$:

$\beta
f_{\textrm{position}}=-\frac{\log{Z}}{Nd}\sim\frac{\log\rho}{d}\sim\frac{1}{2}\log\frac{d}{2\pi
e\ell^{2}},$ (11)

where we’ve neglected terms that vanish at large $d$ (recall that
$\hat{\phi}\equiv 2^{d}\phi/d$ is kept $O(1)$ at large $d$ and $\phi\sim(\pi
e\ell^{2}/2d)^{d/2}\rho$ — see App. A). Note that $\beta
f_{\textrm{position}}$ is equivalently the contribution of the first term in
Eq. (6) — for a translation-invariant classical system without replicas,
$\rho(\textbf{x})$ is simply a constant $\rho$ and so
$\int\textrm{d}\textbf{x}\rho(\textbf{x})\log{\rho(\textbf{x})}=N\log{\rho}$.
Also note that this contribution depends logarithmically on $d$. Usually the
free energy is proportional to the number of degrees of freedom, but in the
limit we are taking there is this additional term going as $\log d$.

Now include interactions (still without replicas). The correction coming from
the virial expansion is

$\beta
f_{\textrm{virial}}=-\frac{\rho}{2d}\int\textrm{d}\textbf{x}\left(e^{-\beta
V(\textbf{x})}-1\right)=-\frac{\hat{\phi}}{2}\int\textrm{d}he^{h}\left(e^{-\beta\overline{v}(h)}-1\right),$
(12)

where the final expression follows from writing
$V(\textbf{x})=\bar{v}(d(|\textbf{x}|/\ell-1))$, substituting
$h=d(|\textbf{x}|/\ell-1)$, and using that
$|\textbf{x}|^{d-1}\sim\ell^{d-1}e^{h}$ at large $d$. For the case of hard
spheres ($\overline{v}(h)=\infty$ for $h<0$ and $\overline{v}(h)=0$ for
$h>0$), the integral in Eq. (12) evaluates to $-1$ and thus

$\beta f_{\textrm{virial}}=\frac{\hat{\phi}}{2}.$ (13)

Now introduce replicas. As mentioned in Sec. 1, it is convenient to think of
each set of $m$ replicas as forming a single “molecule”. In calculating the
entropic contribution to the free energy (first term of Eq. (6)), we still
obtain a term $\beta f_{\textrm{position}}$ coming from the overall density of
molecules. Yet now the internal entropy of each molecule is non-trivial. As
noted above, the other terms of Eq. (6) depend only on the covariance matrix
$\ell^{2}\alpha/d^{2}$ of the coordinates
$\textbf{x}_{a}-\overline{\textbf{x}}$. Thus in minimizing the free energy, we
should take $\rho(\textbf{x})$ to be the distribution with maximum entropy
subject to fixed covariance, namely Gaussian. Using the standard result for
the entropy of a Gaussian, we obtain a “molecular shape” contribution to the
free energy given by

$\beta f_{\textrm{molecular}}=-\frac{1}{2}\log{\textrm{Det}^{\prime}\frac{2\pi
e\ell^{2}\alpha}{d^{2}}},$ (14)

where $\textrm{Det}^{\prime}$ indicates that the zero eigenvalue of $\alpha$
is omitted. Given the form in Eq. (8), this evaluates to

$\beta f_{\textrm{molecular}}=-\frac{m-1}{2}\log{\frac{\pi
e\ell^{2}A}{d^{2}}}.$ (15)

Calculating the interaction term in the replicated setting requires more
thought. Referring to Eq. (6), we see that we need to average
$\exp[-\beta\sum_{a}V(\textbf{x}_{a}-\textbf{y}_{a})]$ over the coordinates of
two molecules x and y. The relative coordinates
$\textbf{x}_{a}-\overline{\textbf{x}}$ and
$\textbf{y}_{a}-\overline{\textbf{y}}$ are Gaussian-distributed with
covariance matrix $\ell^{2}\alpha/d^{2}$, while $\overline{\textbf{x}}$ and
$\overline{\textbf{y}}$ are uniformly distributed by translation invariance.
Thus define $\textbf{u}_{a}\equiv\textbf{x}_{a}-\overline{\textbf{x}}$ and
$\textbf{v}_{a}\equiv\textbf{y}_{a}-\overline{\textbf{y}}$, as well as
$\textbf{r}\equiv\overline{\textbf{x}}-\overline{\textbf{y}}$. The interaction
term then amounts to

$\beta
f_{\textrm{virial}}=-\frac{\rho}{2d}\int\textrm{d}\textbf{r}\left(\big{<}e^{-\beta\sum_{a}V(\textbf{u}_{a}-\textbf{v}_{a}+\textbf{r})}\big{>}-1\right),$
(16)

where the brackets denote an average over the aforementioned Gaussians u and
v.

For the radially symmetric interactions considered here, so that
$V(\textbf{u}_{a}-\textbf{v}_{a}+\textbf{r})$ depends only on
$|\textbf{u}_{a}-\textbf{v}_{a}+\textbf{r}|$, we can substantially simplify
further. Writing the separation as

$|\textbf{u}_{a}-\textbf{v}_{a}+\textbf{r}|=\sqrt{|\textbf{r}|^{2}+2\textbf{r}\cdot(\textbf{u}_{a}-\textbf{v}_{a})+|\textbf{u}_{a}-\textbf{v}_{a}|^{2}}\sim|\textbf{r}|+\frac{\textbf{r}\cdot(\textbf{u}_{a}-\textbf{v}_{a})}{|\textbf{r}|}+\frac{|\textbf{u}_{a}-\textbf{v}_{a}|^{2}}{2|\textbf{r}|},$
(17)

the term $\textbf{r}\cdot(\textbf{u}_{a}-\textbf{v}_{a})/|\textbf{r}|$ is a
Gaussian with covariance $2\ell^{2}\alpha_{ab}/d^{2}$, and the term
$|\textbf{u}_{a}-\textbf{v}_{a}|^{2}/2|\textbf{r}|$ is deterministic and equal
to $\ell^{2}\alpha_{aa}/d|\textbf{r}|$ at large $d$. Thus

$\big{<}e^{-\beta\sum_{a}V(\textbf{u}_{a}-\textbf{v}_{a}+\textbf{r})}\big{>}=\big{(}\textrm{Det}4\pi\alpha\big{)}^{-\frac{1}{2}}\int\prod_{a}\textrm{d}z_{a}e^{-\frac{1}{4}\sum_{ab}\alpha_{ab}^{-1}z_{a}z_{b}}\prod_{a}e^{-\beta\overline{v}(h+\alpha_{aa}+z_{a})},$
(18)

where we have again substituted $h=d(|\textbf{r}|/\ell-1)$. Using the identity
[33]

$\big{(}\textrm{Det}4\pi\alpha\big{)}^{-\frac{1}{2}}\int\prod_{a}\textrm{d}z_{a}e^{-\frac{1}{4}\sum_{ab}\alpha_{ab}^{-1}z_{a}z_{b}}\prod_{a}e^{-\beta\overline{v}(h+\alpha_{aa}+z_{a})}=e^{\sum_{ab}\alpha_{ab}\frac{\partial^{2}}{\partial
h_{a}\partial
h_{b}}}\prod_{a}e^{-\beta\overline{v}(h_{a})}\bigg{|}_{h_{a}=h+\alpha_{aa}},$
(19)

we can thus write the interaction term as

$\beta
f_{\textrm{virial}}=-\frac{\hat{\phi}}{2}\int\textrm{d}he^{h}\left[e^{\sum_{ab}\alpha_{ab}\frac{\partial^{2}}{\partial
h_{a}\partial
h_{b}}}\prod_{a}e^{-\beta\overline{v}(h_{a})}\bigg{|}_{h_{a}=h+\alpha_{aa}}-1\right].$
(20)

Given the form of $\alpha$ in Eq. (8), $\beta f_{\textrm{virial}}$ further
simplifies to

$\displaystyle\beta f_{\textrm{virial}}$
$\displaystyle=-\frac{\hat{\phi}}{2}\int\textrm{d}he^{h}\left[e^{-\frac{A}{2m}\frac{\textrm{d}^{2}}{\textrm{d}h^{2}}}\prod_{a}e^{\frac{A}{2}\frac{\partial^{2}}{\partial
h_{a}^{2}}}e^{-\beta\overline{v}(h_{a})}\bigg{|}_{h_{a}=h+\frac{(m-1)A}{2m}}-1\right]$
(21)
$\displaystyle=-\frac{\hat{\phi}}{2}\int\textrm{d}h\left[\left(e^{-\frac{A}{2m}\frac{\textrm{d}^{2}}{\textrm{d}h^{2}}}e^{h}\right)\prod_{a}e^{\frac{A}{2}\frac{\partial^{2}}{\partial
h_{a}^{2}}}e^{-\beta\overline{v}(h_{a})}\bigg{|}_{h_{a}=h+\frac{(m-1)A}{2m}}-e^{h}\right]$
$\displaystyle=-\frac{\hat{\phi}}{2}\int\textrm{d}he^{h}\left[\left(e^{\frac{A}{2}\frac{\textrm{d}^{2}}{\textrm{d}h^{2}}}e^{-\beta\overline{v}(h+\frac{A}{2})}\right)^{m}-1\right].$

All in all, the replicated free energy of the classical liquid is given by the
sum of $f_{\textrm{position}}$ (Eq. (11)), $f_{\textrm{molecular}}$ (Eq.
(15)), and $f_{\textrm{virial}}$ (Eq. (21)):

$\beta f_{\textrm{classical}}=\frac{1}{2}\log\frac{d}{2\pi
e\ell^{2}}-\frac{m-1}{2}\log{\frac{\pi
e\ell^{2}A}{d^{2}}}-\frac{\hat{\phi}}{2}\int\textrm{d}he^{h}\left[\left(e^{\frac{A}{2}\frac{\textrm{d}^{2}}{\textrm{d}h^{2}}}e^{-\beta\overline{v}(h+\frac{A}{2})}\right)^{m}-1\right].$
(22)

###  2.2 Quantum free energy

Moving on to the quantum system, it is important to balance two
interpretations. In one, we have $m$ replicas, each completing periodic motion
on the thermal circle with period $\beta$ under the influence of kinetic and
potential energies. Following this approach, one could recalculate all of the
above terms while keeping track of additional imaginary-time indices. The only
wholly new contribution to the free energy is the kinetic term (second term of
the first line in Eq. (6)), which evalutes to

$\beta
f_{\textrm{kinetic}}=\frac{M\ell^{2}}{2d^{2}}\sum_{a}\int_{0}^{\beta}\textrm{d}\tau\frac{\partial^{2}}{\partial\tau_{1}\partial\tau_{2}}\alpha_{aa}(\tau_{1},\tau_{2})\Big{|}_{\tau_{1},\tau_{2}=\tau}=-\frac{\hat{M}m\beta}{4}G^{\prime\prime}(0).$
(23)

The other interpretation, however, is to divide imaginary time into
$\beta/\Delta\tau$ slices and treat the system as having $m\beta/\Delta\tau$
replicas, with the kinetic term simply as an interaction between them. With
this interpretation, it is natural that $\alpha_{ab}$ should be promoted to
$\alpha_{ab}(\tau_{1},\tau_{2})$ but otherwise be treated as before. We can
then easily write down how the classical free energy should be modified (and
we have confirmed that a more mechanical derivation gives the same result).

The term $f_{\textrm{position}}$ is the same as in the classical case, and the
term $f_{\textrm{molecular}}$ simply gains a contribution coming from non-zero
frequencies (note that the zero-frequency contribution is unmodified since we
have fixed $G_{0}=0$):

$\beta f_{\textrm{molecular}}=-\frac{m-1}{2}\log{\frac{\pi
e\ell^{2}A}{d^{2}}}-\frac{m}{2}\sum_{\omega\neq 0}\log{G_{\omega}}.$ (24)

The interaction term becomes

$\displaystyle\beta f_{\textrm{virial}}$
$\displaystyle=-\frac{\hat{\phi}}{2}\int\textrm{d}he^{h}\left[e^{\sum_{ab}\int\textrm{d}\tau_{1}\textrm{d}\tau_{2}\alpha_{ab}(\tau_{12})\frac{\partial^{2}}{\partial
h_{a}(\tau_{1})\partial
h_{b}(\tau_{2})}}\prod_{a}e^{-\int\textrm{d}\tau\overline{v}(h_{a}(\tau))}\bigg{|}_{h_{a}(\tau)=h+\alpha_{aa}(0)}-1\right]$
(25)
$\displaystyle=-\frac{\hat{\phi}}{2}\int\textrm{d}he^{h}\left[\left(e^{\int\textrm{d}\tau_{1}\textrm{d}\tau_{2}\frac{A+G(\tau_{12})}{2}\frac{\partial^{2}}{\partial
h(\tau_{1})\partial
h(\tau_{2})}}e^{-\int\textrm{d}\tau\overline{v}(h(\tau))}\bigg{|}_{h(\tau)=h+\frac{A+G(0)}{2}}\right)^{m}-1\right].$

Thus the replicated free energy of the quantum liquid, i.e., Eq. (6), is

$\displaystyle\beta f$ $\displaystyle=\frac{1}{2}\log\frac{d}{2\pi
e\ell^{2}}-\frac{m-1}{2}\log{\frac{\pi
e\ell^{2}A}{d^{2}}}+\frac{m}{2}\sum_{\omega\neq
0}\left(\frac{\hat{M}}{2}\omega^{2}G_{\omega}-\log{G_{\omega}}\right)$ (26)
$\displaystyle\qquad-\frac{\hat{\phi}}{2}\int\textrm{d}he^{h}\left[\left(e^{\int\textrm{d}\tau_{1}\textrm{d}\tau_{2}\frac{A+G(\tau_{12})}{2}\frac{\partial^{2}}{\partial
h(\tau_{1})\partial
h(\tau_{2})}}e^{-\int\textrm{d}\tau\overline{v}(h(\tau))}\bigg{|}_{h(\tau)=h+\frac{A+G(0)}{2}}\right)^{m}-1\right].$

We can rewrite Eq. (26) using some convenient notation. Let
$\langle\,\cdot\,\rangle_{r,GA}$ denote the average over a periodic Gaussian
process $r(\tau)$ with covariance matrix $\langle
r(\tau_{1})r(\tau_{2})\rangle_{r,GA}=A+G(\tau_{12})$. By Eq. (19), the term in
parentheses in the second line of Eq. (26) is
$\langle\exp[-\int\textrm{d}\tau\overline{v}(h+r(\tau))]\rangle_{r,GA}$
evaluated at $h+(A+G(0))/2$. Denote this quantity by $p_{GA}(h)$:

$p_{GA}(h)\equiv\Big{<}e^{-\int\textrm{d}\tau\overline{v}(h+r(\tau))}\Big{>}_{r,GA},$
(27)

so that we can write

$\beta f=\frac{1}{2}\log\frac{d}{2\pi e\ell^{2}}-\frac{m-1}{2}\log{\frac{\pi
e\ell^{2}A}{d^{2}}}+\frac{m}{2}\sum_{\omega\neq
0}\left(\frac{\hat{M}}{2}\omega^{2}G_{\omega}-\log{G_{\omega}}\right)-\frac{\hat{\phi}}{2}\int\textrm{d}he^{h-\frac{A+G(0)}{2}}\Big{[}p_{GA}(h)^{m}-1\Big{]}.\;\;$
(28)

An identity that will often be useful is that

$p_{GA}(h)=e^{\frac{A}{2}\frac{\textrm{d}^{2}}{\textrm{d}h^{2}}}p_{G0}(h)=\int\frac{\textrm{d}z}{\sqrt{2\pi
A}}e^{-\frac{z^{2}}{2A}}p_{G0}(h+z).$ (29)

###  2.3 Saddle points

Our goal is to extremize Eq. (28) in the limit $m\rightarrow 1$. As described
in Sec. 1.2, the FPP $\Phi$ in particular is the first-order coefficient of
the free energy: $f\sim f_{0}+(m-1)\Phi$. However, care must be taken in
performing the extremization. If the zeroth-order term $f_{0}$ has non-trivial
dependence on a parameter, then extremizing $f_{0}$ alone determines that
parameter to leading order. Yet if $f_{0}$ is independent of it, then
extremizing $\Phi$ instead determines the parameter to leading order. In our
case, we have that (neglecting overall constants)

$\displaystyle\beta f_{0}$ $\displaystyle=\frac{1}{2}\sum_{\omega\neq
0}\left(\frac{\hat{M}}{2}\omega^{2}G_{\omega}-\log{G_{\omega}}\right)-\frac{\hat{\phi}}{2}\int\textrm{d}he^{h-\frac{G(0)}{2}}\Big{[}p_{G0}(h)-1\Big{]},$
(30) $\displaystyle\beta\Phi$
$\displaystyle=-\frac{1}{2}\log{A}+\frac{1}{2}\sum_{\omega\neq
0}\left(\frac{\hat{M}}{2}\omega^{2}G_{\omega}-\log{G_{\omega}}\right)-\frac{\hat{\phi}}{2}\int\textrm{d}he^{h-\frac{A+G(0)}{2}}p_{GA}(h)\log{p_{GA}(h)}.$
(31)

Note that $f_{0}$ is independent of $A$. Thus we should extremize $f_{0}$ with
respect to $G_{\omega}$ and extremize $\Phi$ with respect to $A$.

Setting $\partial f_{0}/\partial G_{\omega}=0$ gives the following saddle-
point equations:

$\displaystyle 0$
$\displaystyle=\frac{\hat{M}\omega^{2}}{4}-\frac{1}{2G_{\omega}}+\frac{\hat{\phi}}{4\beta}\int\textrm{d}he^{h-\frac{G(0)}{2}}\Big{[}p_{G0}(h)-1\Big{]}$
(32)
$\displaystyle\qquad+\frac{\hat{\phi}}{4G_{\omega}}\int\textrm{d}he^{h-\frac{G(0)}{2}}p_{G0}(h)-\frac{\hat{\phi}}{4G_{\omega}^{2}}\int\textrm{d}he^{h-\frac{G(0)}{2}}\Big{<}r_{\omega}r_{-\omega}e^{-\int\textrm{d}\tau\overline{v}(h+r(\tau))}\Big{>}_{r,G0}.$

While intimidating, Eq. (32) can be solved numerically by iteration as
described in Sec. 4.

Once Eq. (32) is solved, setting $\partial\Phi/\partial A=0$ gives the
following saddle-point equation for $A$:

$0=-\frac{1}{2A}+\frac{\hat{\phi}}{4}\int\textrm{d}he^{h-\frac{A+G(0)}{2}}p_{GA}(h)\log{p_{GA}(h)}-\frac{\hat{\phi}}{4}\int\textrm{d}he^{h-\frac{A+G(0)}{2}}p_{GA}^{\prime\prime}(h)\Big{[}\log{p_{GA}(h)}+1\Big{]}.$
(33)

We can integrate by parts to write this more compactly as

$\frac{1}{A}=\frac{\hat{\phi}}{2}\int\textrm{d}he^{h-\frac{A+G(0)}{2}}\frac{p_{GA}^{\prime}(h)^{2}}{p_{GA}(h)}.$
(34)

Keeping in mind that $G_{\omega}$ has already been determined, Eq. (34) is a
single equation for $A$ that can simply be plotted to determine graphically
whether solutions exist.

##  3 High Temperatures

Here we at last specialize to hard spheres, corresponding to
$\overline{v}(h)=\infty$ for $h<0$ and $\overline{v}(h)=0$ for $h>0$. This
gives $p_{GA}(h)$ as expressed in Eq. (27) a particularly simple
interpretation: it is the probability that the Gaussian process $r(\tau)$
never crosses below $-h$, which is equivalently the probability that the
minimum of the process is greater than $-h$. By symmetry, $p_{GA}(h)$ is also
the probability that the maximum is less than $h$, and therefore
$p_{GA}^{\prime}(h)$ is the probability density function for
$\max_{\tau}r(\tau)$.

As a warm-up, consider the saddle-point equations at infinite temperature,
which corresponds to the classical limit. Then $G_{\omega}=0$, there is no
$\tau$-dependence to account for, and $r$ is simply a Gaussian with variance
$A$:

$p_{0A}(h)=\int_{-h}^{\infty}\frac{\textrm{d}r}{\sqrt{2\pi
A}}e^{-\frac{r^{2}}{2A}}.$ (35)

We insert this expression into Eq. (34). A graphical inspection as in Fig. 2
demonstrates that solutions exist for $\hat{\phi}>\hat{\phi}_{c}\approx
4.8067$. This is the glass phase of the system.

Figure 2: Eq. (34) (left-hand minus right-hand side) as a function of $A$ for
different values of $\hat{\phi}$, specifically in the infinite-temperature
(i.e., classical) limit. At the critical value $\hat{\phi}_{c}\approx 4.8067$,
the curve acquires a zero near $A\approx 1.15$, indicating that the system is
in a glass phase.

At large but finite temperatures, $r(\tau)$ is no longer constant. Since
fluctuations in time can only increase the probability that
$\max_{\tau}r(\tau)$ exceeds some value, $p_{GA}(h)<p_{0A}(h)$. We will show
that $p_{GA}(h)$ can be well-approximated by $p_{0A}(h-\delta R)$ at high
temperatures for some $\delta R>0$. Our strategy is to determine $p_{G0}(h)$
and then convolve with a Gaussian of variance $A$ to obtain $p_{GA}(h)$ (see
Eq. (29)).

Figure 3: Comparison of $\hat{\phi}_{c}(\beta)$ estimated analytically at high
temperatures (orange curve — see Sec. 3) to that obtained numerically (blue
curve — see Sec. 4).

The first thing to note is that at high temperatures,

$G_{\omega}\sim\frac{2}{\hat{M}\omega^{2}},$ (36)

keeping in mind that $G_{0}=0$. First and foremost, this is confirmed by a
full numerical solution of the saddle-point equations as discussed in Sec. 4 —
Fig. 3 shows that the critical density $\hat{\phi}_{c}(\beta)$ obtained by
assuming Eq. (36) matches the numerical result as $\beta\rightarrow 0$. It can
also be justified by returning to Eq. (30) for $f_{0}$ — whereas each of the
first two terms has a derivative going as $O(\beta^{-2})$, the last term is
exponential in $G$ (recall that $p_{G0}(h)$ can be written as an exponential
of second derivatives) and so its derivative is at most only $O(\beta^{-1})$
(the factor of $\beta^{-1}$ because
$G(\tau)=\beta^{-1}\sum_{\omega}e^{-i\omega\tau}G_{\omega}$).

In any case, note that Eq. (36) is the covariance for Brownian motion, i.e.,
the action governing the distribution of $r(\tau)$ is simply
$\frac{\hat{M}}{4}\int\textrm{d}\tau|\frac{\textrm{d}r(\tau)}{\textrm{d}\tau}|^{2}$,
just as in Brownian motion. However, there are two features that make
$r(\tau)$ not literal Brownian motion: $r(\tau)$ is required to be periodic,
and the fact that $G_{0}=0$ means that $\int\textrm{d}\tau r(\tau)$ is fixed
to be zero.

Evaluating $p_{G0}(h)$ exactly for such a process is quite difficult, but we
can approximate it at small $\beta$. Note that when $G_{\omega}=0$,
$p_{00}(h)$ is simply a step function — 0 for $h<0$ and 1 for $h>0$. At small
but non-zero $G_{\omega}$, such as for small $\beta$, $p_{G0}(h)$ is still
(identically) zero for $h<0$ and still approaches 1 as $h$ increases. We will
approximate it as a step function that transitions at $h=\delta R$, where the
latter can be computed from $p_{G0}(h)$ by $\delta
R=\int_{0}^{\infty}\textrm{d}h(1-p_{G0}(h))$. One can show that this gives the
same answer to leading order when convolved against a slower-varying function,
such as we then do to calculate $p_{GA}(h)$.

Evaluating $\int_{0}^{\infty}\textrm{d}h(1-p_{G0}(h))$ is in fact much more
feasible. Integrate by parts to obtain
$\int_{0}^{\infty}\textrm{d}hhp_{G0}^{\prime}(h)$, and since
$p_{G0}^{\prime}(h)$ is the distribution for $\max_{\tau}r(\tau)$ as noted
above, this is simply the expected value of $\max_{\tau}r(\tau)\equiv
r_{\textrm{max}}$. The condition that $\int\textrm{d}\tau r(\tau)=0$ is
inconvenient, but we can replace it with the condition that $r(0)=0$ for the
purpose of calculating $\langle r_{\textrm{max}}\rangle_{r,G0}$: first note
that $\langle
r_{\textrm{max}}\rangle_{r,G0}=\langle(r_{\textrm{max}}-r(0))\rangle_{r,G0}$,
and then by translation invariance we can shift the trajectories being
averaged over to each have $r(0)=0$ (a more formal way of stating this is that
each trajectory with $\int\textrm{d}\tau r(\tau)=0$ is in one-to-one
correspondence with a trajectory having $r(0)=0$ and the same probability
under the Brownian measure). Thus we are left only with evaluating the average
of $r_{\textrm{max}}$ under periodic Brownian motion starting and ending at
$r(0)=0$ (sometimes called a “Brownian bridge”).

Figure 4: An example of a trajectory drawn from a period-1 Brownian bridge
(green) and its reflected counterpart used to calculate the average of
$r_{\textrm{max}}$ (red). At the first time where the green trajectory crosses
$r=0.2$ (middle horizontal line), we reflect about $r=0.2$ — while the green
trajectory ends at $r=0$ (lower line), the red trajectory ends at $r=0.4$
(upper line).

We can compute the distribution of $r_{\textrm{max}}$ under a Brownian bridge
using the method of images. The bridge is equivalently Brownian motion
conditioned on the event $r(0)=r(\beta)$. For any path in this class that
first crosses $h$ at some time $\tau_{0}$ (thus going from 0 to $h$ to 0 as
time increases from 0 to $\tau_{0}$ to $\beta$), there is an equally likely
path that is mirrored after $\tau_{0}$, ending at $2h$ (see Fig. 4).
Similarly, any path in Brownian motion that ends at $2h$ can be mapped onto a
periodic path by reflecting at the first time where the path crosses $h$. Thus
the probability of a trajectory crossing $h$ and returning to zero equals the
probability of a trajectory ending at $2h$, and therefore the cumulative
distribution function for $r_{\textrm{max}}$ (under the Brownian bridge) is
$1-\exp[-\hat{M}h^{2}/\beta]$. The expected value of $r_{\textrm{max}}$ is
then

$\delta
R=\int_{0}^{\infty}\textrm{d}h\frac{2\hat{M}h^{2}}{\beta}e^{-\frac{\hat{M}h^{2}}{\beta}}=\sqrt{\frac{\pi\beta}{4\hat{M}}}.$
(37)

As discussed above, we take $p_{G0}(h)$ to be a step function at $\delta R$.
Thus

$p_{GA}(h)\sim\int_{-h+\delta R}^{\infty}\frac{\textrm{d}z}{\sqrt{2\pi
A}}e^{-\frac{z^{2}}{2A}}=p_{0A}(h-\delta R).$ (38)

As a result, the saddle-point equation for $A$ amounts to

$\frac{1}{A}=\frac{\hat{\phi}}{2}\int\textrm{d}he^{h-\frac{A}{2}}\frac{p_{0A}^{\prime}(h-\delta
R)^{2}}{p_{0A}(h-\delta R)}=\frac{\hat{\phi}e^{\delta
R}}{2}\int\textrm{d}he^{h-\frac{A}{2}}\frac{p_{0A}^{\prime}(h)^{2}}{p_{0A}(h)}.$
(39)

We see that, compared to the classical case, the effect of quantum
fluctuations is simply to increase the filling fraction by an amount
$e^{\delta R}=1+O(\beta^{1/2})$. Thus the glass phase occurs at a slightly
smaller density, namely

$\hat{\phi}_{c}(\beta)\sim\left(1-\sqrt{\frac{\pi\beta}{4\hat{M}}}\right)\hat{\phi}_{c}(0).$
(40)

It is interesting to note that this correction, although small, is non-
analytic in $\beta$ and in particular is larger than what the first analytic
correction would be.

##  4 Numerical Results

The most demanding aspect of determining the phase diagram is solving Eq. (32)
for $G_{\omega}$, but even this can be done numerically as follows. Define the
“self-energy”

$\Sigma_{\omega}\equiv-\hat{\phi}\frac{\partial}{\partial
G_{\omega}}\int\textrm{d}he^{h-\frac{G(0)}{2}}\Big{[}p_{G0}(h)-1\Big{]},$ (41)

which corresponds to (twice) the last three terms of Eq. (32). The “solution”
for $G_{\omega}$ is then simply

$G_{\omega}=\left[\frac{\hat{M}\omega^{2}}{2}+\Sigma_{\omega}\right]^{-1}.$
(42)

We solve the saddle-point equations iteratively: use an initial guess for
$G_{\omega}$ to evaluate $\Sigma_{\omega}$ via Eq. (41), then evaluate an
updated $G_{\omega}$ via Eq. (42), which gives an updated $\Sigma_{\omega}$
via Eq. (41), and so on. While there is no guarantee that such a
straightforward iteration must converge, we find that in practice it works
quite well, especially when the result for $G_{\omega}$ at slightly higher
temperature and/or lower filling fraction is used as an initial guess.

Fig. 5 plots $G(\tau)$ at a rather low temperature for various values of
$\hat{\phi}$ (still specializing to hard spheres). One can confirm that at
$\hat{\phi}=0$, the Fourier transform of $G_{\omega}$ is quadratic in $\tau$.
We see in Fig. 5 that $G(\tau)$ begins to better resemble an exponential
instead as the filling fraction increases.

Of particular importance for the saddle-point equations is the extreme-value
statistics of a Gaussian process with covariance $G(\tau_{12})$. When
$G(\tau)$ decays exponentially over an interval $\gamma^{-1}$ which is much
smaller than $\beta$, this is essentially the extreme-value statistics of
$\gamma\beta$ independent Gaussians with variance $G(0)$, which is known to
follow the Gumbel distribution:
$\textrm{Pr}[r_{\textrm{max}}<x]=\exp[-e^{-(x-\mu)/\sigma}]$ for certain $\mu$
and $\sigma$. Fig. 6 shows the distribution of $r_{\textrm{max}}$ for
$\hat{\phi}=0$ and $\hat{\phi}=100$ alongside best-fit Gumbel distributions.
Although there are systematic deviations, we see that the Gumbel distribution
gives a reasonably good fit in both cases (even for $\hat{\phi}=0$, where
$G(\tau)$ is far from an exponential).

Most importantly, once $G(\tau)$ is computed for a given choice of $\beta$ and
$\hat{\phi}$, we can determine whether any finite value of $A$ solves Eq.
(34). If so, then the system is in the glass phase. If not, then the system is
in the liquid phase. The resulting phase diagram of quantum hard spheres is
given in Fig. 1. First recall that the numerical results agree quite well with
those of the high-temperature expansion at small $\beta$ (see Fig. 3). More
broadly, we find that the phase boundary $\hat{\phi}_{c}(\beta)$ continues to
decrease monotonically at larger $\beta$. This is in contrast with the results
of MCT [31], but we again stress that these conclusions are not inconsistent,
since MCT relies on a different set of approximations than the high-
dimensional limit considered here.

Figure 5: $G(\tau)$ for various values of $\hat{\phi}$ at $\beta/\hat{M}=10$.

Figure 6: Distribution of $r_{\textrm{max}}$ for two choices of $\hat{\phi}$,
representative of small and large filling fractions, both at
$\beta/\hat{M}=10$. Orange curves are best-fit Gumbel distributions. (Left)
Small $\hat{\phi}$, namely the non-interacting limit. (Right) Large
$\hat{\phi}$.

##  5 Conclusion

To summarize, in this work we have calculated the phase boundary between the
liquid and glass phases of quantum hard spheres, specifically in the infinite-
dimensional limit, where a formally exact calculation of the Franz-Parisi
potential is possible. We have determined the leading-order correction to the
critical filling fraction $\hat{\phi}_{c}$ analytically at small inverse
temperature $\beta$, finding that quantum fluctuations give the spheres an
effective radius that is enhanced by a multiple of the thermal de Broglie
wavelength (and thus $\hat{\phi}_{c}$ decreases as $\beta$ increases). We have
also computed the phase boundary numerically at lower temperatures, confirming
that $\hat{\phi}_{c}$ continues to decrease monotonically.

While it is undeniably valuable to have controlled calculations of the phase
diagram in certain limits, such as we have given here, many open questions and
directions for future research remain. First and foremost is investigating the
extent to which our results are relevant for actual three-dimensional quantum
glasses. As we have noted, the phase diagram of infinite-dimensional quantum
hard spheres stands in some contrast to the predictions of MCT [31]. Neither
calculation can claim to give a fully accurate description of physical
glasses, however — for one thing, the lifetimes of metastable states are
finite outside of mean-field theory and so the liquid-glass phase boundary is
(at least in principle) ill-defined. It would be of great interest to
understand which vestiges of the phase diagram remain in realistic systems,
and in particular whether the critical filling fraction is indeed monotonic.

Another open question is how the phase diagram is modified for soft-sphere
interaction potentials (i.e., potentials that remain finite for separations
less than $\ell$). While we expect it to be qualitatively similar to a certain
degree, one new ingredient that soft-sphere potentials add is quantum
tunneling through the energy barriers repelling spheres. One might imagine
that this effect suppresses glassiness by allowing the system to better escape
from metastable states, and thus there may be an interesting competition
between tunneling and the increased effective radius due to quantum
fluctuations.

Lastly, it will be important to understand how the tendency to form a glass
competes with other types of order at low temperature — this question is
relevant for classical systems as well, but it is especially important here
since our calculation suggests that any (hard-sphere) liquid will become a
glass at sufficiently low temperature. Many substances in nature of course do
not form glasses, instead entering phases such as solids or superfluids.
Accounting for such effects is beyond the scope of the calculation presented
here, but if a suitable generalization can be developed, then it would be
quite interesting to study the interplay between these different types of
order.

### Acknowledgements

We would like to thank L. Foini, D. Reichman, P. Urbani, and F. Zamponi for
valuable discussions and encouragement. This material is based upon work
supported by the JQI fellowship (M.W.) and the National Science Foundation
Graduate Research Fellowship Program under Grant No. DGE 2236417 (R.B.). V.G.
was supported by the U.S. Department of Energy, Office of Science, Basic
Energy Sciences under Award No. DE-SC0001911. We gratefully acknowledge
support from Joint Quantum Institute and from AFOSR under grant number
FA9550-19-1-0360. The numerical data were obtained using the Zaratan high-
performance computing cluster at the University of Maryland, College Park.

##  Appendix A Life at large $d$

The volume of a radius-$\ell/2$ sphere in $d$ dimensions is given by

$V_{d}=\frac{\pi^{\frac{d}{2}}}{\Gamma\big{(}\frac{d}{2}+1\big{)}}\left(\frac{\ell}{2}\right)^{d}\sim\frac{1}{\sqrt{\pi
d}}\left(\frac{\pi e\ell^{2}}{2d}\right)^{\frac{d}{2}},$ (43)

where $\Gamma$ denotes the gamma function. Note that the volume shrinks even
faster than exponentially, in particular faster than $\ell^{d}$. This means
that in large dimensions, the most obvious packing of spheres — a simple cubic
lattice — only fills an infinitesimal fraction of space.

More efficient packings exist with far larger filling fractions — this is in
fact a well-studied problem [43, 44, 45]. Hermann Minkowski proved the first
real lower bound on the maximal filling fraction, that in any dimension a
fraction of $2^{-d}$ is possible. The proof is elementary: given a packing, if
we double the radius of each sphere ($r\rightarrow 2r$), then either every
point in space becomes covered by an expanded sphere or there is a point whose
distance from the center of every sphere is at least $2r$. In the former case,
the original filling fraction must be at least $2^{-d}$. In the latter case,
the packing must not have been optimal (because a sphere of radius $r$ can be
inserted at the aforementioned point).

More advanced proofs have established that packings exist with filling
fractions $\phi$ at various multiples of $d2^{-d}$. Thus it is natural to
write $\phi=d2^{-d}\hat{\phi}$ and consider $\hat{\phi}$ finite at large $d$.
It is interesting to note that the analysis in this paper can apply to
arbitrarily large $\hat{\phi}$, even though packings at sufficiently large
$\hat{\phi}$ haven’t been rigorously proven to exist (although they have been
established for $\hat{\phi}$ well beyond where we locate the liquid-glass
phase boundary [45]).

The search for upper bounds on allowed filling fractions has been even more
difficult. The best known bounds scale as $2^{-cd}$ for a constant $c\approx
0.599$, with hopes to possibly improve to $c\approx 0.779$, but anything
further would require substantial breakthroughs [46].

Another consequence of the scaling of $V_{d}$ in Eq. (43) is that we can
ignore the effects of particle exchange and in particular the possibility of
Bose-Einstein condensation. The densities that we consider in this paper are
of order $d/2^{d}V_{d}$, which grows predominantly as $d^{d/2}$. On the other
hand, the density scale for condensation is of order
$\lambda_{T}^{-d}\sim(M/\beta)^{d/2}$, where $\lambda_{T}$ denotes the thermal
de Broglie wavelength. Since we take $M$ to be of order $d^{2}$, the
condensation density grows as $d^{d}$, much larger than the critical density
for glass formation.

It is also worth asking how many spheres lie within a distance
$\ell(1+\delta)$ from the center of a given sphere (note that two spheres are
in contact when they are separated by distance $\ell$). Multiplying the
density by the volume of a radius-$\ell$ sphere (not $\ell/2$), we find that
the number scales as $de^{d\delta}$. Thus the number of spheres whose edges
are within $O(1/d)$ is $O(d)$, i.e., large. This motivates the interpretation
of the infinite-dimensional liquid as a mean-field system — in general,
systems in which each degree of freedom interacts with a large number of
neighbors are well-described by mean-field theory.

## References

  * [1] L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, and W. van Saarloos, Dynamical Heterogeneities in Glasses, Colloids, and Granular Media.  Oxford University Press, 07 2011.
  * [2] K. Binder and W. Kob, Glassy Materials and Disordered Solids: An Introduction to Their Statistical Mechanics.  World Scientific Publishing Company, 2005.
  * [3] W. Götze, Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory.  Oxford University Press, 2008.
  * [4] V. Lubchenko and P. G. Wolynes, “Theory of structural glasses and supercooled liquids,” Annual Review of Physical Chemistry, vol. 58, no. 1, pp. 235–266, 2007.
  * [5] F. Morone, F. Caltagirone, E. Harrison, and G. Parisi, “Replica theory and spin glasses,” 2014.
  * [6] T. Castellani and A. Cavagna, “Spin-glass theory for pedestrians,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2005, no. 05, p. P05012, 2005.
  * [7] H. Yoshino, “Replica theory of the rigidity of structural glasses,” The Journal of Chemical Physics, vol. 136, no. 21, 2012.
  * [8] A. Altieri and M. Baity-Jesi, “An introduction to the theory of spin glasses,” 2023.
  * [9] M. Mezard, G. Parisi, and M. Virasoro, Spin Glass Theory and Beyond: An Introduction to the Replica Method and Its Applications.  World Scientific Lecture Notes in Physics, World Scientific Publishing Company, 1987.
  * [10] P. Charbonneau, E. Marinari, G. Parisi, F. Ricci-Tersenghi, G. Sicuro, and F. Zamponi, eds., Spin Glass Theory and Far Beyond - Replica Symmetry Breaking After 40 Years.  World Scientific Publishing Company, 2023.
  * [11] D. R. Reichman and P. Charbonneau, “Mode-coupling theory,” Journal of Statistical Mechanics: Theory and Experiment, 2005.
  * [12] G. Szamel, “Mode-coupling theory and beyond: A diagrammatic approach,” Progress of Theoretical and Experimental Physics, 2012.
  * [13] W. Götze, Complex Dynamics of Glass-Forming Liquids.  OUP Oxford, 2009.
  * [14] S. P. Das, “Mode-coupling theory and the glass transition in supercooled liquids,” Rev. Mod. Phys., vol. 76, pp. 785–851, 2004.
  * [15] T. E. Markland, J. A. Morrone, B. J. Berne, K. Miyazaki, E. Rabani, and D. R. Reichman, “Quantum fluctuations can promote or inhibit glass formation,” Nature Physics, vol. 7, no. 2, pp. 134–137, 2011.
  * [16] C. Liu, G. Biroli, D. R. Reichman, and G. Szamel, “Dynamics of liquids in the large-dimensional limit,” Physical Review E, vol. 104, no. 5, 2021.
  * [17] M. Cardenas, S. Franz, and G. Parisi, “Constrained Boltzmann-Gibbs measures and effective potential for glasses in hypernetted chain approximation and numerical simulations,” The Journal of Chemical Physics, vol. 110, no. 3, pp. 1726–1734, 1999.
  * [18] M. Sellitto and F. Zamponi, “A thermodynamic description of colloidal glasses,” EPL, vol. 103, no. 4, p. 46005, 2013.
  * [19] T. Maimbourg, J. Kurchan, and F. Zamponi, “Solution of the dynamics of liquids in the large-dimensional limit,” Phys. Rev. Lett., vol. 116, p. 015902, 2016.
  * [20] E. Agoritsas, T. Maimbourg, and F. Zamponi, “Out-of-equilibrium dynamical equations of infinite-dimensional particle systems i. the isotropic case,” Journal of Physics A: Mathematical and Theoretical, vol. 52, no. 14, p. 144002, 2019.
  * [21] M. Mezard and G. Parisi, “A tentative replica study of the glass transition,” Journal of Physics A: Mathematical and General, vol. 29, no. 20, pp. 6515–6524, 1996.
  * [22] G. Parisi and F. Zamponi, “Mean-field theory of hard sphere glasses and jamming,” Reviews of Modern Physics, vol. 82, no. 1, pp. 789–845, 2010\.
  * [23] F. Zamponi, “Theory of simple glasses,” arXiv:1212.0390, 2012.
  * [24] L. Berthier, H. Jacquin, and F. Zamponi, “Microscopic theory of the jamming transition of harmonic spheres.,” Physical review. E, Statistical, nonlinear, and soft matter physics, vol. 84 5 Pt 1, p. 051103, 2011.
  * [25] L. Foini, G. Semerjian, and F. Zamponi, “Quantum Biroli-Mézard model: Glass transition and superfluidity in a quantum lattice glass model,” Phys. Rev. B, vol. 83, p. 094513, 2011.
  * [26] S. Franz, T. Maimbourg, G. Parisi, and A. Scardicchio, “Impact of jamming criticality on low-temperature anomalies in structural glasses,” Proceedings of the National Academy of Sciences, vol. 116, no. 28, pp. 13768–13773, 2019.
  * [27] C. Artiaco, F. Balducci, G. Parisi, and A. Scardicchio, “Quantum jamming: Critical properties of a quantum mechanical perceptron,” Physical Review A, vol. 103, no. 4, 2021.
  * [28] P. Urbani, “Quantum exploration of high-dimensional canyon landscapes,” arXiv:2301.11236, 2023.
  * [29] L. F. Cugliandolo and M. Mueller, “Quantum glasses – a review,” arXiv:2208.05417, 2022.
  * [30] E. B. Rozenbaum, L. A. Bunimovich, and V. Galitski, “Early-time exponential instabilities in nonchaotic quantum systems,” Physical Review Letters, vol. 125, no. 1, 2020.
  * [31] T. Markland, J. Morrone, B. Berne, K. Miyazaki, E. Rabani, and D. Reichman, “Quantum fluctuations can promote or inhibit glass formation,” Nature Phys., vol. 7, 2010.
  * [32] G. Parisi and F. Zamponi, “Mean-field theory of hard sphere glasses and jamming,” Rev. Mod. Phys., vol. 82, pp. 789–845, 2010.
  * [33] G. Parisi, P. Urbani, and F. Zamponi, Theory of Simple Glasses: Exact Solutions in Infinite Dimensions.  Cambridge University Press, 2020.
  * [34] G. Parisi and F. Slanina, “Toy model for the mean-field theory of hard-sphere liquids,” Physical Review E, vol. 62, no. 5, pp. 6554–6559, 2000.
  * [35] R. Mari and J. Kurchan, “Dynamical transition of glasses: From exact to approximate,” The Journal of Chemical Physics, vol. 135, no. 12, 2011.
  * [36] S. Franz and G. Parisi, “Recipes for metastable states in spin glasses,” Journal de Physique I, vol. 5, no. 11, pp. 1401–1415, 1995.
  * [37] M. Mezard and A. Montanari, Information, Physics, and Computation.  Oxford University Press, 2009.
  * [38] K. H. Fischer and J. A. Hertz, Spin Glasses.  Cambridge University Press, 1991.
  * [39] J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids.  Elsevier, 2006.
  * [40] M. Kardar, Statistical Physics of Particles.  Cambridge University Press, 2007.
  * [41] F. Schreiber, F. Zanini, and F. Roosen-Runge, “Virial expansion – a brief introduction,” 2011.
  * [42] J. Mayer and M. Mayer, Statistical Mechanics.  J. Wiley, 1966.
  * [43] C. A. Rogers, Packing and Covering.  Cambridge Tracts in Mathematics and Mathematical Physics, Cambridge University Press, 1964.
  * [44] C. A. Rogers, “Existence theorems in the geometry of numbers,” Annals of Mathematics, vol. 48, no. 4, pp. 994–1002, 1947.
  * [45] A. Venkatesh, “A Note on Sphere Packings in High Dimension,” International Mathematics Research Notices, vol. 2013, no. 7, pp. 1628–1642, 2012.
  * [46] H. Cohn and Y. Zhao, “Sphere packing bounds via spherical codes,” Duke Mathematical Journal, vol. 163, no. 10, 2014.