Lots of filling out the details of the complexity calculation.

1 files changed, 201 insertions, 103 deletions

diff --git a/2-point.tex b/2-point.tex
index bb1b625..380d677 100644
--- a/2-point.tex
+++ b/2-point.tex
@@ -37,6 +37,8 @@
 
 \cite{Ros_2020_Distribution, Ros_2019_Complex, Ros_2019_Complexity}
 
+\section{Model}
+
 The mixed $p$-spin models are defined by the Hamiltonian
 \begin{equation} \label{eq:hamiltonian}
   H(\mathbf s)=-\sum_p\frac1{p!}\sum_{i_1\cdots i_p}^NJ^{(p)}_{i_1\cdots i_p}s_{i_1}\cdots s_{i_p}
@@ -48,7 +50,7 @@ the energy is typically extensive. The overbar denotes an average
 over the coefficients $J$. The factors $a_p$ in the variances are freely chosen
 constants that define the particular model. For instance, the `pure'
 $p$-spin model has $a_{p'}=\delta_{p'p}$. This class of models encompasses all
-statistically isotropic gaussian random Hamiltonians defined on the
+statistically isotropic Gaussian random Hamiltonians defined on the
 hypersphere.
 
 The covariance between the energy at two different points is a function of the overlap, or dot product, between those points, or
@@ -59,52 +61,92 @@ where the function $f$ is defined from the coefficients $a_p$ by
 \begin{equation}
   f(q)=\frac12\sum_pa_pq^p
 \end{equation}
-In this paper, we will focus on models with a replica symmetric complexity.
+In this paper, we will focus on models with a replica symmetric complexity, but
+many of the intermediate formulae are valid for arbitrary replica symmetry
+breakings.
+
+To enforce the spherical constraint at stationary points, we make use of a Lagrange multiplier $\omega$. This results in the extremal problem
+\begin{equation}
+  H(\mathbf s)+\frac\omega2(\|\mathbf s\|^2-N)
+\end{equation}
+The gradient and Hessian at a stationary point are then
+\begin{align}
+  \nabla H(\mathbf s,\omega)=\partial H(\mathbf s)+\omega\mathbf s
+  &&
+  \operatorname{Hess}H(\mathbf s,\omega)=\partial\partial H(\mathbf s)+\omega I
+\end{align}
+where $\partial=\frac\partial{\partial\mathbf s}$ will always denote the derivative with respect to $\mathbf s$.
 
 We introduce the Kac--Rice \cite{Kac_1943_On, Rice_1944_Mathematical} measure
 \begin{equation}
-  d\nu_H(\mathbf s)
-  =d\mathbf s\,\delta\big(\nabla H(\mathbf s)\big)\,
-  \big|\det\operatorname{Hess}H(\mathbf s)\big|
+  d\nu_H(\mathbf s,\omega)
+  =2\,d\mathbf s\,d\omega\,\delta(\|\mathbf s\|^2-N)\,
+  \delta\big(\nabla H(\mathbf s,\omega)\big)\,
+  \big|\det\operatorname{Hess}H(\mathbf s,\omega)\big|
 \end{equation}
 which counts stationary points of the function $H$. More interesting is the
 measure conditioned on the energy density $E$ and stability $\mu$,
 \begin{equation}
-  d\nu_H(\mathbf s\mid E,\mu)
-  =d\nu_H(\mathbf s)\,
+  d\nu_H(\mathbf s,\omega\mid E,\mu)
+  =d\nu_H(\mathbf s,\omega)\,
   \delta\big(H(\mathbf s)-NE\big)\,
-  \delta\big(N\mu-\operatorname{Tr}\operatorname{Hess}H(\mathbf s)\big)
+  \delta\big(N\mu-\operatorname{Tr}\operatorname{Hess}H(\mathbf s,\omega)\big)
 \end{equation}
 
+\section{Complexity}
 
 We want the typical number of stationary points with energy density
-$E_2$ and stability $\mu_2$ that lie a fixed overlap $q$ from a reference
-stationary point of energy density $E_1$ and stability $\mu_1$.
-\begin{align*}
+$E_1$ and stability $\mu_1$ that lie a fixed overlap $q$ from a reference
+stationary point of energy density $E_0$ and stability $\mu_0$.
+\begin{equation} \label{eq:complexity.definition}
   \Sigma_{12}
-  &=\frac1N\overline{\int\frac{d\nu_H(\mathbf s_0\mid E_0,\mu_0)}{\int d\nu_H(\mathbf s_0'\mid E_0,\mu_0)}\,
-  \log\bigg(\int d\nu_H(\mathbf s_1\mid E_1,\mu_1)\,\delta(Nq-\mathbf s_0\cdot\mathbf s_1)\bigg)}
-\end{align*}
-\begin{align*}
+  =\frac1N\overline{\int\frac{d\nu_H(\pmb\sigma,\varsigma\mid E_0,\mu_0)}{\int d\nu_H(\pmb\sigma',\varsigma'\mid E_0,\mu_0)}\,
+  \log\bigg(\int d\nu_H(\mathbf s,\omega\mid E_1,\mu_1)\,\delta(Nq-\pmb\sigma\cdot\mathbf s)\bigg)}
+\end{equation}
+Both the denominator and the logarithm are treated using the replica trick, which yields
+\begin{equation}
   \Sigma_{12}
-  &=\frac1N\lim_{n\to0}\lim_{m\to-1}\overline{\int d\nu_H(\mathbf s_0\mid E_0,\mu_0)\left(\int d\nu_H(\mathbf s_0'\mid E_0,\mu_0)\right)^m\,
-  \frac\partial{\partial n}\bigg(\int d\nu_H(\mathbf s_1\mid E_1,\mu_1)\,\delta(Nq-\mathbf s_0\cdot \mathbf s_1)\bigg)^n}\\
-  &=\frac1N\lim_{n\to0}\lim_{m\to0}\frac\partial{\partial n}\overline{\int\left(\prod_{b=1}^md\nu_H(\pmb\sigma_b\mid E_0,\mu_0)\right)\left(\prod_{a=1}^nd\nu_H(\mathbf s_a\mid E_1,\mu_1)\,\delta(Nq-\pmb \sigma_1\cdot \mathbf s_a)\right)}
-\end{align*}
+  =\frac1N\lim_{n\to0}\lim_{m\to0}\frac\partial{\partial n}\overline{\int\left(\prod_{b=1}^md\nu_H(\pmb\sigma_b,\varsigma_b\mid E_0,\mu_0)\right)\left(\prod_{a=1}^nd\nu_H(\mathbf s_a,\omega_a\mid E_1,\mu_1)\,\delta(Nq-\pmb \sigma_1\cdot \mathbf s_a)\right)}
+\end{equation}
+Note that because of the structure of \eqref{eq:complexity.definition},
+$\pmb\sigma_1$ is special among the set of $\pmb\sigma$ replicas, since only it
+is constrained to lie a given overlap from the $\mathbf s$ replicas. This
+replica asymmetry will be important later.
+
+\subsection{The Hessian factors}
+
+The double partial derivatives of the energy are Gaussian with the variance
 \begin{equation}
-  \overline{\big|\det\operatorname{Hess}H(s)\big|\,\delta\big(N\mu-\operatorname{Tr}\operatorname{Hess}H(s)\big)}
-  =e^{N\int d\lambda\,\rho(\lambda+\mu)\log|\lambda|}\delta(N\mu-s\cdot\partial H)
+  \overline{(\partial_i\partial_jH(\mathbf s))^2}=\frac1Nf''(1)
 \end{equation}
+which means that the matrix of partial derivatives belongs to the GOE class. Its spectrum is given by the Wigner semicircle
 \begin{equation}
   \rho(\lambda)=\begin{cases}
     \frac2{\pi}\sqrt{1-\big(\frac{\lambda}{\mu_\text m}\big)^2} & \lambda^2\leq\mu_\text m^2 \\
     0 & \text{otherwise}
   \end{cases}
 \end{equation}
+with radius $\mu_\text m=\sqrt{4f''(1)}$. Since the Hessian differs from the
+matrix of partial derivatives by adding the constant diagonal matrix $\omega
+I$, it follows that the spectrum of the Hessian is a Winger semicircle shifted
+by $\omega$, or $\rho(\lambda+\omega)$.
+
+The average over factors depending on the Hessian alone can be made separately
+from those depending on the gradient or energy, since for random Gaussian
+fields the Hessian is independent of these \cite{Bray_2007_Statistics}. In
+principle the fact that we have conditioned the Hessian to belong to stationary
+points of certain energy, stability, and proximity to another stationary point
+will modify its statistics, but these changes will only appear at subleading
+order in $N$ \cite{Ros_2019_Complexity}. At leading order, the various expectations factorize, each yielding
+\begin{equation}
+  \overline{\big|\det\operatorname{Hess}H(\mathbf s,\omega)\big|\,\delta\big(N\mu-\operatorname{Tr}\operatorname{Hess}H(\mathbf s,\omega)\big)}
+  =e^{N\int d\lambda\,\rho(\lambda+\mu)\log|\lambda|}\delta(N\mu-N\omega)
+\end{equation}
+Therefore, all of the Lagrange multipliers are fixed identically to the stabilities $\mu$. We define the function
 \begin{equation}
   \begin{aligned}
     \mathcal D(\mu)
-    &=\int d\lambda\,\rho(\lambda+\mu)\ln|\lambda| \\
+    &=\int d\lambda\,\rho(\lambda+\mu)\log|\lambda| \\
     &=\begin{cases}
       \frac12+\log\left(\frac12\mu_\text m\right)+\frac\mu{\mu_\text m}\left(\frac\mu{\mu_\text m}-\sqrt{\big(\frac\mu{\mu_\text m}\big)^2-1}\right)
       -\log\left(\frac{\mu}{\mu_\text m}-\sqrt{\big(\frac\mu{\mu_\text m}\big)^2-1}\right) & \mu>\mu_\text m \\
@@ -115,28 +157,92 @@ stationary point of energy density $E_1$ and stability $\mu_1$.
     \end{cases}
   \end{aligned}
 \end{equation}
+and the full factor due to the Hessians is
+\begin{equation}
+  e^{Nm\mathcal D(\mu_0)+Nn\mathcal D(\mu_1)}\left[\prod_a^m\delta(N\mu_0-N\varsigma_a)\right]\left[\prod_a^n\delta(N\mu_1-N\omega_a)\right]
+\end{equation}
 
-\begin{align}
-  \mathcal Q_{00}=\begin{bmatrix}
-    \hat\beta_0\\\hat\mu_0\\C^{00}\\R^{00}\\D^{00}
-  \end{bmatrix}
-  &&
-  \mathcal Q_{11}=\begin{bmatrix}
-    \hat\beta_1\\\hat\mu_1\\C^{11}\\R^{11}\\D^{11}
-  \end{bmatrix}
-  &&
-  \mathcal Q_{01}=\begin{bmatrix}
-    \hat\mu_{01}\\C^{01}\\R^{01}\\R_{10}\\D^{01}
-  \end{bmatrix}
+\subsection{The other factors}
+
+Having integrated over the Lagrange multipliers using the $\delta$ functions
+resulting from the average of the Hessians, the remaining part of the integrand
+has the form
+\begin{equation}
+  e^{
+    -Nm\hat\beta_0E_0-Nn\hat\beta_1E_1
+    -\sum_a^m\left[(\pmb\sigma_a\cdot\hat{\pmb\sigma}_a)\mu_0
+      -\frac12\hat\mu_0(N-\pmb\sigma_a\cdot\pmb\sigma_a)
+    \right]
+      -\sum_a^n\left[(\mathbf s_a\cdot\hat{\mathbf s}_a)\mu_1
+      -\frac12\hat\mu_1(N-\mathbf s_a\cdot\mathbf s_a)
+      -\frac12\hat\mu_{12}(Nq-\pmb\sigma_1\cdot\mathbf s_a)
+    \right]
+    +\int d\mathbf t\,\mathcal O(\mathbf t)H(\mathbf t)
+  }
+\end{equation}
+where we have introduced the linear operator
+\begin{equation}
+  \mathcal O(\mathbf t)
+  =\sum_a^m\delta(\mathbf t-\pmb\sigma_a)\left(
+    i\hat{\pmb\sigma}_a\cdot\partial_{\mathbf t}+\hat\beta_0
+    \right)
+    +
+    \sum_a^n\delta(\mathbf t-\mathbf s_a)\left(
+    i\hat{\mathbf s}_a\cdot\partial_{\mathbf t}+\hat\beta_1
+    \right)
+\end{equation}
+We have written the $H$-dependant terms in this strange form for the ease of taking the average over $H$: since it is Gaussian-correlated, it follows that
+\begin{equation}
+  \overline{e^{\int d\mathbf t\,\mathcal O(\mathbf t)H(\mathbf t)}}
+  =e^{\frac12\int d\mathbf t\,d\mathbf t'\,\mathcal O(\mathbf t)\mathcal O(\mathbf t')\overline{H(\mathbf t)H(\mathbf t')}}
+  =e^{N\frac12\int d\mathbf t\,d\mathbf t'\,\mathcal O(\mathbf t)\mathcal O(\mathbf t')f\big(\frac{\mathbf t\cdot\mathbf t'}N\big)}
+\end{equation}
+It remains only to apply the doubled operators to $f$ and then evaluate the simple integrals over the $\delta$ measures. We do not include these details, which are standard.
+
+\subsection{Hubbard--Stratonovich}
+
+Having expanded this expression, we are left with an argument in the exponential which is a function of scalar products between the fields $\mathbf s$, $\hat{\mathbf s}$, $\pmb\sigma$, and $\hat{\pmb\sigma}$. We will change integration coordinates from these fields to matrix fields given by the scalar products, defined as
+\begin{align} \label{eq:fields}
+  C^{00}_{ab}=\frac1N\pmb\sigma_a\cdot\pmb\sigma_b &&
+  R^{00}_{ab}=-i\frac1N{\pmb\sigma}_a\cdot\hat{\pmb\sigma}_b &&
+  D^{00}_{ab}=\frac1N\hat{\pmb\sigma}_a\cdot\hat{\pmb\sigma}_b \\
+  C^{01}_{ab}=\frac1N\pmb\sigma_a\cdot\mathbf s_b &&
+  R^{01}_{ab}=-i\frac1N{\pmb\sigma}_a\cdot\hat{\mathbf s}_b &&
+  R^{10}_{ab}=-i\frac1N\hat{\pmb\sigma}_a\cdot{\mathbf s}_b &&
+  D^{01}_{ab}=\frac1N\hat{\pmb\sigma}_a\cdot\hat{\mathbf s}_b \\
+  C^{11}_{ab}=\frac1N\mathbf s_a\cdot\mathbf s_b &&
+  R^{11}_{ab}=-i\frac1N{\mathbf s}_a\cdot\hat{\mathbf s}_b &&
+  D^{11}_{ab}=\frac1N\hat{\mathbf s}_a\cdot\hat{\mathbf s}_b
 \end{align}
+We insert into the integral the product of $\delta$ functions enforcing these
+definitions, integrated over the new matrix fields, which is equivalent to
+multiplying by one. Once this is done, the many scalar products appearing
+throughout can be replaced by the matrix fields, and the original vector fields
+can be integrated over. Conjugate matrix field integrals created when the
+$\delta$ functions are promoted to exponentials can be evaluated by saddle
+point in the standard way, yielding an effective action depending on the above
+matrix fields alone.
+
+We will always assume that the square matrices $C^{00}$, $R^{00}$, $D^{00}$,
+$C^{11}$, $R^{11}$, and $D^{11}$ are hierarchical matrices, with each set of
+three sharing the same hierarchical structure. In particular, we immediately
+define $c_\mathrm d^{00}$, $r_\mathrm d^{00}$, $d_\mathrm d^{00}$, $c_\mathrm d^{11}$, $r_\mathrm d^{11}$, and
+$d_\mathrm d^{11}$ as the value of the diagonal elements of these matrices,
+respectively. Note that $c_\mathrm d^{00}=c_\mathrm d^{11}=1$ due to the spherical constraint.
+
+Defining the `block' fields $\mathcal Q_{00}=(\hat\beta_0, \hat\mu_0, C^{00},
+R^{00}, D^{00})$, $\mathcal Q_{11}=(\hat\beta_1, \hat\mu_1, C^{11}, R^{11},
+D^{11})$, and $\mathcal Q_{01}=(\hat\mu_{01},C^{01},R^{01},R^{10},D^{01})$
+the resulting complexity is
 \begin{equation}
   \Sigma_{01}
   =\frac1N\lim_{n\to0}\lim_{m\to0}\frac\partial{\partial n}\int d\mathcal Q_{00}\,d\mathcal Q_{11}\,d\mathcal Q_{01}\,e^{Nm\mathcal S_0(\mathcal Q_{00})+Nn\mathcal S_1(\mathcal Q_{00},\mathcal Q_{11},\mathcal Q_{01})}
 \end{equation}
+where
 \begin{equation}
   \begin{aligned}
     &\mathcal S_0(\mathcal Q_{00})
-    =-\hat\beta_0E_0-r^{00}_d\mu_0-\frac12\hat\mu_0(1-c^{00}_d)+\mathcal D(\mu_0)\\
+    =-\hat\beta_0E_0-r^{00}_\mathrm d\mu_0-\frac12\hat\mu_0(1-c^{00}_\mathrm d)+\mathcal D(\mu_0)\\
     &\quad+\frac1m\bigg\{
       \frac12\sum_{ab}^m\left[
         \hat\beta_1^2f(C^{00}_{ab})-(2\hat\beta_1R^{00}_{ab}+D^{00}_{ab})f'(C^{00}_{ab})+(R_{ab}^{00})^2f''(C_{ab}^{00})
@@ -144,11 +250,11 @@ stationary point of energy density $E_1$ and stability $\mu_1$.
     \bigg\}
   \end{aligned}
 \end{equation}
-
+is the action for the ordinary, one-point complexity, and remainder is given by
 \begin{equation}
   \begin{aligned}
     &\mathcal S(\mathcal Q_{00},\mathcal Q_{11},\mathcal Q_{01})
-    =-\hat\beta_1E_1-\mu_1r^{11}_d-\frac12\hat\mu_1(1-c^{11}_d) \\
+    =-\hat\beta_1E_1-r^{11}_\mathrm d\mu_1-\frac12\hat\mu_1(1-c^{11}_\mathrm d)+\mathcal D(\mu_1) \\
     &\quad+\frac1n\sum_b^n\left\{-\frac12\hat\mu_{12}(q-C^{01}_{1b})+\sum_a^m\left[
       \hat\beta_0\hat\beta_1f(C^{01}_{ab})-(\hat\beta_0R^{01}_{ab}+\hat\beta_1R^{10}_{ab}+D^{01}_{ab})f'(C^{01}_{ab})+R^{01}_{ab}R^{10}_{ab}f''(C^{01}_{ab})
   \right]\right\}
@@ -174,10 +280,15 @@ stationary point of energy density $E_1$ and stability $\mu_1$.
     \bigg\}
   \end{aligned}
 \end{equation}
+Because of the structure of this problem in the twin limits of $m$ and $n$ to
+zero, the parameters $\mathcal Q_{00}$ can be evaluated at a saddle point of
+$\mathcal S_0$ alone. This means that these parameters will take the same value
+they take when the ordinary, 1-point complexity is calculated.
 
+The $m\times n$ matrices $C^{01}$, $R^{01}$, $R^{10}$, and $D^{01}$ are
+expected to have the following form at the saddle point:
 \begin{align}
-  C^{01}
-  =
+  C^{01}=
   \begin{subarray}{l}
     \hphantom{[}\begin{array}{ccc}\leftarrow&n&\rightarrow\end{array}\hphantom{\Bigg]}\\
   \left[
@@ -189,8 +300,7 @@ stationary point of energy density $E_1$ and stability $\mu_1$.
     \end{array}
   \right]\begin{array}{c}
     \\\uparrow\\m-1\\\downarrow
-  \end{array}\\
-  \vphantom{\begin{array}{c}n\end{array}}
+  \end{array}
   \end{subarray}
   &&
   R^{01}
@@ -217,8 +327,10 @@ stationary point of energy density $E_1$ and stability $\mu_1$.
     0&\cdots&0
   \end{bmatrix}
 \end{align}
+where only the first row is nonzero as a result of the sole linear term
+proportional to $C_{1b}^{01}$ in the action.
 
-The inverse of block hierarchical matrix is still a block hierarchical matrix, since (dropping the superscripts for clarity)
+The inverse of block hierarchical matrix is still a block hierarchical matrix, since
 \begin{equation}
   \begin{bmatrix}
     C^{00}&iR^{00}\\iR^{00}&D^{00}
@@ -231,72 +343,35 @@ The inverse of block hierarchical matrix is still a block hierarchical matrix, s
 \end{equation}
 Because of the structure of the 01 matrices, the volume element will depend only on the diagonal if this matrix. If we write
 \begin{align}
-  \tilde c_d^{00}&=[(C^{00}D^{00}+R^{00}R^{00})^{-1}C^{00}]_{11} \\
-  \tilde r_d^{00}&=[(C^{00}D^{00}+R^{00}R^{00})^{-1}R^{00}]_{11} \\
-  \tilde d_d^{00}&=[(C^{00}D^{00}+R^{00}R^{00})^{-1}D^{00}]_{11}
+  \tilde c_\mathrm d^{00}&=[(C^{00}D^{00}+R^{00}R^{00})^{-1}C^{00}]_{\text d} \\
+  \tilde r_\mathrm d^{00}&=[(C^{00}D^{00}+R^{00}R^{00})^{-1}R^{00}]_{\text d} \\
+  \tilde d_\mathrm d^{00}&=[(C^{00}D^{00}+R^{00}R^{00})^{-1}D^{00}]_{\text d}
 \end{align}
-
-In the replica symmetric case,
-\begin{align}
-  \tilde c_d^{00}=\frac1{(r^{00}_d)^2+d^{00}_d} &&
-  \tilde r_d^{00}=\frac{r^{00}_d}{(r^{00}_d)^2+d^{00}_d} &&
-  \tilde d_d^{00}=\frac{d^{00}_d}{(r^{00}_d)^2+d^{00}_d}
-\end{align}
-
+then the result is
 \begin{equation}
-  \begin{bmatrix}
-    q^2\tilde d_d^{00}+2qr_{10}\tilde r^{00}_d-r_{10}^2\tilde d^{00}_d
-    &
-    i\left[d_{01}(r_{10}\tilde c^{00}_d-q\tilde r^{00}_d)+r_{01}(r_{10}\tilde r^{00}_d+q\tilde d^{00}_d)\right]
-    \\
-    i\left[d_{01}(r_{10}\tilde c^{00}_d-q\tilde r^{00}_d)+r_{01}(r_{10}\tilde r^{00}_d+q\tilde d^{00}_d)\right]
-    &
-    d_{01}^2\tilde c^{00}_d+2r_{01}d_{01}\tilde r^{00}_d-r_{01}^2\tilde d^{00}_d
-  \end{bmatrix}
+  \begin{aligned}
+     &   \begin{bmatrix}
+          C^{01}&iR^{01}\\iR^{10}&D^{01}
+        \end{bmatrix}^T
+        \begin{bmatrix}
+          C^{00}&iR^{00}\\iR^{00}&D^{00}
+        \end{bmatrix}^{-1}
+        \begin{bmatrix}
+          C^{01}&iR^{01}\\iR^{10}&D^{01}
+        \end{bmatrix} \\
+    &\qquad=\begin{bmatrix}
+      q^2\tilde d_\mathrm d^{00}+2qr_{10}\tilde r^{00}_\mathrm d-r_{10}^2\tilde d^{00}_\mathrm d
+      &
+      i\left[d_{01}(r_{10}\tilde c^{00}_\mathrm d-q\tilde r^{00}_\mathrm d)+r_{01}(r_{10}\tilde r^{00}_\mathrm d+q\tilde d^{00}_\mathrm d)\right]
+      \\
+      i\left[d_{01}(r_{10}\tilde c^{00}_\mathrm d-q\tilde r^{00}_\mathrm d)+r_{01}(r_{10}\tilde r^{00}_\mathrm d+q\tilde d^{00}_\mathrm d)\right]
+      &
+      d_{01}^2\tilde c^{00}_\mathrm d+2r_{01}d_{01}\tilde r^{00}_\mathrm d-r_{01}^2\tilde d^{00}_\mathrm d
+    \end{bmatrix}
+  \end{aligned}
 \end{equation}
 where each block is a constant $n\times n$ matrix.
 
-In the twin limits of $m$ and $n$ to zero, the saddle point conditions for the variables involving only the reference critical point (those in $\mathcal Q_{00}$) reduce to the ordinary, 1-point conditions. With a replica-symmetric ansatz, these conditions are
-\begin{align}
-  \hat\beta_0
-  &=-\frac{(\epsilon_0+\mu_0)f'(1)+\epsilon_0f''(1)}{f(1)\big(f'(1)+f''(1)\big)-f'(1)^2}\\
-  r_d^{00}
-  &=\frac{\mu_0f(1)+\epsilon_0f'(1)}{f(1)\big(f'(1)+f''(1)\big)-f'(1)^2} \\
-  d_d^{00}
-  &=\frac1{f'(1)}
-  -\left(
-    \frac{\mu_0f(1)+\epsilon_0f'(1)}{f(1)\big(f'(1)+f''(1)\big)-f'(1)^2}
-    \right)^2
-\end{align}
-
-\begin{align*}
-  &
-  \begin{bmatrix}
-    \tilde c&\tilde r\\\tilde r&\tilde d
-  \end{bmatrix}
-  =
-  \begin{bmatrix}
-    q&ir_{10}\\ir_{01}&d_{01}
-  \end{bmatrix}
-  \begin{bmatrix}
-    1&ir_{0}\\
-    ir_{0}&d_{0}
-  \end{bmatrix}^{-1}
-  \begin{bmatrix}
-    q&ir_{01}\\ir_{10}&d_{01}
-  \end{bmatrix}\\
-  &=
-  \frac1{r_{0}^2+d_{0}}\begin{bmatrix}
-    q^2d_{0}+2qr_{0}r_{10}-r_{10}^2
-    &
-    i\left[d_{01}(r_{10}-r_0q)+r_{01}(r_0r_{10}+d_0q)\right]
-    \\
-    i\left[d_{01}(r_{10}-r_0q)+r_{01}(r_0r_{10}+d_0q)\right]
-    &
-    d_{01}^2+2r_{0}r_{01}d_{01}-d_{0}r_{01}^2
-  \end{bmatrix}
-\end{align*}
-This matrix with modify the diagonal of the RS matrix for the second spin.
 
 Define $\tilde C=C-\tilde c$, $\tilde R=R-\tilde r$, $\tilde D=D-\tilde d$. Then
 \begin{align*}
@@ -351,6 +426,29 @@ Solving, we get
     \right)
 \end{align*}
 
+\section{Replica symmetric case}
+
+With a replica-symmetric ansatz, these conditions are
+\begin{align}
+  \hat\beta_0
+  &=\frac{(\epsilon_0+\mu_0)f'(1)+\epsilon_0f''(1)}{f(1)\big(f'(1)+f''(1)\big)-f'(1)^2}\\
+  r_\mathrm d^{00}
+  &=\frac{\mu_0f(1)+\epsilon_0f'(1)}{f(1)\big(f'(1)+f''(1)\big)-f'(1)^2} \\
+  d_\mathrm d^{00}
+  &=\frac1{f'(1)}
+  -\left(
+    \frac{\mu_0f(1)+\epsilon_0f'(1)}{f(1)\big(f'(1)+f''(1)\big)-f'(1)^2}
+    \right)^2
+\end{align}
+
+In the replica symmetric case,
+\begin{align}
+  \tilde c_\mathrm d^{00}=\frac1{(r^{00}_\mathrm d)^2+d^{00}_\mathrm d} &&
+  \tilde r_\mathrm d^{00}=\frac{r^{00}_\mathrm d}{(r^{00}_\mathrm d)^2+d^{00}_\mathrm d} &&
+  \tilde d_\mathrm d^{00}=\frac{d^{00}_\mathrm d}{(r^{00}_\mathrm d)^2+d^{00}_\mathrm d}
+\end{align}
+
+
 \begin{equation}
   \hat\beta_2E_2-r_{22}^{(0)}\mu_2\frac12\left\{
     \hat\beta_2^2\big(f(1)-f(q_{22}^{(0)})\big)