\section*{Results}
\emph{Phase diagram --} 
Using an approach similar to that used for solving the SK model, the exact $d\rightarrow} \def\la{\left\langle} \def\ra{\right\rangle\io$ solution for 
$d$-dimensional identical hard spheres of unit diameter
can be formulated in terms of a caging order parameter $\D(y)$~\cite{Pa83,MPV87}. 
This functional order parameter, which encodes the width $\D$ of metabasins on a (properly defined) scale $y$, 
is obtained by numerically solving a set of integro-differential equations (Supplementary Note~1)
and then used to calculate the theoretical liquid-glass phase diagram (Fig.~\ref{fig:PD}). 
The theory predicts that a compressed liquid falls out of equilibrium and becomes a glass at a pressure that
depends on the compression rate. Once in a glass state, further compression results in a quick increase of the system pressure $p$, and upon jamming $p\rightarrow} \def\la{\left\langle} \def\ra{\right\rangle \io$~\cite{PZ10}. The final jamming density depends
on compression speed, hence defining a J-line of jammed states~\cite{PZ10}. Two examples of glass compression obtained using an approximate state following are reported in Fig.~\ref{fig:PD}~\cite{KZ10} . 

Independently of compression rate, the glass basin in which the system is initially trapped
undergoes a Gardner transition~\cite{Ga85}, at a line computed in Ref.~\cite{KPUZ13}.
Our key result is that, at pressures above this line, \emph{basins transform
into metabasins that contain a collection of marginally stable glasses}, 
a phenomenon that is described by a non-trivial caging order parameter $\D(y)$ as in the SK model~\cite{Pa83}.
Finding the solution that describes the marginal phase allows us to delimit the marginal phase boundary (Supplementary Note~1) to within the Gardner transition
line of Ref.~\cite{KPUZ13}, the J-line, and the ``threshold'' line determined following the prescription of Ref.~\cite{Ri13}. 
The fact that within this region at least one eigenvalue of the stability matrix in the free energy space vanishes
confirms that this phase is indeed marginally stable (Supplementary Note~1)~\cite{DK83}. We also find that,
while the radius of the innermost fractal basins shrinks to zero as a power-law $\D_{\mathrm{EA}} \sim p^{-\k}$ (see below), 
the radius of the largest metabasins remains of order one.  
Close to jamming the total entropy of a group of metabasins of width $\Delta$ grows as $\Delta^{1/\k}$, 
 hence the basins have a phase space structure whose fractal dimension is $2/\k$ 
(see Supplementary Note~1 for a more detailed discussion).
The {\it marginal} phase is thus also {\it fractal}.

The existence of the marginal phase can be qualitatively tested by molecular dynamics (MD) numerical simulations
in finite $d$ (Supplementary Note~2) by considering the outcome of a slow compression from the liquid up to 
jamming~\cite{CIPZ11,CCPZ12}.
Jammed systems are isostatic, 
and thus particles have an average of 2$d$ force-bearing neighbors~\cite{Mo98,TW99,Ro00,OLLN02,BW09b,CCPZ12}, 
which is much smaller than the $\mathcal{O}(e^d)$ neighbors that isotropically cage a particle in an equilibrated dense liquid. 
Because the identity of the force-bearing neighbors at jamming uniquely characterizes the state,
their emergence sensitively depends on the landscape structure (Fig.~\ref{fig:basins}). 
In the simple basin scenario, force-bearing neighbors at jamming should be fully determined immediately upon leaving the equilibrium liquid; in the meta/subbasin scenario, that determination should only occur once sufficiently deep in the glass for transitions between subbasins to be fully suppressed; in a fractal phase, 
by contrast, the contacts should be gradually determined as jamming is approached. 
To test this scenario, we consider a glass configuration at pressure $p_{\mathrm{init}}$. Starting from this configuration, we perform several
independent compressions up to $p_\mathrm{f} = 10^{10}$ and for each compressed configuration we measure the force network. 
We obtain a set of contact variables $f^{(a)}_{ij}$, which are set to unity if particles $i$ and $j$ 
form a force-bearing contact in configuration $a$ and
to zero otherwise. The average of $\langle f^{(a)}_{ij} f^{(b)}_{ij} \rangle$ over pairs $ab$ of compressed configurations and over contacts $ij$ provides
a measure of similarity between the force networks. The fact that this quantity 
increases smoothly upon increasing $p_\mathrm{init}$ indicates that
the force network is only partially encoded in the initial configuration, in support
of the fractal landscape scenario (Fig.~\ref{fig:MDSA}a and Supplementary Note~2). 

\begin{figure}
\center{
\includegraphics[width=0.5\columnwidth]{MSDA.eps}}
\caption{
(a) Overlap between the force network edges $f_{ij}={0,1}$
that connect two particles $i$ and $j$ in two glass configurations $a$ and $b$ (at pressure $p_\mathrm{f} = 10^{10}$),
obtained by independent compression of the same initial configuration at initial pressure $p_\mathrm{init}$ (Supplementary Note~2). 
(b) Time-evolution of the mean-square displacement $\Delta(t)$ for glasses at $p=10^{2}$, $10^3$, $10^4$, $10^5$, $10^6$, $10^7$, and $10^{8}$ in $d$=4. The solid line indicates the ballistic $dt^2$ behavior. The long-time value is the cage size $d \, \Delta_{\mathrm{EA}}$.
(c) The pressure evolution of the cage size $\Delta_{\mathrm{EA}}$ in  various dimension closely follows a power-law $\sim p^{-\k}$ with 
$\k=1.41575$ as predicted by the theory.
(d) Analytical results for the order parameter $\D(y)$ at $\widehat} \newcommand{\wt}{\widetilde\f=10$.
Increasing the cutoff $y_{\rm max}\sim p$ indicates that the scaling regime $\D(y) \sim y^{-\k}$ extends to all
$y$. (Supplementary Note~1)}
\label{fig:MDSA}
\end{figure}	

\emph{Criticality of the jamming transition --} 
The equations that describe the marginal phase are formulated in terms of the caging order parameter $\Delta(y)$ and the pair correlation function $g(r)$, which also encodes the probability distribution of forces in the packing. 
Upon approaching the J-line, i.e., as $p\rightarrow} \def\la{\left\langle} \def\ra{\right\rangle\io$, these equations develop a scaling regime (Fig.~\ref{fig:MDSA}d) that is characterized by three main critical exponents: $\theta$ for the weak forces, $\alpha$ for the quasi-contacts, and $\kappa$ for $\D$ itself.
A (non-trivial) generalization of the approach developed for the SK model~\cite{Pa06} allows us to obtain theoretical values for these exponents (Supplementary Note~1). Interestingly, the condition that fixes their value is precisely equivalent to the marginal stability condition. 
The theory therefore predicts that the criticality of the jamming transition directly follows from its location inside the marginal phase.

A striking signature of marginality is the scaling of the innermost basin width captured by the Edwards--Anderson cage size 
$\D_{\mathrm{EA}} \sim p^{-\kappa}$. Although $\kappa=3/2$ was proposed in earlier studies~\cite{WNW05,WSNW05,BW09b,IBB13},
the theory predicts a slightly smaller $\kappa = 1.41574$ that is in remarkable agreement with our numerical results 
(Fig.~\ref{fig:MDSA}c).
Because single-particle caging by immediate neighbors (a simple Einstein model for glasses), 
would give $\k=2$~\cite{IBB13}, $\k < 2$ implies that fluctuations near jamming are divergently larger than for independent vibrations, in support of
their cooperative nature~\cite{BW09b,IBB13}. 
Note that if one ignores the fractal phase, an explicit computation erroneously gives $\k=1$~\cite{KPUZ13}.
Also, note that the exponent $\k$ controls the fractal dimension of the basins, as discussed above. 

The pair correlation function $g(r)$ also bears a signature of the criticality at the jamming transition.
The theory predicts, consistently with the analysis of~\cite{DTS05}, that when $p \rightarrow} \def\la{\left\langle} \def\ra{\right\rangle \io$, $g(r)$ develops an isostatic contact peak
characterized by a scaling function $\FF(\l) \equiv g(r)/g(1)$ for $\l=(r-1)p$ (Supplementary Note~1 and Fig.~\ref{fig:gr}). 
It also predicts that the scaling
function of the contact peak decays as $\FF(\l)\sim \l^{-2-\th}$ at large $\l$.
The distribution $P(f)$ of inter-particle forces in the packing, which is related to the scaling function of the contact peak by
$\FF(\l)=\int_0^\io df \, f\, P(f) \, e^{-\l f}$~\cite{DTS05,PZ10,CCPZ12}, thus also decays as a power law $P(f)\sim f^\theta$
at small forces. Note that, as observed in~\cite{Wy12}, this phenomenon is closely related to what happens to the distribution of  
frozen fields in the SK model~\cite{SD84}, 
which is thought to explain the Coulomb gap in interacting electron systems~\cite{MP07}.
Beyond the contact peak, the slower decay of pair correlation function follows another power-law $g(r) \sim (r-1)^{-\a}$ that describes the abundance
of quasi-contacts.
These scalings of $g(r)$ are crucial for determining the mechanical stability of packings~\cite{Wy12,LDW13}. Perturbing a packing breaks
some contacts with small forces, while also forming new contacts from what were previously quasi-contacts. Based on this observation a scaling relation for mechanical stability $\a = 1/(2+\th)$ can be derived~\cite{Wy12}. Remarkably, the exponents {\it predicted} by our theory, $\a=0.41269$ and $\th=0.42311$, satisfy this scaling relation to within numerical precision.
Prior estimates of these exponents were also obtained by numerical simulation. The quasi-contact exponent $\a$ has been measured by several groups in dimension $d$ ranging
from $2$ to $13$, all obtaining 
roughly $\a\approx 0.4$~\cite{DTS05,SDST06,CCPZ12,LDW13,AST13}, the most precise estimates being $\a=0.41(3)$~\cite{CCPZ12}.
The weak force exponent $\th$ is, however, more difficult to measure,
and values spanning the interval $\th \approx 0.2 \div 0.45$ have been reported~\cite{CCPZ12,LDW13}. 
Although the existence of a second exponent $\theta'<\theta$ has been shown to affect the tail of $P(f)$~\cite{LDW13}, 
its role in determining $\FF(\l)$ and its large-dimensional scaling remains to be clarified.
Additional numerical simulations are thus needed to test the theory more stringently. 


A prediction for the force distribution $P(f)$ at jamming is also available from the theory.
But because this function is not completely determined by the scaling regime, it must be obtained by solving
the full equations that describe the marginal phase. 
Numerically, the function $\FF(\l)=\int_0^\io df \, f\, P(f) \, e^{-\l f}$ is much easier to measure than $P(f)$ because it only depends on structural information, while, in
hard spheres, forces must be determined from the collision dynamics. A measure of $\FF(\l)$ thus also provides a more precise way to measure $\th$.
The theoretical prediction for the scaling function $\FF(\l)$ and therefore for $\th$ are
tested against numerical simulations in Fig.~\ref{fig:gr}, with very good agreement.

\begin{figure}
\center{
\includegraphics[width=0.5\columnwidth]{gr.eps}}
\caption{
Scaling of the contact peak of the pair correlation $g(r)$. The theory predicts that in the limit $p\rightarrow} \def\la{\left\langle} \def\ra{\right\rangle\io$ one has
$g(r)/g(1) = \FF(\l)$ with $\l = (r-1) p$. Points are numerical data for $g(r)$ in $d=4$ obtained at several pressures 
$p=10^{4},10^5, \cdots, 10^{12}$;
the full line is the theoretical prediction. 
}
\label{fig:gr}
\end{figure}	



\section*{Discussion} 
We have described the marginal phase that is present below
the Gardner transition for hard spheres in $d=\io$.
Using this result we have shown that the jamming transition happens inside the marginal phase and
its criticality in low-dimensional systems
is perfectly described by our approach.
This analysis opens the way for analytically determining many other properties of jammed packings,
such as their shear modulus~\cite{BW06,Yo12}, 
and the properties of avalanches~\cite{LMW10}. 
It also offers a pathway for understanding other glass properties. The nature of aging, for instance, undergoes an abrupt change
at a temperature lower than the glass transition~\cite{St77}, which could signal the crossing between the stable and the fractal glass regimes~\cite{Ri13}.
The observation of dynamical heterogeneities in low temperature glasses~\cite{VKBZ02,BDC08} at timescales much shorter than the inter-basin relaxation 
(a fact that does not fit in the two-timescale picture)  may also be taken as a  signature of
the fractal phase, wherein the dynamical correlation length diverges.
Finally, note that in the $d=\infty$ limit, the barriers controlling the inter-basin relaxation should scale as $d$, 
and the largest separating subbasins as $d^{1/3}$, suggesting the existence of quasi-localized excitations~\cite{XVLN10} 
(stringlike, in three dimensions~\cite{WWW}).



\section*{Methods}

Results are based on the combination of analytical and numerical methods. Analytical results come from the exact solution of hard spheres 
in the limit $d\rightarrow} \def\la{\left\langle} \def\ra{\right\rangle\io$, which, for convenience, is obtained 
using the replica method, but any other method would give the same result.
The fractal phase is described by a function $\D(y)$ for $y\in [1,y_{\rm max}]$, as in the SK model~\cite{Pa83}. The cutoff $y_{\rm max} \sim p$ diverges with pressure.
With these definitions $\D_{\rm EA} = \D(y_{\rm max})$ is the mean square displacement in the smallest subbasins, where
$\D(y)$ can be computed by numerically solving a set of coupled integro-differential equations obtained from the replica approach (Supplementary Note~1). Numerical results are obtained by standard event-driven 
molecular dynamics simulations in $d=3$ to 8~\cite{SDST06,CCPZ12}. Compressions are made using
the Lubachevsky-Stillinger algorithm~\cite{SDST06} (Supplementary Note~2)




\emph{Acknowledgments --}
We thank Silvio Franz for an extremely useful discussion. We also thank Ludovic Berthier, Carolina Brito, Florent Krzakala, Edan Lerner, Federico Ricci-Tersenghi, Tommaso Rizzo, Hajime Yoshino and Matthieu Wyart for many interesting exchanges related to this work. Financial support was provided by
the European Research Council through ERC grant agreement
no. 247328 and ERC grant NPRGGLASS. PC acknowledges support from the Alfred P. Sloan Foundation.

