Finite-Size Relaxational Dynamics of a Spike Random Matrix Spherical Model

Abstract

We present a thorough numerical analysis of the relaxational dynamics of the Sherrington–Kirkpatrick spherical model with an additive non-disordered perturbation for large but finite sizes N. In the thermodynamic limit and at low temperatures, the perturbation is responsible for a phase transition from a spin glass to a ferromagnetic phase. We show that finite-size effects induce the appearance of a distinctive slow regime in the relaxation dynamics, the extension of which depends on the size of the system and also on the strength of the non-disordered perturbation. The long time dynamics are characterized by the two largest eigenvalues of a spike random matrix which defines the model, and particularly by the statistics concerning the gap between them. We characterize the finite-size statistics of the two largest eigenvalues of the spike random matrices in the different regimes, sub-critical, critical, and super-critical, confirming some known results and anticipating others, even in the less studied critical regime. We also numerically characterize the finite-size statistics of the gap, which we hope may encourage analytical work which is lacking. Finally, we compute the finite-size scaling of the long time relaxation of the energy, showing the existence of power laws with exponents that depend on the strength of the non-disordered perturbation in a way that is governed by the finite-size statistics of the gap.

1. Introduction

Quenched random interactions are known to be at the origin of complex behavior in many body systems, both in thermodynamics and dynamics [1,2,3]. The understanding of the properties of this kind of systems in the thermodynamic limit has been steadily growing over the last 40 years or so, mainly through research into solutions of mean field, fully connected models, and also from numerical simulations of finite dimensional ones. On the other side, the behavior of systems composed of a large but finite number of degrees of freedom is much less understood. This case is relevant in many applications in many branches of science, e.g., optimization and inference algorithms [4], biological populations [5], and neural networks [6], to cite but a few. Powerful techniques, such as the saddle point method, are not so useful for studying systems far from the thermodynamic limit. There are a few classes of models in which both the thermodynamics and the dynamical behavior can be solved exactly and still show interesting non-trivial properties qualitatively similar to those of more complex systems. Well-known examples include systems in which the degrees of freedom obey a spherical constraint. One of these models with quenched random pairwise interactions, the Spherical Sherrington–Kirkpatrick model (SSK), allows an exact solution of its thermodynamic properties using tools from Random Matrix Theory [7], without the need to use the more involved replica formalism, necessary for the Ising case. The Langevin relaxational dynamics of the SSK model were solved in [8], where it was shown that the long time relaxation is slow, e.g., with the energy density decaying with a power law in time after a quench from a high temperature initial state to a temperature below the spin glass transition temperature. Interesting out-of-equilibrium features of the dynamics, such as the phenomenon of aging, are present in the model and were completely characterized. The solution of both the thermodynamics and the dynamics of the model was possible due to the knowledge of the spectral properties of the random interactions matrix. In the case of the Gaussian Orthogonal Ensemble, the relevant information is in the Wigner semi-circle density of eigenvalues. When the size of the matrix, $N\times N$, is large but finite, the situation changes. The support of the eigenvalue density is not limited anymore, and the probability distribution of the eigenvalues on the soft edge is given by the celebrated Tracy–Widom $\beta$ distributions [9,10], where $\beta =1,2,4$ refer to the orthogonal, unitary, and symplectic ensembles, respectively. Fluctuations of the free energy of the SSK model were studied, e.g., in [11,12,13,14,15]. These fluctuations are governed by the statistics of the largest eigenvalue of the GOE interaction matrix. The finite N fluctuations in the Langevin dynamics were studied in [16,17]. In these works, a new algebraic (power law) scaling regime was found and characterized, not present in the $N\to \infty$ regime. In the dynamical context, besides the relevance of the largest eigenvalue, which is directly proportional to the ground state energy, the gap between the two largest eigenvalues is also a fundamental quantity to compute the relevant time/size scalings of the long time relaxation.

An interesting related model is the SSK model supplemented with an additive Curie–Weiss term in the Hamiltonian, or equivalently, where the original random interaction matrix is perturbed by a rank one matrix which has the effect of shifting the average value of the random matrix elements from zero to a non zero value. The thermodynamics of this model were also solved in the original work by Kosterlitz et al. [7]. At low enough temperatures, the model presents a phase transition from a spin glass to a ferromagnetic phase at a critical value of the relative strength between the random interactions and the Curie–Weiss one. In the mathematics literature, these kinds of random matrices with finite rank perturbations are called “spike random matrices”. There is a large body of work devoted to the study of the spectral properties of spike random matrices [18,19,20,21,22,23,24,25]. Of special interest to the physics community is the result, originally presented in [19], of a sharp phase transition in the statistics of the largest eigenvalue of particular classes of spike random matrices. Of course, this phase transition has an immediate interpretation in the context of the thermodynamic and also dynamic transitions in the SSK model and related ones. Recently, a renewed interest in the statistical behavior of spike random matrices is manifested in several works focusing on different applications, e.g., overlaps between eigenvectors of correlated spike random matrices [26], analysis of optimal learning rates in non-convex optimization [27], low-rank matrix estimation [28], limits of detection of planted states [29], and ruggedness of complex energy landscapes [30]. From a dynamical perspective, an understanding of spike random matrix models may shed light on problems such as the feasibility of identifying a deterministic signal in a random environment, the reconstruction of hidden patterns in a complex landscape, or the efficiency of search algorithms.

In this work, we perform a numerical study of the statistics of the largest eigenvalues and the gap between the two largest ones in spike random matrices from the GOE ensemble. With the information gained, we then describe the relaxation of the excess energy from the ground state of the spike SSK model following a quench from a high-temperature initial state directly to zero temperature. We show that, as is the case in the standard SSK model, the relaxation shows a new scaling regime, present when the system transitions from the spin glass to the ferromagnetic phases, i.e., a critical scaling regime. In this critical sector, we show that the scaling relaxation behavior is governed by a one-parameter scaling function which depends on the relative strength of the random and Curie–Weiss term and also on the size of the system.

The paper is organized as follows: in Section 2 we introduce the model studied and summarize some known results which will be useful later; in Section 3 we present a numerical study of the statistics of the two largest eigenvalues and the gap of a rank one spike GOE matrix; in Section 4, we present our results on the relaxation dynamics of the model, both in the thermodynamic limit and for large but finite system sizes. Finally, in Section 5, we present a brief discussion of the work and our conclusions.

2. The Model

The spherical Sherrington–Kirkpatrick (SSK) model with a Curie–Weiss (CW) perturbation is described by the following Hamiltonian:

$$\begin{array}{cc}\hfill \mathcal{H}[\overrightarrow{S},z]& ={\mathcal{H}}_{SSK}\left(\overrightarrow{S}\right)+{\mathcal{H}}_{CW}\left(\overrightarrow{S}\right)+\frac{z}{2}\left(\overrightarrow{S}{}^2-N\right)\hfill \\ \hfill & =-\frac{1}{2}\overrightarrow{S}{}^T\mathit{M}\overrightarrow{S}+\frac{z}{2}\left(\overrightarrow{S}{}^2-N\right)\hfill \\ \hfill & =-\frac{1}{2\sqrt{N}}\sum _{i\ne j}^N{J}_{ij}{s}_i{s}_j-\frac{\theta }{2N}\sum _{i\ne j}^N{s}_i{s}_j+\frac{z}{2}\left(\sum _i^N{s}_i^2-N\right),\hfill \end{array}$$

(1)

where z is a Lagrange multiplier which enforces the spherical constraint:

$$\overrightarrow{S}{}^2=\sum _i^N{s}_i^2=N.$$

(2)

In the previous expressions, the spin variables $s_i\in [-\sqrt{N},\sqrt{N}]$ are described by a N component vector $\overrightarrow{S}=(s_1,\dots ,s_N)$. The coupling constants $J_{ij}$ are chosen form a real and symmetric random matrix from the Gaussian Orthogonal Ensemble (GOE), $\mathit{J}={\left\{J_{ij}\right\}}_{(i,j)\in {[1,N]}^2}$, with zero mean and variance $J^2$. $\theta \in {\mathbb{R}}^{+}$ measures the intensity of the deterministic perturbation. Thus, $\mathit{M}$ is a real symmetric $N\times N$ spike random matrix, whose off-diagonal elements are Gaussian distributed with the following mean and variance:

$$p\left(M_{ij}\right)=\mathcal{N}\left[\mu =\frac{\theta }{N},{\sigma }^2=\frac{J^2}{N}\right].$$

(3)

The diagonal elements are zero. The Hamiltonian can be rewritten decomposing $\overrightarrow{S}$ as a linear combination of the eigenvectors $\left\{{\overrightarrow{V}}_{\mu }\right\}$ of the coupling matrix $\mathit{M}$, with ${\overrightarrow{V}}_{\mu }·\overrightarrow{V_{\nu }}={\delta }_{\mu \nu }$. Therefore, with the following notation $s_{\mu }=\overrightarrow{S}·{\overrightarrow{V}}_{\mu }$ for the projections of $\overrightarrow{S}$ on the eigenvectors of $\mathit{M}$, the Hamiltonian becomes:

$$\mathcal{H}[\overrightarrow{S},z]=-\frac{1}{2}\sum _{\mu =1}^N({\lambda }_{\mu }-z)s_{\mu }^2-\frac{z}{2}N,$$

(4)

with ${\left\{{\lambda }_{\mu }\right\}}_{\mu \in [1,N]}$ being the set of N eigenvalues of $\mathit{M}$, with associated eigenvectors $\left\{{\overrightarrow{V}}_{\mu }\right\}$. The eigenvalues are organized such that $\mathrm{Max}\left[{\lambda }_{\mu }\right]={\lambda }_1>{\lambda }_2>\dots >{\lambda }_N=\mathrm{Min}\left[{\lambda }_{\mu }\right]$. In the large N limit, the eigenvalue density distribution of $\mathit{M}$ is given by [31]:

$$\rho \left(\lambda \right)=\left\{\begin{array}{cc}{\rho }_W\left(\lambda \right)\hfill & ,\theta \le J\hfill \\ {\rho }_W\left(\lambda \right)+\frac{1}{N}\delta \left[\lambda -\left(\theta +J^2/\theta \right)\right]\hfill & ,\theta >J\hfill \end{array}\right.,$$

(5)

where ${\rho }_W\left(\lambda \right)$ is the Wigner semicircle law:

$${\rho }_W\left(\lambda \right)=\left\{\begin{array}{cc}\frac{{\left(4J^2-{\lambda }^2\right)}^{1/2}}{2\pi J^2}\hfill & ,\left|\lambda \right|<2J\hfill \\ 0\hfill & ,\left|\lambda \right|>2J\hfill \end{array}\right..$$

(6)

The result (5) means that, if $\theta \le J$, the spectrum of $\mathit{M}$ is given by the Wigner law, corresponding to the GOE ensemble. Otherwise, if $\theta >J$, the largest eigenvalue ${\lambda }_1$ detaches from the Wigner semicircle, becoming an outlier with a delta peak at $\lambda =\theta +J^2{\theta }^{-1}$.

The following set of Langevin equations governs the overdamped dynamics of the model:

$$\begin{array}{cc}\hfill \frac{\partial s_i\left(t\right)}{\partial t}& =-{\delta }_{s_i}\mathcal{H}[\overrightarrow{S},z]+{\xi }_i\left(t\right)\hfill \\ \hfill & =\sum _{j\ne i}^N{M}_{ij}{s}_j\left(t\right)-z(t,\left\{s_{\mu }\left(0\right)\right\})s_i\left(t\right)+{\xi }_i\left(t\right),\forall i\in [1,N],\hfill \end{array}$$

(7)

where ${\xi }_i\left(t\right)$ represents a Gaussian white noise with zero mean and variance $\langle {\xi }_i\left(t\right){\xi }_j\left(t^{\prime }\right)\rangle =2T{\delta }_{ij}\delta (t-t^{\prime })$ and T is the temperature of a thermal bath. As in [17], here we are interested in the zero temperature limit of the Langevin equations, which, in the eigenbasis of $\mathit{M}$, reads:

$$\frac{\partial s_{\mu }\left(t\right)}{\partial t}=\left[{\lambda }_{\mu }-z(t,\left\{s_{\mu }\left(0\right)\right\})\right]s_{\mu }\left(t\right),\forall \mu \in [1,N].$$

(8)

At long times, the system must fall in a stable or metastable state of the free energy which, at $T=0$, reduces to the Hamiltonian $\mathcal{H}[\overrightarrow{S},z]$. The asymptotic stationary state will depend on the initial conditions $\left\{s_{\mu }\left(0\right)\right\}$, as stated explicitly in (8). By setting ${lim}_{t\to \infty }{\partial }_t{s}_{\mu }\left(t\right)=0$ we obtain the criteria:

$${\delta }_{s_{\mu }}\mathcal{H}[\overrightarrow{S},z]=-({\lambda }_{\mu }-z)s_{\mu }=0,\forall \mu \in [1,N],$$

(9)

complemented by the spherical constraint in Equation (2). This system of equations admits the $2N$ solutions:

$$\overrightarrow{S}=\pm \sqrt{N}{\overrightarrow{V}}_{\mu }\mathrm{and}z={\lambda }_{\mu }\forall \mu \in [1,N].$$

(10)

Their stability is determined by the Hessian ${\delta }_{s_{\mu }}{\delta }_{s_{\nu }}\mathcal{H}[\overrightarrow{S},z]=-({\lambda }_{\nu }-z){\delta }_{\mu \nu }$. Taking a given metastable state $\overrightarrow{S}=\sqrt{N}{\overrightarrow{V}}_{\mu }$, the local landscape has $N-\mu$ stable directions, $\mu -1$ unstable directions and a marginal flat one. The energy of each of these configurations is equal to $-{\lambda }_{\mu }N/2$. Thus, the system should always equilibrate in one of the solutions $\pm \sqrt{N}{\overrightarrow{V}}_1$, as they are the only stable ones. The ground state energy density is then simply given by $e_{\mathrm{eq}}=-{\lambda }_1/2$.

Our primary interest here is to describe the behavior, at long times, of the excess energy density, $\Delta e(t,N)=e(t,N)-e_{\mathrm{eq}}\left(N\right)$, for arbitrary system sizes N. To this end, recalling results in [16,17], it can be shown that the time-dependent Lagrange multiplier has a simple relation with the energy density, $z(t,N)=-2e(t,N)$, leading to the exact expression for the excess energy density:

$$\Delta e(t,N)=\frac{{\lambda }_1}{2}+e(t,N)=\frac{1}{2}\frac{{\sum }_{\mu =2}^N{s}_{\mu }^2\left(0\right)({\lambda }_1-{\lambda }_{\mu })e^{2({\lambda }_{\mu }-{\lambda }_1)t}}{s_1^2\left(0\right)+{\sum }_{\mu =2}^N{s}_{\mu }^2\left(0\right)e^{2({\lambda }_{\mu }-{\lambda }_1)t}}.$$

(11)

Because of the dependence of the above expression on the relaxation rates ${\lambda }_{\mu }-{\lambda }_1$, one expects that the late dynamics of the model will be dominated by the gap, $g={\lambda }_1-{\lambda }_2$, between the two largest eigenvalues of the random matrix $\mathit{M}$:

$$\Delta e(t,N)\stackrel{t\to \infty }{\to }\frac{1}{2}\frac{s_2^2\left(0\right)}{s_1^2\left(0\right)}g{e}^{-2gt}.$$

(12)

In the previous expression, there are two sources of fluctuations: the statistics of the gap and the initial conditions. The vector of initial conditions $\overrightarrow{S}\left(0\right)=\left\{s_{\mu }\left(0\right)\right\}$ can be written in the basis of eigenvectors of the $\mathit{M}$ matrix in the form $\overrightarrow{S}\left(0\right)={\sum }_{\nu }{c}_{\nu }{\overrightarrow{V}}_{\nu }$, where the coefficients $c_{\nu }$ satisfy the condition $S^2\left(0\right)={\sum }_{\nu \eta }{c}_{\nu }{c}_{\eta }{\overrightarrow{V}}_{\nu }·{\overrightarrow{V}}_{\eta }={\sum }_{\nu }{c}_{\nu }^2=N$. In the present work, we are primarily interested in a flat distribution on the basis of eigenvectors, that is $c_{\nu }=1$ for all $\nu$, which can be associated to thermal equilibrium at a very high temperature. It corresponds to:

$$s_{\mu }\left(0\right)=s_{\mu }^{\mathrm{flat}}\left(0\right)=1\forall \mu \in [1,N].$$

(13)

The relaxation dynamics of the model depend on the form of the eigenvalue density (5). The effect of the delta contribution is to induce a phase transition when the intensity of the Curie–Weiss term attains the value $\theta =J$ in the thermodynamic limit. Although the CW term remains weaker than the random couplings’ intensity, $\theta <J$, the system behaves like the pure SSK model, relaxing towards a disordered ground state $\overrightarrow{S}$ with a characteristic slow dynamics, as described in [8,16,17]. At finite temperatures, the thermodynamics corresponds to a spin glass phase, originally described in [7]. On the other hand, if the perturbation is strong enough, $\theta >J$, the largest eigenvalue detaches from the bulk of the spectrum, inducing a fast relaxation towards a ferromagnetic ground state, where all the spin variables $s_i$ align in the same direction. For not too high temperatures, exactly at $\theta =J$, the system goes through a continuous phase transition between a disordered spin glass phase and a ferromagnetically ordered one, in the thermodynamic limit [7]. In [13], finite-size fluctuations of the free energy of the model at both sides of the spin glass-ferromagnetic transition were characterized. Here, we are interested in characterizing the finite-size fluctuations of the relaxation dynamics, following a quench from an infinite temperature initial state down to zero temperature for different values of the CW perturbation intensity.

Considering random initial conditions as given by (13), at long times, the behavior of the average excess energy is given by:

$$\mathbb{E}[\Delta e(t,N)]\stackrel{t\to \infty }{\to }\frac{1}{2}\mathbb{E}\left[g{e}^{-2gt}\right].$$

(14)

At present, the statistical properties of the gap $g={\lambda }_1-{\lambda }_2$ are not known. In order to describe its approximate behavior, in the following section, we will pursue a thorough numerical investigation of the statistics of the two largest eigenvalues and the gap of the spike random matrix $\mathit{M}$, for large but finite system size N.

3. Statistics of the Two Largest Eigenvalues and the Gap for Finite-Size Spike Matrices

In the limit $N\to \infty$, the distribution of eigenvalues is given by (5). When N is finite, the border of the spectrum shows finite-size fluctuations. For spike matrices belonging to the complex Wishart ensemble a phase transition was identified in the behavior of ${\lambda }_1$, as the mean value of the elements of the random matrix changes [19,20]. A similar behavior for real Wishart matrices was conjectured in [19] and subsequently confirmed by several approaches (see, e.g., [23,24] and references therein). Extensions for the GOE and other Gaussian ensembles were considered in [24]. Its connection with the thermodynamic phase transition in the SSK model is immediate because the free energy of the model (which reduces to the average Hamiltonian at zero temperature) is proportional to ${\lambda }_1$. Then, when considering finite-size fluctuations at $T=0$, three regimes are of interest: a sub-critical regime when $\theta \ll J$, a critical one when $\theta \sim J$, and a super-critical one when $\theta \gg J$. The fluctuations of ${\lambda }_1$ in the sub-critical and super-critical regimes have been considered in several works [12,13,21,22,25]. Nevertheless, results on the critical regime are scarce [23,24]. The following results are known: fixing $J=1$, as long as $\theta <1$, the perturbation has little effect on the behavior of ${\lambda }_1$. In this case its distribution is described by the GOE Tracy–Widom (TW) distribution [10,13,25]:

$$N^{-2/3}\rho \left[N^{2/3}({\lambda }_1-\mathbb{E}\left[{\lambda }_1\right])\right]⟹TW,\theta <1.$$

(15)

Instead, when $\theta >1$, ${\lambda }_1$ becomes an isolated eigenvalue as it goes away from the support of the semicircle, being freer to fluctuate around the expected value, $\mathbb{E}\left[{\lambda }_1\right]=\theta +{\theta }^{-1}$. In this case, the fluctuation of ${\lambda }_1$ is of order $\mathcal{O}\left(N^{-1/2}\right)$, described by a normal distribution [13,25]:

$$N^{-1/2}\rho \left[N^{1/2}({\lambda }_1-\mathbb{E}\left[{\lambda }_1\right])\right]⟹\mathcal{N}\left[0,2(1-{\theta }^{-2})\right],\theta >1.$$

(16)

Figure 1 shows the behavior of the probability density distribution of ${\lambda }_1$, collected from an ensemble of ${10}^4$ spike random matrices of size $N=100$, for several values of $\theta$, shown in the color scale to the right of the figure. The distributions are centered at zero, $\overline{{\lambda }_1}$ stands for the ensemble average and ${\sigma }_{\theta }$ is the predicted standard deviation given in (16), applied only in its valid interval $\theta >1$. The Tracy–Widom and the normal distribution, shown in continuous and dashed lines, are properly scaled. The plots of the TW distributions were done by using the publicly available package in https://github.com/yymao/TracyWidom/ (accessed on 10 February 2022). This package uses interpolation tables from [32,33]. In agreement with the results above, for $\theta \ll 1$, the pdf of the largest eigenvalue is well described by a TW distribution. At the other end, when $\theta \gg 1$, a normal distribution with the theoretically predicted behavior is observed. A crossover behavior at intermediate values of $\theta$ is also observed. This is the critical regime. At present, there are a few results on the behavior of the largest eigenvalue of spike random matrices in the critical regime [23,24], from which we have been able to describe the scaling of the expectation value of ${\lambda }_1$ and ${\lambda }_2$, as will be shown later.

3.1. Expectation Value of ${\lambda }_1$

3.1.1. Sub-Critical Regime, $\theta \ll 1$

When $\theta <1$ and for large but finite N, ${\lambda }_1$ is expected to behave as ${\lambda }_1=2+\xi N^{-2/3}$, where $\xi$ is a random variable described by the GOE TW distribution, with expected value $\mathbb{E}\left[\xi \right]=-1.21$ [10,34]. Then, the expected value of ${\lambda }_1$ would behave as $\mathbb{E}\left[{\lambda }_1\right]\approx 2-1.21N^{-2/3}$. Nevertheless, the previous result is valid when the diagonal elements of the random matrix are non-null. In the present case, the matrix $\mathit{M}$ is traceless, with all the diagonal elements equal to zero. In Figure 2a, we can see that the semicircle moves approximately linearly to the left as the perturbation intensity increases, while in Figure 2b it is clear that bigger matrix sizes N suffer smaller shifts. This is a consequence of the traceless character of the matrix. Upon changing the average value of its elements, $\theta /N$, the eigenvalues will have to rescale their expected values in order to satisfy the condition that they must add up to zero. This is analog to a center of mass conservation of the eigenvalue density.

Then, because for finite N the weight of each eigenvalue is $1/N$, the expected value of ${\lambda }_1$ should approximately be given by:

$$\mathbb{E}\left[{\lambda }_1\right]\approx \left(2-1.21N^{-2/3}\right)\left(1-\frac{1}{N}\right),\theta \ll 1.$$

(17)

Note that because this correction acts equally on every eigenvalue, it will have no effect on the gap $g={\lambda }_1-{\lambda }_2$, which is the relevant quantity for the long time dynamics. Figure 3 shows the deviation of the numerical average of ${\lambda }_1$ from the theoretical prediction (17). The improvement of the collapse after the inclusion of the center of mass conservation effect is evident in the right panel. It is also possible to note that the collapse breaks down for $\theta >0.6$ for the sizes considered when the system begins to cross over to the critical regime.

3.1.2. Super-Critical Regime, $\theta \gg 1$

This is the regime in which ${\lambda }_1$ becomes isolated from the bulk. In this case the finite N fluctuations are predicted to be Gaussian, given by Equation (16). As in the $\theta <1$ case, the center of mass conservation must be obeyed. We found that it amounts to a shift of the large N result $(\theta +{\theta }^{-1})$ by the appropriate weight factor $1/N$. Then, for $\theta \gg 1$, the expected value of ${\lambda }_1$ will be given by:

$$\mathbb{E}\left[{\lambda }_1\right]\approx (\theta +{\theta }^{-1})\left(1-\frac{1}{N}\right),\theta \gg 1.$$

(18)

In Figure 4, the effect of the center of mass correction for $\theta >1$ can be appreciated. In this case, besides N, there is a dependence on $\theta$, evident in the left panel of the figure. Upon considering the center of mass correction, the result agrees well with Equation (18), when $\theta \gg 1$.

3.1.3. Critical Regime, $\theta \sim 1$

In reference [24], the statistics of the largest eigenvalue of spike real Gaussian random matrices in the critical regime are considered. The critical regime is defined for fixed values of the parameter $\omega =N^{1/3}(\theta -1)\in (-\infty ,\infty ]$, and was described originally for the spike complex Wishart ensemble in [19], where the phenomenon of the phase transition in the statistics of the largest eigenvalue was identified. In Theorem 1.5 of [24], it is shown that, in the critical regime, the eigenvalues of spike Gaussian random matrices are given in terms of the eigenvalues of the stochastic Airy operator ${\mathcal{H}}_{\beta ,\omega }$ with suitable boundary conditions. For finite $\omega$, the statistics of the largest eigenvalue ${\lambda }_1$ is described by a “one parameter family of deformations of the Tracy–Widom($\beta$)” distributions, interpolating between the usual values of $\beta =1,2,4$. As a consequence, in the critical regime, one expects the ${\lambda }_1$ fluctuations to be approximately described by the TW distribution, but not exactly, with a difference that depends on the value of $\omega$. In Figure 5, we show the (numerical) standard deviation of ${\lambda }_1$ for different system sizes. In the left panel, the raw data are shown as a function of $\theta$. In the right panel, a data collapse is shown, with the scaling variable $\omega$ as defined above, assuming fluctuations to scale with $N^{2/3}$, as would be expected for a perfect Tracy–Widom behavior. Although the collapse is good for $\omega <0$ and performs better in the whole interval as the size N grows, the quality of the collapse decays as $\omega$ grows. According to the results in [24], a continuous change in the exponent, away from $2/3$, should be expected.

3.2. Expectation Value of ${\lambda }_2$

As ${\lambda }_1$ jumps outside the semicircle of the Wigner law, it is expected that ${\lambda }_2$ will take its place at the soft edge of the eigenvalue density function. In particular, it is expected that ${\lambda }_2$ will show fluctuations given by the Tracy–Widom distribution. Then, the expectation value of ${\lambda }_2$ should behave as:

$$\mathbb{E}\left[{\lambda }_2\right]\approx (2-1.21N^{-2/3})\left(1-\frac{1}{N}\right)-\frac{(\theta +{\theta }^{-1})}{N},\theta \gg 1,$$

(19)

where the first term corresponds to $\mathbb{E}\left[{\lambda }_1\right]$ in the sub-critical regime, Equation (17), and the second is the correction due to the center of mass conservation when the largest eigenvalue has detached from the bulk. The behavior of the shift of the numerical average ${\overline{\lambda }}_2$ from the expectation given by Equation (19) is shown in Figure 6. In the left panel, only the first term on the right-hand side of (19) is shown, while the right panel shows the full expression after taking into account the center of the mass conservation term. A progressive good collapse can be seen in the $\theta \gg 1$ regime, as N grows.

In Figure 7, we show a data collapse of the fluctuations of ${\lambda }_2$. The collapse is good for the largest sizes, in agreement with theoretical expectations.

3.3. Statistics of Small Gaps $g={\lambda }_1-{\lambda }_2\ll 1$

The behavior of the gap between the two largest eigenvalues, $g={\lambda }_1-{\lambda }_2$, will also depend on the regime considered. In the sub-critical regime, the effect of the deterministic perturbation is negligible, and one expects that the statistics of the gap will be governed by the results of reference [35]. In turn, this will lead to power law time/size scalings, as studied in [16,17]. In the super-critical regime, the two largest eigenvalues become approximately independent random variables. From Equations (15) and (16), the Tracy–Widom character of the fluctuations of ${\lambda }_2$, going to zero as $N^{-2/3}$, are sub-leading with respect to the Gaussian fluctuations of ${\lambda }_1$, which decay as $N^{-1/2}$. Then, in the regime $\theta \gg 1$, one can expect (and we have numerically verified) the fluctuations of small gaps to be Gaussian, like those of ${\lambda }_1$. In turn, this will reflect in exponential time relaxations of observables, such as the energy gap, a typical behavior of ferromagnetic phases. More interesting is the intermediate, critical regime. In this regime, in which ${\lambda }_1$ and ${\lambda }_2$ are strongly correlated, the statistical behavior of the gap is not known. With the aim of describing the long time behavior of the energy gap, we have pursued a numerical characterization of the statistics of the gap in the small gap regime, $g\ll 1$, relevant to the long time relaxation dynamics. In Figure 8, the distribution of the gap is shown in double logarithmic scale for an ensemble of random spike matrices of size $N=1000$ and different values of $\theta \ge 1$. In all cases, it can be seen that the behavior is algebraic for small g.

Thus, for the small gaps regime, we expect that the pdf of the gap will approximately behave as:

$$f\left(g\right)\sim b(\theta ,N)g^{a(\theta ,N)},$$

(20)

where $a(\theta ,N)$ and $b(\theta ,N)$ are parameters to be determined. We numerically adjusted those parameters to fit the data in the interval $g<\mathbb{E}\left[g\right]-{\sigma }_g$. In Figure 9 we show the behavior of the gap exponent with $\theta$, for different system sizes. The left panel shows that $a(\theta ,N)$ is a constant equal to one in the sub-critical regime, in agreement with the results for the pure SSK model [35]. For $\theta >1$, the numerical analysis suggests a linear behavior, with a slope dependent with N. A very good data collapse is obtained as a function of the critical scaling variable $N^{1/3}(\theta -1)$, as can be seen in the right panel of Figure 9. Then, for $\theta >1$, the gap exponent behaves approximately as:

$$a(\theta ,N)\approx 1+c_1{N}^{1/3}(\theta -1)-c_2,$$

(21)

where $c_1\approx 2.1$ and $c_2\approx 5.5$ are fit parameters.

Figure 10 shows the behavior of the gap parameter $b(\theta ,N)$. In the left panel, it can be seen that it presents a minimum at $\theta >1$, which moves towards $\theta =1$ as the system size grows. Then, the critical sector can be identified as the region around the minimum. A good finite-size collapse around the minima is obtained with

$$b(\theta ,N)=N^{4/3}h\left[N^{1/3}(\theta -1)\right],$$

(22)

where $h\left(x\right)$ is an unknown scaling function, as shown in the right panel of the figure.

With these results, it is now possible to compute the finite-size behavior of the time-dependent excess energy, which is the subject of the next section.

4. Long Time Decay of the Excess Energy

4.1. $N\to \infty$ Limit

For random initial conditions given by (13), the Lagrange multiplier is given by [16,17]:

$$z\left(t\right)=\frac{1}{2}\frac{d}{dt}\ln \left[\frac{1}{N}\sum _{\mu =1}^N{e}^{2{\lambda }_{\mu }t}\right]\stackrel{N\to \infty }{\to }\frac{1}{2}\frac{d}{dt}ln\left[{\int }_{-\infty }^{\infty }d\lambda \rho \left(\lambda \right)e^{2\lambda t}\right],$$

(23)

where the density of eigenvalues, $\rho \left(\lambda \right)$, is given by (5). Performing the integrations, the exact solution is given by:

$$z\left(t\right)=\frac{I_1\left(4t\right)}{2t}+e^{2(\theta +{\theta }^{-1})t},\theta >1,$$

(24)

where $I_1\left(x\right)$ is a modified Bessel function of the first kind. In the long time regime, the above expression has the asymptotic behavior:

$$z\left(t\right)\to \left(\theta +{\theta }^{-1}\right)+\left[2-\left(\theta +{\theta }^{-1}\right)\right]\frac{e^{[4-2(\theta +{\theta }^{-1})]t}}{4\sqrt{2\pi }{t}^{3/2}},\theta >1.$$

(25)

Remembering that the Lagrange multiplier is proportional to the energy density, $z\left(t\right)=-2e\left(t\right)$ and that, when $\theta >1$, the largest eigenvalue is given by ${\lambda }_1=\theta +{\theta }^{-1}$, we obtain for the long time behavior of the average excess energy the result:

$$\underset{N\to \infty }{lim}\mathbb{E}[\Delta e(t,\theta )]=\left\{\begin{array}{cc}\epsilon \left(t\right)\hfill & \mathrm{for}\theta \le 1\hfill \\ \frac{\left[\left(\theta +{\theta }^{-1}\right)-2\right]}{8\sqrt{2\pi }}\frac{e^{[4-2(\theta +{\theta }^{-1})]t}}{t^{3/2}}\hfill & \mathrm{for}\theta >1,\hfill \end{array}\right.$$

(26)

where $\epsilon \left(t\right)=3/8t$ is the known result for the pure SSK model [8]. We note that, as expected, in the $\theta >1$ ferromagnetic regime the asymptotic relaxation, when $t\gg {[4-2(\theta +{\theta }^{-1})]}^{-1}$, is exponential, faster than in the spin glass phase of the model.

4.2. Finite System Size

With the results obtained for the pdf of the gap in the small gap regime for finite N, we can compute the late time behavior of the average excess energy, as given by Equation (14):

$$\mathbb{E}[\Delta e(t,\theta ,N)]={\int }_0^{\infty }g{e}^{-2gt}f\left(g\right)dg.$$

(27)

The exact behavior of $f\left(g\right)$ for large gaps is not known. Nevertheless, an inspection of the numerical behavior of $f\left(g\right)$ (see, e.g., Figure 8) shows that there is a scale $r_N$ such that for $g\gg r_N$ the gap distribution density shows a fast monotonic decrease towards zero. Then we can write:

$$\begin{array}{ccc}\hfill \mathbb{E}[\Delta e(t,\theta ,N\left)\right]& \approx & b(\theta ,N){\int }_0^{r_N}dgg{e}^{-2gt}{g}^{a(\theta ,N)}+{\int }_{r_N}^{\infty }dgg{e}^{-2gt}h(g/r_N)\hfill \\ & =& b(\theta ,N)t^{-(2+a)}{\int }_0^{r_{N}t}dx{x}^{1+a}{e}^{-2x}+r_N^2{\int }_1^{\infty }dxx{e}^{-2\left(r_{N}t\right)x}h\left(x\right),\hfill \end{array}$$

(28)

where $h\left(x\right)$ is a monotonic decreasing function of its argument. At long times, $r_{N}t\gg 1$, the second integral becomes negligible, and then:

$$\mathbb{E}[\Delta e(t,\theta ,N)]\approx b(\theta ,N)t^{-(2+a)}\left[2^{-(2+a)}\left(\Gamma (2+a)-\Gamma (2+a,2r_{N}t)\right)\right],$$

(29)

with the limit ${lim}_{r_{N}t\to \infty }\Gamma (2+a,2r_{N}t)=1$. From (22), (29) and the behavior observed in the numerical results of Figure 11, the long time behavior of the average excess energy, in the critical regime, can be described as follows:

$$\mathbb{E}[\Delta e(t,\theta ,N)]\sim \left\{\begin{array}{cc}\upsilon \left(t\right)\hfill & \mathrm{for}t<N^{1/3}\hfill \\ N^{1/3}{f}_a\left(t{N}^{-1/3}\right)\hfill & \mathrm{for}t>N^{1/3}\hfill \end{array}\right.$$

(30)

The function $\upsilon \left(t\right)$ represents the time relaxation of the excess energy in the thermodynamic limit, given by (26) when $\theta >1$ and $f_a\left(x\right)=c/x^{2+a}$ is a one-parameter scaling function, with a given by (21) and c a constant. It is to be noted the difference between the previous results and those for the pure SSK model [16,17]. When $\theta >1$ and $t<N^{1/3}$, the relaxation is exponential, given by $\upsilon \left(t\right)$ instead of the power law in the pure case. Additionally, the scaling exponent $1/3$ of the algebraic regime differs from the $2/3$ of the pure SSK model. In the limit $a\to 1$, the expression leads to the expected behavior for the pure model $\mathbb{E}[\Delta e]\sim N^{4/3}/t^3$ [16,17].

In Figure 11, we can see data collapses of the average excess energy of the spike SSK model, for different fixed values of the parameter $a(\theta ,N)$. They show good agreement with the previous results. Similarly to the pure SSK model, there is an algebraic scaling regime not present in the $N\to \infty$ limit. In this case, the exponent $2+a(\theta ,N)$ depends on both the spike intensity $\theta$ and the size of the system N, according to (21). In the figures, particular combinations of $\theta$ and N were chosen in order to keep the value of $a(\theta ,N)$ approximately constant. As the exponent grows, so does the slope of the power law. The relaxation becomes faster, but a power law regime can be identified by values of the parameter a as large as $a=6$, as shown in the figure.

5. Discussion and Conclusions

In the first part of this work, we have presented a numerical study of the statistics of the two largest eigenvalues and the gap for random matrices from the Gaussian Orthogonal Ensemble perturbed by a deterministic rank one matrix. The largest eigenvalue of such spike random matrices is known to go through a phase transition as the intensity of the deterministic perturbation attains a critical value [19,24]. The statistics in the sub-critical and super-critical regimes are well described in the literature, while results for the critical regime are scarce [23,24]. Our numerical analysis of the average values of the two largest eigenvalues confirmed the analytical results from the literature, after the inclusion of additional size effects due to the traceless character of the matrices considered in this work; these results are of interest in physics models. In the critical regime, we showed results on the fluctuations of ${\lambda }_1$ and ${\lambda }_2$ in good agreement with available theoretical predictions on the existence of a critical scaling regime where the fluctuations are described by a one-parameter family of scaling functions, which can be seen as continuous deformations of the Tracy–Widom distribution [23,24]. Although our results are compatible with that conclusions, more work is needed to extend and clarify the interpretation of the results from the mathematical literature in the physical models context.

For the statistics of the gap, at present, there are no known analytical results. Then, we pursued a numerical characterization of the small gap regime of the probability density function, which is the relevant regime for the long time behavior of physical observables. We show evidence that the pdf of the gap has a power law behavior for small gaps. The exponent of the power law depends on both the intensity of the deterministic perturbation $\theta$ and the system size N, in the form which defines the critical sector of the model, $a(\theta ,N)\sim N^{1/3}(\theta -1)$, when $\theta >1$. After the characterization of the pdf of the gap, we described the long time decay of the average excess energy of the Spherical Sherrington–Kirkpatrick model with a Curie–Weiss perturbation term. We first obtained the analytical result in the large N limit, showing that, as expected, the relaxation is exponential for $\theta >1$. We then considered the large but finite N behavior. The most interesting and new regime to describe is near the phase transition between the spin glass and ferromagnetic phases. In this critical regime, using the results obtained for the gap pdf, we showed the existence of a sector with power law relaxation as a function of the scaling variable $t{N}^{-1/3}$.

The results for the gap pdf and the excess energy relaxation are our main new results, not previously reported in the literature. Being mainly of a numerical character, we expect that they will motivate us to pursue analytical approaches to the computation of the gap probability distribution function, which has been shown to be a relevant random variable to describe the late time dynamics of spherical models with pairwise interactions.