Finite-size relaxational dynamics of a spike random matrix spherical model

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 is characterized by the two largest eigenvalues of a spike random matrix which defines the model, and particularly by the statistics of the gap between them. We characterize the finite size statistics of the two largest eignevalues 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 strenght of the non-disordered perturbation, in a way which is governed by the finite size statistics of the gap.


I. INTRODUCTION
Quenched random interactions are known to be at the origin of complex behavior in many body systems, both in the thermodynamics and dynamics as well [1][2][3]. The understanding of the properties of this kind of systems in the thermodynamic limit has been steadily growing in the last 40 years or so, mainly through 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], neural networks [6], to cite but a few. Powerful techniques, like the saddle point method, are not so useful for studying systems far from the thermodynamic limit. There are a few class of models in which both the thermodynamics and the dynamical behavior as well can be solved exactly and still they show interesting non-trivial properties qualitatively similar to more complex systems. Well known examples are 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 was 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, like the phenomenon of aging, is present in the model and were completely characterized. The solution of both the thermodynamics and the dynamics of the model were 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 × 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 β distributions [9,10], where β = 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]. 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 [15,16]. In these works, a new algebraic (power law) scaling regime was found and characterized, not present in the N → ∞ regime. In the dynamical context, besides the relevance of the largest eigenvalue, which is directly proportional to the ground state energy, also the gap between the two largest eigenvalues is a fundamental quantity to compute the relevant time/size scalings of the long time relaxation. A related interesting 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 was 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, this kind 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 [17][18][19][20][21][22][23][24]. Of special interest for the physics community is the result, originally presented in [18], 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 in different applications, e.g overlaps between eigenvectors of correlated spike random matrices [25], analysis of optimal learning rates in non-convex optimization [26], low-rank matrix estimation [27], limits of detection of planted states [28], ruggedness of complex energy landscapes [29]. From a dynamical perspective, understanding of spike random matrix models may shed light on problems like the feasibility of identifying a deterministic signal in a random environment, the reconstruction of hidden patterns in a complex landscape or the efficency 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 II we introudce the model studied and summarize some known results which will be useful later; in Section III 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 IV 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 V, we make a brief discussion of the work and present our conclusions.

II. THE MODEL
The spherical Sherrington-Kirkpatrick (SSK) model with a Curie-Weiss (CW) perturbation is described by the following Hamiltonian: where z is a Lagrange multiplier which enforces the spherical constraint: In the previous expressions, the spin variables s i ∈ [− √ N , √ N ] are described by a N component vector ⃗ S = (s 1 , ..., s N ). The coupling constants J ij are chosen form a real and symmetric random matrix from the Gaussian Orthogonal Ensemble (GOE), J = {J ij } (i,j)∈[1,N ] 2 , with zero mean and variance J 2 . θ ∈ R + measures the intensity of the deterministic perturbation. Thus, M is a real symmetric N × N spike random matrix, whose off-diagonal elements are Gaussian distributed with the following mean and variance: The diagonal elements are zero. The Hamiltonian can be rewritten decomposing ⃗ S as a linear combination of the eigenvectors { ⃗ V µ } of the coupling matrix M , with ⃗ V µ · ⃗ V ν = δ µν . Therefore, with the following notation s µ = ⃗ S · ⃗ V µ for the projections of ⃗ S on the eigenvectors of M , the Hamiltonian becomes: with {λ µ } µ∈ [1,N ] being the set of N eigenvalues of M , with associated eigenvectors { ⃗ V µ }. The eigenvalues are organised such that Max[λ µ ] = λ 1 > λ 2 > ... > λ N = Min[λ µ ]. In the large N limit, the eigenvalue density distribution of M is given by [30]: where ρ W (λ) is the Wigner semicircle law: The result (5) means that, if θ ≤ J, the spectrum of M is given by the Wigner law, corresponding to the GOE ensemble. Otherwise, if θ > J, the largest eigenvalue λ 1 detaches from the Wigner semicircle, becoming an outlier with a delta peak at λ = θ + J 2 θ −1 .
The overdamped dynamics of the model is governed by the set of Langevin equations : where ξ i (t) represents a Gaussian white noise with zero mean and variance ⟨ξ i (t)ξ i (t ′ )⟩ = 2T δ ij δ(t − t ′ ) and T is the temperature of a thermal bath. As in [16], here we are interested in the zero temperature limit of the Langevin equations which, in the eigenbasis of M , reads: 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 H[ ⃗ S, z]. The asymptotic stationary state will depend on the initial conditions {s µ (0)}, as stated explicitely in (8). By setting lim t→∞ ∂ t s µ (t) = 0 we obtain the criteria: complemented by the spherical constraint in eq. (2). This system of equations admits the 2N solutions: Their stability is determined by the Hessian Taking a given metastable state ⃗ S = √ N ⃗ V µ , the local landscape has N − µ stable directions, µ − 1 unstable directions and a marginal flat one. The energy of each of these configurations is equal to −λ µ N/2. Thus, the system should always equilibrate in one of the solutions ± √ N ⃗ V 1 , as they are the only stable ones. The ground state energy density is then simply given by e eq = −λ 1 /2. Our primary interest here is to describe the behaviour, at long times, of the excess energy density, ∆e(t, N ) = e(t, N ) − e eq (N ), for arbitrary system sizes N . To this end, recalling results in [15,16], 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: Because of the dependence of the above expression on the relaxation rates λ µ − λ 1 , one expects that the late dynamics of the model will be dominated by the gap, g = λ 1 − λ 2 , between the two largest eigenvalues of the random matrix M : In the previous expression there are two sources of fluctuations: the statistics of the gap and the initial conditions. The vector of initial conditions ⃗ S(0) = {s µ (0)} can be written in the basis of eigenvectors of the M matrix in the form ⃗ In the present work, we are primarily interested in a flat distribution on the basis of eigenvectors, that is c ν = 1 for all ν, which can be associated to thermal equilibrium at a very high temperature. It corresponds to: The relaxation dynamics of the model depends on the form of the eigenvalue density (5). The effect of the delta contribution is to induce a phase transtion when the intensity of the Curie-Weiss term attains the value θ = J in the thermodynamic limit. While the CW term remains weaker than the random couplings intensity, θ < J, the system behaves like the pure SSK model, relaxing towards a disordered ground state ⃗ S with a characteristic slow dynamics, as described in [8,15,16]. 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, θ > 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 θ = J, the system goes through a continuous phase transtion between a disordered spin glass phase and a ferromagnetically ordered one, in the thermodynamic limit [7]. In [12] finite size fluctuations of the free energy of the model at both sides of the spin glass-ferromagnetic transition where 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 givem by (13), at long times, the behavior of the average excess energy is given by: At present, the statistical properties of the gap g = λ 1 − λ 2 are not known. In order to describe its approximate behaviour, in the following we will pursue a thorough numerical investigation of the statistics of the two largest eigenvalues and the gap of the spike random matrix M , for large but finite system size N .

III. STATISTICS OF THE TWO LARGEST EIGENVALUES AND THE GAP FOR FINITE SIZE SPIKE MATRICES
In the limit N → ∞ 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 λ 1 , as the mean value of the elements of the random matrix changes [18,19]. A similar behavior for real Wishart matrices was conjectured in [18] and subsequently confirmed by several approaches (see e.g. [22,23] and references therein). Extensions for the GOE and other Gaussian ensembles were considered in [23]. 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 λ 1 . Then, when considering finite size fluctutations at T = 0, three regimes are of interest: a sub-critical regime when θ ≪ J, a critical one when θ ∼ J and a super-critical one when θ ≪ J. The fluctuations of λ 1 in the sub-critical and super-critical regimes have been considered in several works [11,12,20,21,24]. Nevertheless, results on the critical regime are scarce [22,23]. The following results are known: fixing J = 1, as long as θ < 1, the perturbation has little effect on the behavior of λ 1 . In this case its distribution is described by the GOE Tracy-Widom (TW) distribution [10,12,24]: Instead, when θ > 1, λ 1 becomes an isolated eigenvalue as it goes away from the support of the semicircle, being freer to fluctuate around the expected value, E[λ 1 ] = θ + θ −1 . In this case, the fluctuation of λ 1 is of order O(N −1/2 ), described by a normal distribution [12,24]: Figure 1 shows the behavior of the probability density distribution of λ 1 , collected from an ensemble of 10 4 spike random matrices of size N = 100, for several values of θ, shown in color scale to the right of the figure. The distributions are centered at zero, λ 1 stands for the ensemble average and σ θ is the predicted standard deviation given in (16), applied only in its valid interval θ > 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/. This package uses interpolation tables from [31,32].
In agreement with the results above, for θ ≪ 1 the pdf of the largest eigenvalue is well described by a TW distribution. At the other end, when θ ≫ 1 a normal distribution with the theoretically predicted behavior is observed. It is also observed a crossover behavior at intermediate values of θ. This is the critical regime. At present, there are a few results on the behaviour of the largest eigenvalue of spike random matrices in the critical regime [22,23], from which we have been able to describe the scaling of the expectation value of λ 1 and λ 2 , as will be shown later. When θ < 1 and for large but finite N , λ 1 is expected to behave as λ 1 = 2 + ξ N −2/3 , where ξ is a random variable described by the GOE TW distribution, with expected value E[ξ] = −1.21 [10,33]. Then, the expected value of λ 1 would behave as E[λ 1 ] ≈ 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 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, θ/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 λ 1 should approximately be given by: Note that, because this correction acts equally on every eigenvalue, it will have no effect in the gap g = λ 1 − λ 2 , which is the relevant quantity for the long time dynamics. Figure 3 shows the deviation of the numerical average of λ 1 from the theoretical prediction (17). The improvement of the collapse after 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 θ > 0.6, for the sizes considered, when the system begins to cross over to the critical regime.

Super-critical regime, θ ≫ 1
This is the regime in which λ 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 θ < 1 case, the center of mass conservation must be obeyed. We found that it amounts to a shift of the large N result (θ + θ −1 ) by the appropriate weigth factor 1/N . Then, for  (17), for θ < 1. The left panel shows results without considering the correction due to the conservation of the center of mass. In the right panel, after inclusion of the correction, the data shows a good collapse for growing sizes N and sufficiently small values of θ. θ ≫ 1, the expected value of λ 1 will be given by: In Figure 4 the effect of the center of mass correction for θ > 1 can be appreciated. In this case, besides N , there is a dependence on θ, evident in the left panel of the figure. Upon considering the center of mass correction, the result agrees well with equation (18), when θ ≫ 1.

Critical regime, θ ∼ 1
In reference [23] the statistics of the largest eigenvalue of spike real Gaussian random matrices in the critical regime is considered. The critical regime is defined for fixed values of the parameter ω = N 1/3 (θ − 1) ∈ (−∞, ∞], and was described originally for the spike complex Wishart ensemble in [18], where the phenomenon of the phase transition in the statistics of the largest eigenvalue was indentified. In Theorem 1.5 of [23] 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 H β,ω with suitable boundary conditions. For finite ω, the statistics of the largest eigenvalue λ 1 is described by a "one parameter family of deformations of the Tracy-Widom(β)" distributions, interpolating between the usual values of β = 1, 2, 4. As a consequence, in the critical regime, one expects the λ 1 fluctuations to be approximately described by the TW distribution, but not exactly, with a difference that depends on the value of ω. In Figure 5 we show the (numerical) standard deviation of λ 1 for different system sizes. In the left panel the raw data is shown as a function of θ. In the right panel a data collapse is shown, with the scaling variable ω as defined above, assuming fluctuations to scale with N 2/3 , as would be expected for a perfect Tracy-Widom behavior. While the collapse is good for ω < 0 and performs better in the whole interval as the size N grows, the quality of the collapse decays as ω grows. According to the results in [23], a continuous change in the exponent, away from 2/3, should be expected.

B. Expectation value of λ2
As λ 1 jumps outside the semicircle of the Wigner law, it is expected that λ 2 will take its place at the soft edge of the eigenvalue density function. In particular, it is expected that λ 2 will show fluctuations given by the Tracy-Widom distribution. Then, the expectation value of λ 2 should behave as: where the first term corresponds to E[λ 1 ] 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 λ 2 from the expectation given by eq. (19) is shown in Figure 6. In the left panel only the first term on the righthand side of (19) is shown, while the right panel shows the full expression, after taking into account the center of mass conservation term. A progressive good collapse can be seen, in the θ ≫ 1 regime, as N grows. In Figure 7 we show a data collpase of the fluctuations of λ 2 . The collapse is good for the largest sizes, in The behavior of the gap between the two largest eigenvalues, g = λ 1 −λ 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 [34]. In turn, this will lead to power law time/size scalings, as studied in [15,16]. In the super-critical regime, the two largest eigenvalues become approximately independent random variables. From eq. (16), λ 1 shows Gaussian flucutations which grow with θ. Then, in this regime, the fluctuations of the gap are expected to be Gaussian also, similar to what happens with λ 1 . In turn, this will reflect in exponential time relaxations of observables, like the energy gap, a typical behavior of ferromagnetic phases. More interesting is the intermediate, critical regime. In this regime, in which λ 1 and λ 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 ≪ 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 θ ≥ 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: where a(θ, N ) and b(θ, N ) are parameters to be determined. We numerically adjusted those parameters to fit the data in the interval g < E[g] − σ g . In Figure 9 we show the behavior of the gap exponent with θ, for different system sizes. The left panel shows that a(θ, N ) is a constant equal to one in the sub-critical regime, in agreement with the results for the pure SSK model [34]. For θ > 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 (θ − 1), as can be seen in the right panel of Figure 9. Then, for θ > 1, the gap exponent behaves approximately as: where c 1 ≈ 2.1 and c 2 ≈ 5.5 are fit parameters. With these results, now it is possible to compute the finite size behavior of the time dependent excess energy, which is the subject of the next section.

IV. LONG TIME DECAY OF THE EXCESS ENERGY
For random initial conditions given by (13), the Lagrange multiplier is given by [15,16]: where the density of eigenvalues, ρ(λ), is given by (5). Performing the integrations, the exact solution is given by: where I 1 (x) is a modified Bessel function of the first kind. In the long time regime, the above expression has the asymptotic behavior: Remembering that the Lagrange multiplier is proportional to the energy density, z(t) = −2e(t) and that, when θ > 1, the largest eigenvalue is given by λ 1 = θ + θ −1 , we obtain for the long time behavior of the average excess energy the result: where ε(t) = 3/8t is the known result for the pure SSK model [8]. We note that, as expected, in the θ > 1 ferromagnetic regime the asymptotic relaxation, when t ≫ [4 − 2(θ + θ −1 )] −1 , is exponential, faster than in the spin glass phase of the model.

B. 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 eq. (14): which can be decomposed in the form ∞ 0 = r 0 + ∞ r . In the long time limit the dynamics will be dominated by the sector of small values of g. By a similar analysis to that presented in equations (42)-(43) of [16], it follows that the second integral will be negligible. Then, using (20): which has the exact solution: with the limit lim rt→∞ Γ(2 + a, 2rt) = 1. From the analytical results and numerical analysis, we found that the behavior of the average excess energy, in the regime θ > 1, can be described as follows: The function υ(t) represents the time relaxation of the excess energy in the thermodynamic limit, given by (25) when θ > 1 and f a (x) = 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 [15,16]. When θ > 1 and t < N 1/3 the relaxation is exponential, given by υ(t), instead of the power law in the pure case. Also, the scaling exponent 1/3 of the algebraic regime differs from the 2/3 of the pure SSK model.
In Fig. 10 we can see data collapses of the average excess energy of the spike SSK model, for different fixed values of the parameter a(θ, 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 → ∞ limit. In this case, the exponent 2 + a(θ, N ) depends on both the spike intensity θ and the size of the system N , according to (21). In the figures, particular combinations of θ and N were chosen in order to keep the value of a(θ, 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 = 7, as shown in the figure.  (θ, N ). The slope of the algebraic regime is governed by the one parameter scaling function f a (x).

V. 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 [18,23]. The statistics in the sub-critical and super-critical regimes are well described in the literature while results for the critical regime are scarce [22,23]. Our numerical analysis on the average values of the two largest eigenvalues confirmed analytical results from the literature, after inclusion of additional size effects due to the traceless character of the matrices considered in this work, which are of interest in physics models. In the critical regime, we showed results on the fluctuations of λ 1 and λ 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 [22,23]. While 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 θ and the system size N , in the form which defines the critical sector of the model, a(θ, N ) ∼ N 1/3 (θ − 1), when θ > 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 θ > 1. We then considered the large but finite N behavior. The most interesting and new regime to describe is near the phase transtion 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 tN −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 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.