Merge branch 'master' into apsaps.v3

6 files changed, 466 insertions, 203 deletions

diff --git a/appeal.tex b/appeal.tex
new file mode 100644
index 0000000..6af3109
--- /dev/null
+++ b/appeal.tex
@@ -0,0 +1,77 @@
+\documentclass[a4paper]{letter}
+
+\usepackage[utf8]{inputenc} % why not type "Bézout" with unicode?
+\usepackage[T1]{fontenc} % vector fonts plz
+\usepackage{newtxtext,newtxmath} % Times for PR
+\usepackage[
+  colorlinks=true,
+  urlcolor=purple,
+  linkcolor=black,
+  citecolor=black,
+  filecolor=black
+]{hyperref} % ref and cite links with pretty colors
+\usepackage{xcolor}
+\usepackage[style=phys]{biblatex}
+
+\addbibresource{bezout.bib}
+
+\signature{
+  \vspace{-6\medskipamount}
+  \smallskip
+  Jaron Kent-Dobias \& Jorge Kurchan
+}
+
+\address{
+  Laboratoire de Physique\\
+  Ecole Normale Sup\'erieure\\
+  24 rue Lhomond\\ 
+  75005 Paris
+}
+
+\begin{document}
+\begin{letter}{
+  Editorial Office\\
+  Physical Review Letters\\
+  1 Research Road\\
+  Ridge, NY 11961
+}
+
+\opening{To the editors of Physical Review,}
+
+We wish to appeal your decision on our manuscript \emph{Complex complex
+landscapes}, which received a single referee report.
+
+We believe that the referee's overall criticisms of our paper are not entirely
+justified (and above all, difficult to answer).  We have, however, submitted a
+revised manuscript clarifying the specific aspects that the referee found trying.
+
+The referee seems particularly worried that we have cited articles that are not
+themselves sufficiently cited, so we thought it may be useful propose a set of
+referees that are beyond suspicion of incompetence or uncitedness:
+
+\begin{tabular}{ll}
+  G Ben Arous  & Courant \\ 
+  M Berry      & Bristol\\
+  Y Fyodorov   & King's College London \\
+  Daniel Fisher& Stanford\\
+  T Lubensky   & U Penn\\
+  M Moore      & Manchester \\
+  E Witten     & IAS Princeton
+\end{tabular}
+
+We have also pointed out  some very first results of the geometric implications for a
+random complex landscape of our calculation, something we plan to expend on in
+a future full article.
+
+Let us conclude by remarking that, although it is probably true that this paper will not
+make more than one hundred citations next year, we are  confident that it will still 
+be considered relevant in ten years time.
+
+
+\closing{Sincerely,}
+
+\vspace{1em}
+
+\end{letter}
+
+\end{document}
diff --git a/bezout.tex b/bezout.tex
index 8e354da..1480cce 100644
--- a/bezout.tex
+++ b/bezout.tex
@@ -42,69 +42,91 @@
 
 \maketitle
 
-Spin-glasses have long been considered the paradigm of many variable `complex
-landscapes,' a subject that includes neural networks and optimization problems,
-most notably constraint satisfaction \cite{Mezard_2009_Information}.  The most tractable family of these
-are the mean-field spherical $p$-spin models \cite{Crisanti_1992_The} (for a
-review see \cite{Castellani_2005_Spin-glass}) defined by the energy
+Spin-glasses are the paradigm of many-variable `complex landscapes,' a category
+that also includes neural networks and optimization problems like constraint
+satisfaction \cite{Mezard_2009_Information}.  The most tractable family of
+these are the mean-field spherical $p$-spin models \cite{Crisanti_1992_The}
+(for a review see \cite{Castellani_2005_Spin-glass}) defined by the energy
 \begin{equation} \label{eq:bare.hamiltonian}
   H_0 = \frac1{p!}\sum_{i_1\cdots i_p}^NJ_{i_1\cdots i_p}z_{i_1}\cdots z_{i_p},
 \end{equation}
 where $J$ is a symmetric tensor whose elements are real Gaussian variables and
-$z\in\mathbb R^N$ is constrained to the sphere $z^2=N$. This problem has been
+$z\in\mathbb R^N$ is constrained to the sphere $z^Tz=N$. This problem has been
 studied in the algebra \cite{Cartwright_2013_The} and probability literature
 \cite{Auffinger_2012_Random, Auffinger_2013_Complexity}.  It has been attacked
 from several angles: the replica trick to compute the Boltzmann--Gibbs
 distribution \cite{Crisanti_1992_The}, a Kac--Rice \cite{Kac_1943_On,
 Rice_1939_The, Fyodorov_2004_Complexity} procedure (similar to the
 Fadeev--Popov integral) to compute the number of saddle-points of the energy
-function \cite{Crisanti_1995_Thouless-Anderson-Palmer}, and the
-gradient-descent---or more generally Langevin---dynamics staring from a
-high-energy configuration \cite{Cugliandolo_1993_Analytical}. Thanks to the
-simplicity of the energy, all these approaches yield analytic results in the
-large-$N$ limit.
+function \cite{Crisanti_1995_Thouless-Anderson-Palmer}, and gradient-descent
+(or more generally Langevin) dynamics starting from a high-energy configuration
+\cite{Cugliandolo_1993_Analytical}. Thanks to the simplicity of the energy, all
+these approaches yield analytic results in the large-$N$ limit.
 
 In this paper we extend the study to complex variables: we shall take
 $z\in\mathbb C^N$ and $J$ to be a symmetric tensor whose elements are
 \emph{complex} normal, with $\overline{|J|^2}=p!/2N^{p-1}$ and
 $\overline{J^2}=\kappa\overline{|J|^2}$ for complex parameter $|\kappa|<1$. The
-constraint remains $z^2=N$.
+constraint remains $z^Tz=N$.
 
-The motivations for this paper are of two types. On the practical side, there
+The motivations for this paper are of three types. On the practical side, there
 are indeed situations in which complex variables appear naturally in disordered
 problems: such is the case in which the variables are \emph{phases}, as in
 random laser problems \cite{Antenucci_2015_Complex}. Quiver Hamiltonians---used
 to model black hole horizons in the zero-temperature limit---also have a
 Hamiltonian very close to ours \cite{Anninos_2016_Disordered}.  A second reason
 is that, as we know from experience, extending a real problem to the complex
-plane often uncovers underlying simplicity that is otherwise hidden, sheding
+plane often uncovers underlying simplicity that is otherwise hidden, shedding
 light on the original real problem, e.g., as in the radius of convergence of a
 series.
 
-Deforming an integral in $N$ real variables to a surface of dimension $N$ in
+Finally, deforming an integral in $N$ real variables to a surface of dimension $N$ in
 $2N$-dimensional complex space has turned out to be necessary for correctly
 defining and analyzing path integrals with complex action (see
 \cite{Witten_2010_A, Witten_2011_Analytic}), and as a useful palliative for the
 sign problem \cite{Cristoforetti_2012_New, Tanizaki_2017_Gradient,
-Scorzato_2016_The}.  In order to do this correctly, the features of landscape
-of the action in complex space---like the relative position of its
-saddles---must be understood. Such landscapes are in general not random: here
-we propose to follow the strategy of computer science of understanding the
-generic features of random instances, expecting that this sheds light on the
+Scorzato_2016_The}.  In order to do this correctly, features of landscape
+of the action in complex space---such as the relative position of saddles and the existence of Stokes lines joining them---must be understood. Such landscapes are in general not random: here
+we follow the standard strategy of computer science by understanding the
+generic features of random instances, expecting that this sheds light on
 practical, nonrandom problems.
 
 Returning to our problem, the spherical constraint is enforced using the method
-of Lagrange multipliers: introducing $\epsilon\in\mathbb C$, our energy is
+of Lagrange multipliers: introducing $\epsilon\in\mathbb C$, our constrained
+energy is
 \begin{equation} \label{eq:constrained.hamiltonian}
   H = H_0+\frac p2\epsilon\left(N-\sum_i^Nz_i^2\right).
 \end{equation}
-We choose to constrain our model by $z^2=N$ rather than $|z|^2=N$ in order to
-preserve the analyticity of $H$. The nonholomorphic constraint also has a
-disturbing lack of critical points nearly everywhere: if $H$ were so
-constrained, then $0=\partial^* H=-p\epsilon z$ would only be satisfied for
-$\epsilon=0$.  
-
-The critical points are of $H$ given by the solutions to the set of equations
+One might balk at the constraint $z^Tz=N$---which could appropriately be called
+a \emph{hyperbolic} constraint---by comparison with $z^\dagger z=N$. The reasoning
+behind the choice is twofold.
+
+First, we seek draw conclusions from our model that are applicable to generic
+holomorphic functions without any symmetry. Samples of $H_0$ nearly provide
+this, save for a single anomaly: the value of the energy and its gradient at
+any point $z$ correlate along the $z$ direction, with $\overline{H_0\partial
+H_0}\propto \overline{H_0(\partial H_0)^*}\propto z$. This anomalous direction
+should thus be forbidden, and the constraint surface $z^Tz=N$ accomplishes this.
+
+Second, taking the constraint to be the level set of a holomorphic function
+means the resulting configuration space is a \emph{bone fide} complex manifold,
+and therefore permits easy generalization of the integration techniques
+referenced above. The same cannot be said for the space defined by $z^\dagger
+z=N$, which is topologically the $(2N-1)$-sphere and cannot admit a complex
+structure.
+
+Imposing the constraint with  a holomorphic function
+makes the resulting configuration space a \emph{bone fide} complex manifold, which is, as we mentioned, the
+situation we wish to model. The same cannot be said for the space defined by $z^\dagger
+z=N$, which is topologically the $(2N-1)$-sphere, does not admit a complex
+structure, and thus yields a trivial structure of saddles.
+However, we will introduce the bound $r^2\equiv z^\dagger z/N\leq R^2$
+on the `radius' per spin as a device to classify saddles.   We shall see that this
+`radius' $r$ and its upper bound $R$ are insightful knobs in our present
+problem, revealing structure as they are varied. Note that taking $R=1$ reduces
+the problem to that of the ordinary $p$-spin.
+
+The critical points are of $H$ given by the solutions to
 \begin{equation} \label{eq:polynomial}
   \frac{p}{p!}\sum_{j_1\cdots j_{p-1}}^NJ_{ij_1\cdots j_{p-1}}z_{j_1}\cdots z_{j_{p-1}}
   = p\epsilon z_i
@@ -114,20 +136,16 @@ equations of degree $p-1$, to which one must add the constraint.  In this sense
 this study also provides a complement to the work on the distribution of zeroes
 of random polynomials \cite{Bogomolny_1992_Distribution}, which are for $N=1$
 and $p\to\infty$.  We see from \eqref{eq:polynomial} that at any critical
-point, $\epsilon=H/N$, the average energy.
-
-Since $H$ is holomorphic, any critical point of $\operatorname{Re}H$ is also a
-critical point of $\operatorname{Im}H$. The number of critical points of $H$ is
-therefore the same as that of $\operatorname{Re}H$. From each saddle 
-emerge gradient lines of $\operatorname{Re}H$, which are also ones of constant
-$\operatorname{Im}H$ and therefore constant phase.
+point $\epsilon=H_0/N$, the average energy.
 
-Writing $z=x+iy$, $\operatorname{Re}H$ can be considered a real-valued function
-of $2N$ real variables. Its number of saddle-points is given by the usual
-Kac--Rice formula:
+Since $H$ is holomorphic, any critical point of $\operatorname{Re}H$ is also
+one of $\operatorname{Im}H$, and therefore of $H$ itself. Writing $z=x+iy$ for
+$x,y\in\mathbb R^N$, $\operatorname{Re}H$ can be considered a real-valued
+function of $2N$ real variables. The number of critical points of $H$ is thus given by the
+usual Kac--Rice formula applied to $\operatorname{Re}H$:
 \begin{equation} \label{eq:real.kac-rice}
   \begin{aligned}
-    \mathcal N_J&(\kappa,\epsilon)
+    \mathcal N&(\kappa,\epsilon,R)
       = \int dx\,dy\,\delta(\partial_x\operatorname{Re}H)\delta(\partial_y\operatorname{Re}H) \\
       &\hspace{6pc}\times\left|\det\begin{bmatrix}
           \partial_x\partial_x\operatorname{Re}H & \partial_x\partial_y\operatorname{Re}H \\
@@ -135,14 +153,21 @@ Kac--Rice formula:
         \end{bmatrix}\right|.
   \end{aligned}
 \end{equation}
-The Cauchy--Riemann equations may be used to write this in a manifestly complex
-way.  With the Wirtinger derivative $\partial=\frac12(\partial_x-i\partial_y)$,
-one can write $\partial_x\operatorname{Re}H=\operatorname{Re}\partial H$ and
+This expression is to be averaged over $J$ to give the complexity $\Sigma$ as
+$N \Sigma= \overline{\log\mathcal N}$, a calculation that involves the replica
+trick. Based on the experience from similar problems \cite{Castellani_2005_Spin-glass},  the
+\emph{annealed approximation} $N \Sigma \sim \log \overline{ \mathcal N}$ is
+expected to be exact wherever the complexity is positive.
+
+The Cauchy--Riemann equations may be used to write \eqref{eq:real.kac-rice} in
+a manifestly complex way.  With the Wirtinger derivative
+$\partial\equiv\frac12(\partial_x-i\partial_y)$, one can write
+$\partial_x\operatorname{Re}H=\operatorname{Re}\partial H$ and
 $\partial_y\operatorname{Re}H=-\operatorname{Im}\partial H$. Carrying these
-transformations through, we have
+transformations through, one finds
 \begin{equation} \label{eq:complex.kac-rice}
   \begin{aligned}
-    \mathcal N_J&(\kappa,\epsilon)
+    \mathcal N&(\kappa,\epsilon,r)
       = \int dx\,dy\,\delta(\operatorname{Re}\partial H)\delta(\operatorname{Im}\partial H) \\
                 &\hspace{6pc}\times\left|\det\begin{bmatrix}
             \operatorname{Re}\partial\partial H & -\operatorname{Im}\partial\partial H \\
@@ -155,8 +180,9 @@ transformations through, we have
   \end{aligned}
 \end{equation}
 This gives three equivalent expressions for the determinant of the Hessian: as
-that of a $2N\times 2N$ real matrix, that of an $N\times N$ Hermitian matrix,
-i.e. the norm squared of that of an $N\times N$ complex symmetric matrix.
+that of a $2N\times 2N$ real symmetric matrix, that of the $N\times N$ Hermitian
+matrix $(\partial\partial H)^\dagger\partial\partial H$, or the norm squared of
+that of the $N\times N$ complex symmetric matrix $\partial\partial H$.
 
 These equivalences belie a deeper connection between the spectra of the
 corresponding matrices. Each positive eigenvalue of the real matrix has a
@@ -167,57 +193,13 @@ Hessian is therefore the same as the distribution of singular values of
 $\partial\partial H$, or the distribution of square-rooted eigenvalues of
 $(\partial\partial H)^\dagger\partial\partial H$.
 
-The expression \eqref{eq:complex.kac-rice} is to be averaged over $J$ to give
-the complexity $\Sigma$ as $N \Sigma= \overline{\log\mathcal N} = \int dJ \,
-\log \mathcal N_J$, a calculation that involves the replica trick. In most the
-parameter-space that we shall study here, the \emph{annealed approximation} $N
-\Sigma \sim \log \overline{ \mathcal N} = \log\int dJ \, \mathcal N_J$ is
-exact. 
-
 A useful property of the Gaussian $J$ is that gradient and Hessian at fixed
-$\epsilon$ are statistically independent \cite{Bray_2007_Statistics,
+energy $\epsilon$ are statistically independent \cite{Bray_2007_Statistics,
 Fyodorov_2004_Complexity}, so that the $\delta$-functions and the Hessian may
-be averaged independently. The $\delta$-functions are converted to exponentials
-by the introduction of auxiliary fields $\hat z=\hat x+i\hat y$.  The average
-of those factors over $J$ can then be performed. A generalized
-Hubbard--Stratonovich allows a change of variables from the $4N$ original
-and auxiliary fields to eight bilinears defined by $Na=|z|^2$, $N\hat a=|\hat
-z|^2$, $N\hat c=\hat z^2$, $Nb=\hat z^*z$, and $Nd=\hat zz$ (and their
-conjugates). The result, to leading order in $N$, is
-\begin{equation} \label{eq:saddle}
-    \overline{\mathcal N}(\kappa,\epsilon)
-        = \int da\,d\hat a\,db\,db^*d\hat c\,d\hat c^*dd\,dd^*e^{Nf(a,\hat a,b,\hat c,d)},
-\end{equation}
-where the argument of the exponential is
-\begin{widetext}
-  \begin{equation}
-    f=2+\frac12\log\det\frac12\begin{bmatrix}
-      1 & a & d & b \\
-      a & 1 & b^* & d^* \\
-      d & b^* & \hat c & \hat a \\
-      b & d^* & \hat a & \hat c^*
-    \end{bmatrix}
-    +\int d\lambda\,\rho(\lambda)\log|\lambda|^2
-    +p\operatorname{Re}\left\{
-      \frac18\left[\hat aa^{p-1}+(p-1)|d|^2a^{p-2}+\kappa(\hat c^*+(p-1)b^2)\right]-\epsilon b
-    \right\}.
-  \end{equation}
-  The integral of the distribution $\rho$ of eigenvalues of $\partial\partial
-  H$ comes from the Hessian and is dependant on $a$ alone. This function has an
-  extremum in $\hat a$, $b$, $\hat c$, and $d$ at which its value is
-  \begin{equation} \label{eq:free.energy.a}
-      f(a)=1+\frac12\log\left(\frac4{p^2}\frac{a^2-1}{a^{2(p-1)}-|\kappa|^2}\right)+\int d\lambda\,\rho(\lambda)\log|\lambda|^2
-      -2C_+[\operatorname{Re}(\epsilon e^{-i\theta})]^2-2C_-[\operatorname{Im}(\epsilon e^{-i\theta})]^2,
-  \end{equation}
-\end{widetext}
-where $\theta=\frac12\arg\kappa$ and
-\begin{equation}
-  C_{\pm}=\frac{a^p(1+p(a^2-1))\mp a^2|\kappa|}{a^{2p}\pm(p-1)a^p(a^2-1)|\kappa|-a^2|\kappa|^2}.
-\end{equation}
-This leaves a single parameter, $a$, which dictates the magnitude of $|z|^2$,
-or alternatively the magnitude $y^2$ of the imaginary part. The latter vanishes
-as $a\to1$, where (as we shall see) one recovers known results for the real
-$p$-spin.
+be averaged independently. First we shall compute the spectrum of the Hessian,
+which can in turn be used to compute the determinant. Then we will treat the
+$\delta$-functions and the resulting saddle point equations. The results of
+these calculations begin around \eqref{eq:bezout}.
 
 The Hessian $\partial\partial H=\partial\partial H_0-p\epsilon I$ is equal to
 the unconstrained Hessian with a constant added to its diagonal. The eigenvalue
@@ -229,28 +211,28 @@ Hessian of the unconstrained Hamiltonian is
   =\frac{p(p-1)}{p!}\sum_{k_1\cdots k_{p-2}}^NJ_{ijk_1\cdots k_{p-2}}z_{k_1}\cdots z_{k_{p-2}},
 \end{equation}
 which makes its ensemble that of Gaussian complex symmetric matrices, when the
-direction along the constraint is neglected. Given its variances
-$\overline{|\partial_i\partial_j H_0|^2}=p(p-1)a^{p-2}/2N$ and
+anomalous direction normal to the constraint surface is neglected. Given its variances
+$\overline{|\partial_i\partial_j H_0|^2}=p(p-1)r^{p-2}/2N$ and
 $\overline{(\partial_i\partial_j H_0)^2}=p(p-1)\kappa/2N$, $\rho_0(\lambda)$ is
 constant inside the ellipse
 \begin{equation} \label{eq:ellipse}
-  \left(\frac{\operatorname{Re}(\lambda e^{i\theta})}{a^{p-2}+|\kappa|}\right)^2+
-  \left(\frac{\operatorname{Im}(\lambda e^{i\theta})}{a^{p-2}-|\kappa|}\right)^2
-  <\frac{p(p-1)}{2a^{p-2}}
+  \left(\frac{\operatorname{Re}(\lambda e^{i\theta})}{r^{p-2}+|\kappa|}\right)^2+
+  \left(\frac{\operatorname{Im}(\lambda e^{i\theta})}{r^{p-2}-|\kappa|}\right)^2
+  <\frac{p(p-1)}{2r^{p-2}}
 \end{equation}
 where $\theta=\frac12\arg\kappa$ \cite{Nguyen_2014_The}. The eigenvalue
 spectrum of $\partial\partial H$ is therefore constant inside the same ellipse
 translated so that its center lies at $-p\epsilon$.  Examples of these
 distributions are shown in the insets of Fig.~\ref{fig:spectra}.
 
-The eigenvalue spectrum of the Hessian of the real part is different from the
-spectrum $\rho(\lambda)$ of $\partial\partial H$, but rather equivalent to the
+The eigenvalue spectrum of the Hessian of the real part is not the
+spectrum $\rho(\lambda)$ of $\partial\partial H$, but instead the
 square-root eigenvalue spectrum of $(\partial\partial H)^\dagger\partial\partial H$;
 in other words, the singular value spectrum $\rho(\sigma)$ of $\partial\partial
 H$. When $\kappa=0$ and the elements of $J$ are standard complex normal, this
 is a complex Wishart distribution. For $\kappa\neq0$ the problem changes, and
 to our knowledge a closed form is not in the literature.  We have worked out an
-implicit form for this spectrum using the replica method.
+implicit form for the singular value spectrum using the replica method.
 
 Introducing replicas to bring the partition function into the numerator of the
 Green function \cite{Livan_2018_Introduction} gives
@@ -263,25 +245,24 @@ Green function \cite{Livan_2018_Introduction} gives
       \right]
     \right\},
   \end{equation}
-  with sums taken over repeated Latin indices.  The average is then made over
+  with sums taken over repeated Latin indices. The average is then made over
   $J$ and Hubbard--Stratonovich is used to change variables to the replica matrices
-  $N\alpha_{\alpha\beta}=(\zeta^{(\alpha)})^*\cdot\zeta^{(\beta)}$ and
-  $N\chi_{\alpha\beta}=\zeta^{(\alpha)}\cdot\zeta^{(\beta)}$ and a series of
-  replica vectors. The replica-symmetric ansatz leaves all off-diagonal
-  elements and vectors zero, and
-  $\alpha_{\alpha\beta}=\alpha_0\delta_{\alpha\beta}$,
+  $N\alpha_{\alpha\beta}=(\zeta^{(\alpha)})^\dagger\zeta^{(\beta)}$ and
+  $N\chi_{\alpha\beta}=(\zeta^{(\alpha)})^T\zeta^{(\beta)}$, and a series of
+  replica vectors. The replica-symmetric ansatz leaves all replica vectors
+  zero, and $\alpha_{\alpha\beta}=\alpha_0\delta_{\alpha\beta}$,
   $\chi_{\alpha\beta}=\chi_0\delta_{\alpha\beta}$. The result is
   \begin{equation}\label{eq:green.saddle}
     \overline G(\sigma)=N\lim_{n\to0}\int d\alpha_0\,d\chi_0\,d\chi_0^*\,\alpha_0
     \exp\left\{nN\left[
-      1+\frac{p(p-1)}{16}a^{p-2}\alpha_0^2-\frac{\alpha_0\sigma}2+\frac12\log(\alpha_0^2-|\chi_0|^2)
+      1+\frac{p(p-1)}{16}r^{p-2}\alpha_0^2-\frac{\alpha_0\sigma}2+\frac12\log(\alpha_0^2-|\chi_0|^2)
       +\frac p4\operatorname{Re}\left(\frac{(p-1)}8\kappa^*\chi_0^2-\epsilon^*\chi_0\right)
     \right]\right\}.
     \nonumber % He's too long, and we don't cite him (now)!
   \end{equation}
 \end{widetext}
 
-\begin{figure}[b]
+\begin{figure}
   \centering
 
   \includegraphics{spectra_00.pdf}
@@ -290,25 +271,27 @@ Green function \cite{Livan_2018_Introduction} gives
   \includegraphics{spectra_15.pdf}
 
   \caption{
-    Eigenvalue and singular value spectra of the matrix $\partial\partial H$
-    for $p=3$, $a=\frac54$, and $\kappa=\frac34e^{-i3\pi/4}$ with (a)
-    $\epsilon=0$, (b) $\epsilon=-\frac12|\epsilon_{\mathrm{th}}|$, (c)
+    Eigenvalue and singular value spectra of the Hessian $\partial\partial H$
+    of the $3$-spin model with $\kappa=\frac34e^{-i3\pi/4}$. Pictured
+    distributions are for critical points at `radius' $r=\sqrt{5/4}$ and with
+    energy per spin (a) $\epsilon=0$, (b)
+    $\epsilon=-\frac12|\epsilon_{\mathrm{th}}|$, (c)
     $\epsilon=-|\epsilon_{\mathrm{th}}|$, and (d)
-    $\epsilon=-\frac32|\epsilon_{\mathrm{th}}|$. The shaded region of each inset
-    shows the support of the eigenvalue distribution. The solid line on each
-    plot shows the distribution of singular values, while the overlaid
-    histogram shows the empirical distribution from $2^{10}\times2^{10}$ complex
-    normal matrices with the same covariance and diagonal shift as
-    $\partial\partial H$.
+    $\epsilon=-\frac32|\epsilon_{\mathrm{th}}|$. The shaded region of each
+    inset shows the support of the eigenvalue distribution \eqref{eq:ellipse}.
+    The solid line on each plot shows the distribution of singular values
+    \eqref{eq:spectral.density}, while the overlaid histogram shows the
+    empirical distribution from $2^{10}\times2^{10}$ complex normal matrices
+    with the same covariance and diagonal shift as $\partial\partial H$.
   } \label{fig:spectra}
 \end{figure}
 
 The argument of the exponential has several saddles. The solutions $\alpha_0$
 are the roots of a sixth-order polynomial, and the root with the smallest value
-of $\operatorname{Re}\alpha_0$ in all the cases we studied gives the correct
-solution. A detailed analysis of the saddle point integration is needed to
-understand why this is so. Given such $\alpha_0$, the density of singular
-values follows from the jump across the cut, or
+of $\operatorname{Re}\alpha_0$ gives the correct solution in all the cases we
+studied. A detailed analysis of the saddle point integration is needed to
+understand why this is so. Evaluated at such a solution, the density of
+singular values follows from the jump across the cut, or
 \begin{equation} \label{eq:spectral.density}
   \rho(\sigma)=\frac1{i\pi N}\left(
     \lim_{\operatorname{Im}\sigma\to0^+}\overline G(\sigma)
@@ -318,77 +301,109 @@ values follows from the jump across the cut, or
 Examples can be seen in Fig.~\ref{fig:spectra} compared with numeric
 experiments.
 
-The transition from a one-cut to two-cut singular value spectrum naturally
-corresponds to the origin leaving the support of the eigenvalue spectrum.
-Weyl's theorem requires that the product over the norm of all eigenvalues must
-not be greater than the product over all singular values \cite{Weyl_1912_Das}.
-Therefore, the absence of zero eigenvalues implies the absence of zero singular
-values. The determination of the threshold energy -- the energy at which the
-distribution of singular values becomes gapped -- is then reduced to a
-geometry problem, and yields
+The formation of a gap in the singular value spectrum naturally corresponds to
+the origin leaving the support of the eigenvalue spectrum.  Weyl's theorem
+requires that the product over the norm of all eigenvalues must not be greater
+than the product over all singular values \cite{Weyl_1912_Das}.  Therefore, the
+absence of zero eigenvalues implies the absence of zero singular values. The
+determination of the threshold energy---the energy at which the distribution
+of singular values becomes gapped---is reduced to the geometry problem of
+determining when the boundary of the ellipse defined in \eqref{eq:ellipse}
+intersects the origin, and yields
 \begin{equation} \label{eq:threshold.energy}
   |\epsilon_{\mathrm{th}}|^2
-  =\frac{p-1}{2p}\frac{(1-|\delta|^2)^2a^{p-2}}
+  =\frac{p-1}{2p}\frac{(1-|\delta|^2)^2r^{p-2}}
   {1+|\delta|^2-2|\delta|\cos(\arg\kappa+2\arg\epsilon)}
 \end{equation}
-for $\delta=\kappa a^{-(p-2)}$.
+for $\delta=\kappa r^{-(p-2)}$. Notice that the threshold depends on both the
+energy per spin $\epsilon$ on the `radius' $r$ of the saddle.
+
+We will now address the $\delta$-functions of \eqref{eq:complex.kac-rice}.
+These are converted to exponentials by the introduction of auxiliary fields
+$\hat z=\hat x+i\hat y$.  The average over $J$ can then be performed. A
+generalized Hubbard--Stratonovich allows a change of variables from the $4N$
+original and auxiliary fields to eight bilinears defined by $Nr=z^\dagger z$,
+$N\hat r=\hat z^\dagger\hat z$, $Na=\hat z^\dagger z$, $Nb=\hat z^Tz$, and
+$N\hat c=\hat z^T\hat z$ (and their conjugates). The result, to leading order
+in $N$, is
+\begin{equation} \label{eq:saddle}
+    \overline{\mathcal N}(\kappa,\epsilon,R)
+        = \int dr\,d\hat r\,da\,da^*db\,db^*d\hat c\,d\hat c^*e^{Nf(r,\hat r,a,b,\hat c)},
+\end{equation}
+where the argument of the exponential is
+\begin{widetext}
+  \begin{equation}
+    f=2+\frac12\log\det\frac12\begin{bmatrix}
+      1 & r & b & a \\
+      r & 1 & a^* & b^* \\
+      b & a^* & \hat c & \hat r \\
+      a & b^* & \hat r & \hat c^*
+    \end{bmatrix}
+    +\int d\lambda\,d\lambda^*\rho(\lambda)\log|\lambda|^2
+    +p\operatorname{Re}\left\{
+      \frac18\left[\hat rr^{p-1}+(p-1)|b|^2r^{p-2}+\kappa(\hat c^*+(p-1)a^2)\right]-\epsilon a
+    \right\}.
+  \end{equation}
+  The spectrum $\rho$ is given in \eqref{eq:ellipse} and is dependant on $r$ alone. This function has an
+  extremum in $\hat r$, $a$, $b$, and $\hat c$ at which its value is
+  \begin{equation} \label{eq:free.energy.a}
+      f=1+\frac12\log\left(\frac4{p^2}\frac{r^2-1}{r^{2(p-1)}-|\kappa|^2}\right)+\int d\lambda\,\rho(\lambda)\log|\lambda|^2
+      -2C_+[\operatorname{Re}(\epsilon e^{-i\theta})]^2-2C_-[\operatorname{Im}(\epsilon e^{-i\theta})]^2,
+  \end{equation}
+\end{widetext}
+where $\theta=\frac12\arg\kappa$ and
+\begin{equation}
+  C_{\pm}=\frac{r^p(1+p(r^2-1))\mp r^2|\kappa|}{r^{2p}\pm(p-1)r^p(r^2-1)|\kappa|-r^2|\kappa|^2}.
+\end{equation}
+Notice that level sets of $f$ in energy $\epsilon$ also give ellipses, but of
+different form from the ellipse in \eqref{eq:ellipse}.
 
-Given $\rho$, the integral in \eqref{eq:free.energy.a} may be preformed for
-arbitrary $a$. The resulting expression is maximized for $a=\infty$ for all
-values of $\kappa$ and $\epsilon$. Taking this saddle gives
+This expression is maximized for $r=R$, its value at the boundary, for
+all values of $\kappa$ and $\epsilon$. Evaluating the complexity at this
+saddle, in the limit of unbounded spins, gives
 \begin{equation} \label{eq:bezout}
-  \log\overline{\mathcal N}(\kappa,\epsilon)
+  \lim_{R\to\infty}\log\overline{\mathcal N}(\kappa,\epsilon,R)
   =N\log(p-1).
 \end{equation}
-This is, to this order,  precisely the Bézout bound, the maximum number of
-solutions that $N$ equations of degree $p-1$ may have
-\cite{Bezout_1779_Theorie}. That we saturate this bound is perhaps not
-surprising, since the coefficients of our polynomial equations
-\eqref{eq:polynomial} are complex and have no symmetries. Reaching Bézout in
-\eqref{eq:bezout} is not our main result, but it provides a good check.
-Analogous asymptotic scaling has been found for the number of pure Higgs states
-in supersymmetric quiver theories \cite{Manschot_2012_From}.
-
-More insight is gained by looking at the count as a function of $a$, defined by
-$\overline{\mathcal N}(\kappa,\epsilon,a)=e^{Nf(a)}$.  In the large-$N$ limit,
-this is the cumulative number of critical points, or the number of critical
-points $z$ for which $|z|^2\leq a$. We likewise define the $a$-dependant
-complexity $\Sigma(\kappa,\epsilon,a)=N\log\overline{\mathcal
-N}(\kappa,\epsilon,a)$
+This is, to leading order, precisely the Bézout bound, the maximum number of
+solutions to $N$ equations of degree $p-1$ \cite{Bezout_1779_Theorie}. That we
+saturate this bound is perhaps not surprising, since the coefficients of our
+polynomial equations \eqref{eq:polynomial} are complex and have no symmetries.
+Reaching Bézout in \eqref{eq:bezout} is not our main result, but it provides a
+good check.  Analogous asymptotic scaling has been found for the number of pure
+Higgs states in supersymmetric quiver theories \cite{Manschot_2012_From}.
 
 \begin{figure}[htpb]
   \centering
   \includegraphics{complexity.pdf}
   \caption{
-    The complexity of the 3-spin model at $\epsilon=0$ as a function of
-    $a=|z|^2=1+y^2$ at several values of $\kappa$. The dashed line shows
+    The complexity of the 3-spin model as a function of the maximum `radius'
+    $R$ at zero energy and several values of $\kappa$. The dashed line shows
     $\frac12\log(p-1)$, while the dotted shows $\log(p-1)$.
   } \label{fig:complexity}
 \end{figure}
 
-Everything is analytically tractable for $\epsilon=0$, giving
+For finite $R$, everything is analytically tractable at $\epsilon=0$:
 \begin{equation} \label{eq:complexity.zero.energy}
-  \Sigma(\kappa,0,a)
-  =\log(p-1)-\frac12\log\left(\frac{1-|\kappa|^2a^{-2(p-1)}}{1-a^{-2}}\right).
+  \Sigma(\kappa,0,R)
+  =\log(p-1)-\frac12\log\left(\frac{1-|\kappa|^2R^{-4(p-1)}}{1-R^{-4}}\right).
 \end{equation}
-Notice that the limit of this expression as $a\to\infty$ corresponds with
-\eqref{eq:bezout}, as expected. This is plotted as a function of $a$ for
-several values of $\kappa$ in Fig.~\ref{fig:complexity}. For any $|\kappa|<1$,
-the complexity goes to negative infinity as $a\to1$, i.e., as the spins are
-restricted to the reals.  This is natural, given that the $y$ contribution to
-the volume shrinks to zero as that of an $N$-dimensional sphere $\sum_i y_i^2=N(a-1)$ with volume
-$\sim(a-1)^N$.  However, when the result is analytically continued to
+This is plotted as a function of $R$ for several values of $\kappa$ in
+Fig.~\ref{fig:complexity}. For any $|\kappa|<1$, the complexity goes to
+negative infinity as $R\to1$, i.e., as the spins are restricted to the reals.
+This is natural, since volume of configuration space vanishes in this limit
+like $(R^2-1)^N$. However, when the result is analytically continued to
 $\kappa=1$ (which corresponds to real $J$) something novel occurs: the
-complexity has a finite value at $a=1$. Since the $a$-dependence gives a
-cumulative count, this implies a $\delta$-function density of critical points
-along the line $y=0$.  The number of critical points contained within is
+complexity has a finite value at $R=1$. This implies a $\delta$-function
+density of critical points on the $r=1$ (or $y=0$) boundary.  The number of
+critical points contained there is
 \begin{equation}
-  \lim_{a\to1}\lim_{\kappa\to1}\log\overline{\mathcal N}(\kappa,0,a)
+  \lim_{R\to1}\lim_{\kappa\to1}\log\overline{\mathcal N}(\kappa,0,R)
   = \frac12N\log(p-1),
 \end{equation}
 half of \eqref{eq:bezout} and corresponding precisely to the number of critical
-points of the real $p$-spin model (note the role of conjugation symmetry,
-already underlined in \cite{Bogomolny_1992_Distribution}). The full
+points of the real $p$-spin model. (Note the role of conjugation symmetry,
+already underlined in \cite{Bogomolny_1992_Distribution}.) The full
 $\epsilon$-dependence of the real $p$-spin is recovered by this limit as
 $\epsilon$ is varied.
 
@@ -396,23 +411,35 @@ $\epsilon$ is varied.
   \centering
   \includegraphics{desert.pdf}
   \caption{
-    The minimum value of $a$ for which the complexity is positive as a function
-    of (real) energy $\epsilon$  for the 3-spin model at several values of
-    $\kappa$.
+    The value of bounding `radius' $R$ for which $\Sigma(\kappa,\epsilon,R)=0$ as a
+    function of (real) energy per spin $\epsilon$  for the 3-spin model at
+    several values of $\kappa$. Above each line the complexity is positive and
+    critical points proliferate, while below it the complexity is negative and
+    critical points are exponentially suppressed. The dotted black lines show
+    the location of the ground and highest exited state energies for the real
+    3-spin model.
   } \label{fig:desert}
 \end{figure}
 
+In the thermodynamic limit, \eqref{eq:complexity.zero.energy} implies that most
+critical points are concentrated at infinite radius $r$. For finite $N$ the
+average radius of critical points is likewise finite. By differentiating
+$\overline{\mathcal N}$ with respect to $R$ and normalizing, one obtains the
+distribution of critical points as a function of $r$. This yields an average
+radius proportional to $N^{1/4}$.  One therefore expects typical critical
+points to have a norm that grows modestly with system size.
+
 These qualitative features carry over to nonzero $\epsilon$. In
-Fig.~\ref{fig:desert} we show that for $\kappa<1$ there is always a gap of $a$
-close to one for which there are no solutions. When $\kappa=1$---the analytic
-continuation to the real computation---the situation is more interesting. In
-the range of energies where there are real solutions this gap closes, which is
-only possible if the density of solutions diverges at $a=1$.  Another
-remarkable feature of this limit is that there is still a gap without solutions
-around `deep' real energies where there is no real solution. A moment's thought
-tells us that this is a necessity: otherwise a small perturbation of the $J$s
-could produce an unusually deep solution to the real problem, in a region where
-this should not happen.
+Fig.~\ref{fig:desert} we show that for $\kappa<1$ there is always a gap in $r$
+close to one in which solutions are exponentially suppressed. When
+$\kappa=1$---the analytic continuation to the real computation---the situation
+is more interesting. In the range of energies where there are real solutions
+this gap closes, which is only possible if the density of solutions diverges at
+$r=1$. Outside this range, around `deep' real energies where real solutions are
+exponentially suppressed, the gap remains. A moment's thought tells us that
+this is necessary: otherwise a small perturbation of the $J$s could produce
+an unusually deep solution to the real problem, in a region where this should
+not happen.
 
 \begin{figure}[t]
   \centering
@@ -425,37 +452,38 @@ this should not happen.
   \caption{
     Energies at which states exist (green shaded region) and threshold energies
     (black solid line) for the 3-spin model with
-    $\kappa=\frac34e^{-i3\pi/4}$ and (a) $a=2$, (b) $a=1.325$, (c) $a=1.125$,
-    and (d) $a=1$. No shaded region is shown in (d) because no states exist at
+    $\kappa=\frac34e^{-i3\pi/4}$ and (a) $r=\sqrt2$, (b) $r=\sqrt{1.325}$, (c) $r=\sqrt{1.125}$,
+    and (d) $r=1$. No shaded region is shown in (d) because no states exist at
     any energy.
   } \label{fig:eggs}
 \end{figure}
 
 The relationship between the threshold and ground, or extremal, state energies
-is richer than in the real case. In Fig.~\ref{fig:eggs} these are shown in the
+is richer than in the real case.  In Fig.~\ref{fig:eggs} these are shown in the
 complex-$\epsilon$ plane for several examples. Depending on the parameters, the
-threshold might always come at smaller magnitude than the extremal state, or
-always come at larger magnitude, or cross as a function of complex argument.
-For sufficiently large $a$ the threshold always comes at larger magnitude than
-the extremal state. If this were to happen in the real case, it would likely
-imply our replica symmetric computation is unstable, since having a ground
-state above the threshold implies a ground state Hessian with many negative
-eigenvalues, a contradiction. However, this is not an obvious contradiction in
-the complex case. The relationship between the threshold, i.e., where the gap
-appears, and the dynamics of, e.g., a minimization algorithm or physical
-dynamics, are a problem we hope to address in future work.
-
- This paper provides a first step towards the study of a complex landscape with
- complex variables. The next obvious one is to study the topology of the
+threshold might have a smaller or larger magnitude than the extremal state, or
+cross as a function of complex argument.  For sufficiently large $r$ the
+threshold is always at a larger magnitude. If this were to happen in the real
+case, it would likely imply our replica symmetric computation were unstable,
+since having a ground state above the threshold implies a ground state Hessian
+with many negative eigenvalues, a contradiction. However, this is not an
+contradiction in the complex case, where the energy is not bounded from below.
+The relationship between the threshold, i.e., where the gap appears, and the
+dynamics of, e.g., a minimization algorithm, deformed integration cycle, or
+physical dynamics, are a problem we hope to address in future work.
+
+ This paper provides a first step towards the study of complex landscapes with
+ complex variables. The next obvious step is to study the topology of the
  critical points, the sets reached following  gradient descent (the
  Lefschetz thimbles), and ascent (the anti-thimbles) \cite{Witten_2010_A,
  Witten_2011_Analytic, Cristoforetti_2012_New, Behtash_2017_Toward,
  Scorzato_2016_The}, which act as constant-phase integrating `contours.'
  Locating and counting the saddles that are joined by gradient lines---the
  Stokes points, which play an important role in the theory---is also well within
- reach of the present-day spin-glass literature techniques.    We anticipate
- that the threshold level, where the system develops a mid-spectrum gap, will
- play a crucial role as it does in the real case.
+ reach of the present-day spin-glass literature techniques. We anticipate
+ that the threshold level, where the system develops a mid-spectrum gap, plays
+ a crucial role in determining whether these Stokes points proliferate under
+ some continuous change of parameters.
 
 \begin{acknowledgments}
   We wish to thank Alexander Altland, Satya Majumdar and Gregory Schehr for a useful suggestions.
diff --git a/complexity.pdf b/complexity.pdf
index b68f2cf..cd4864d 100644
--- a/complexity.pdf
+++ b/complexity.pdf
diff --git a/desert.pdf b/desert.pdf
index e19484a..c9e03df 100644
--- a/desert.pdf
+++ b/desert.pdf
diff --git a/referee_respose.txt b/referee_respose.txt
new file mode 100644
index 0000000..3b204c6
--- /dev/null
+++ b/referee_respose.txt
@@ -0,0 +1,158 @@
+----------------------------------------------------------------------
+Response to Referee A -- LZ16835/Kent-Dobias
+----------------------------------------------------------------------
+
+Referee A wrote:
+> The authors consider the mean-field p-spin spherical model with
+> *complex* variables and study the number of saddle points of the
+> energy and the eigenvalue distribution of their Hessian matrix. The
+> main result of the rather technical computation is that in a
+> particular limit (concretely kappa->1) the known results for the real
+> p-spin spherical model are reproduced, the (expected) Bézout bound for
+> the number of solutions of the saddle point equations is reached and
+> that the relationships between the “threshold” and extremal state
+> energies is richer in the complex case than in the real case.
+> 
+> I must admit that I was not able to grasp any far-reaching
+> consequences of the computational tour de force only hinted at in the
+> manuscript, and I fear that a nonexpert reader would also not be able
+> to do so. Two arguments are pushed forward by the authors to justify
+> the dissemination of their results to the broader readership of PRL:
+> One is that there are indeed situations in which complex variables
+> appear naturally in disordered system. The first example the authors
+> mention is a Hamiltonian that could be relevant for with random Laser
+> problems and was analyzed 2015 in PRA, which has up to now 30
+> citations according to Google Scholar, and the second example is a
+> Hamiltonian from sting theory that was analyzed in 2016 in JHEP, which
+> has up to now 31 citations. I do not feel that these two examples
+> prove that the enumeration of saddle points if the p-spin model is
+> important or of broad interest.
+> 
+> The second argument of the authors is that extending a real problem to
+> the complex plane often uncovers underlying simplicity that is
+> otherwise hidden, shedding light on the original real problem. Here I
+> come back to what I already mentioned above: I do not see any
+> simplicity emerging from the present calculation and I also do not see
+> the original problem in a new light. Therefore, I do not think that
+> one of the four PRL criteria is actually fulfilled and I recommend to
+> transfer the manuscript to PRE.
+
+We disagree with the referee's assessment here, as we have also explained in
+our letter to the editors. Something in particular that goes unaddressed is
+another motivation (which in the referee's defense we did not enumerate clearly
+in our draft): that understanding the distribution of complex critical points
+is necessary in the treatment of a large class of integrals involved both in
+the definition of quantum mechanics with a complex action and in ameliorating
+the sign problem in, e.g., lattice QCD.
+
+If the criteria for publication is to be "first past the post" of cited
+citations, one might examine our citations of that literature:
+
+  - Analytic continuation of Chern-Simons theory, E Witten (2011): 444 citations
+
+  - New approach to the sign problem in quantum field theories: High density
+    QCD on a Lefschetz thimble, M Cristoforetti et al (2012): 285 citations
+
+Both works are concerned with the location and relative positions of critical
+points of complex theories. In the resubmitted manuscript we have better
+emphasized this motivation.
+
+Referee A wrote:
+> Although the first part of the manuscript is well written and well
+> understandable (at least for me) from page 2 on it becomes very
+> technical and unreadable for a non-expert. If the reader skips to the
+> results and tries to understand the figures she/he is left with the
+> ubiquitous parameter a, whose physical meaning is hidden deep in the
+> saddle point calculation (“dictates the magnitude of |z|^2” – well,
+> with respect to the solutions of (3): is “a” the average value of the
+> modulus squared of the solution z’s or not?). Similar with epsilon:
+> apparently it is the average energy of the saddle point solution – why
+> not writing so also in the figure captions? The paper would profit a
+> lot from a careful rewriting of at least the result section and to
+> provide figure captions with the physical meaning of the quantities
+> and parameters shown.
+
+We thank the referee for their helpful suggestions with regards to the
+readability of our manuscript. In the resubmitted version, much has been
+rewritten for clarity. We would like to highlight several of the most
+substantive changes:
+
+  - The ubiquitous parameter 'a' was replaced by the more descriptive 'r^2', as
+    it is a sort of radius, along with a new parameter 'R^2' which bounds it.
+    Descriptions in English of these were added to the figure captions.
+
+  - The technical portion of the paper was reordered to connect better with the
+    sections preceding and following it.
+
+  - The location of the results is now indicated before the beginning of the
+    technical portion for readers interested in skipping ahead.
+
+Referee A wrote:
+> A couple of minor, technical, quibbles:
+> 
+> 1) If there is any real world application of a p-spin model with
+> complex variables it will NOT have a spherical constraint. I would
+> suggest to discuss the consequences of this constraint, which is
+> introduced for computational simplicity.
+> 
+> 2) After eq. (2): ”We choose to constrain our model by z^2=N.“ Then it
+> is not a spherical constraint any more – does it have any physical
+> relevance?
+
+We have added a more detailed discussion of the constraint to address these
+confusions, emphasizing its purpose. The new paragraphs are:
+
+> One might balk at the constraint $z^Tz=N$---which could appropriately be
+> called a \emph{hyperbolic} constraint---by comparison with $z^\dagger z=N$.
+> The reasoning behind the choice is twofold.
+>
+> First, we seek draw conclusions from our model that are applicable to generic
+> holomorphic functions without any symmetry. Samples of $H_0$ nearly provide
+> this, save for a single anomaly: the value of the energy and its gradient at
+> any point $z$ correlate along the $z$ direction, with $\overline{H_0\partial
+> H_0}\propto \overline{H_0(\partial H_0)^*}\propto z$. This anomalous
+> direction should thus be forbidden, and the constraint surface $z^Tz=N$
+> accomplishes this.
+>
+> Second, taking the constraint to be the level set of a holomorphic function
+> means the resulting configuration space is a \emph{bone fide} complex
+> manifold, and therefore permits easy generalization of the integration
+> techniques referenced above. The same cannot be said for the space defined by
+> $z^\dagger z=N$, which is topologically the $(2N-1)$-sphere and cannot admit
+> a complex structure.
+>
+> Imposing the constraint with  a holomorphic function makes the resulting
+> configuration space a \emph{bone fide} complex manifold, which is, as we
+> mentioned, the situation we wish to model. The same cannot be said for the
+> space defined by $z^\dagger z=N$, which is topologically the $(2N-1)$-sphere,
+> does not admit a complex structure, and thus yields a trivial structure of
+> saddles.  However, we will introduce the bound $r^2\equiv z^\dagger z/N\leq
+> R^2$ on the `radius' per spin as a device to classify saddles.   We shall see
+> that this `radius' $r$ and its upper bound $R$ are insightful knobs in our
+> present problem, revealing structure as they are varied. Note that taking
+> $R=1$ reduces the problem to that of the ordinary $p$-spin.
+
+Referee A wrote:
+> 3) On p.2: “…a, which dictates the magnitude of |z|^2, or
+> alternatively the magnitude y^2 of the imaginary part. The last part
+> is hard to understand, should be explained.
+
+We thank the referee for pointing out this confusing statement, which was
+unnecessary and removed.
+
+> 4) On p.2: “In most the parameter space we shall study her, the
+> annealed approximation is exact.” I think it is necessary to provide
+> some evidence her, because the annealed approximation is usually a
+> pretty severe approximation.
+
+We have nuanced the statement in question and added a citation to a review
+article which outlines the reasoning for analogous models. The amended sentence
+reads:
+
+> Based on the experience from similar problems \cite{Castellani_2005_Spin-glass},
+> the \emph{annealed approximation} $N\Sigma\sim\log\overline{\mathcal N}$ is
+> expected to be exact wherever the complexity is positive.
+
+Sincerely,
+Jaron Kent-Dobias & Jorge Kurchan
+
diff --git a/threshold.pdf b/threshold.pdf
new file mode 100644
index 0000000..538d3d8
--- /dev/null
+++ b/threshold.pdf