11institutetext: Department of Theoretical Physics, Wrocław University of
Science and Technology, 50-370 Wrocław, Poland 22institutetext: Institut für
Theoretische Physik, Goethe Universität, Max-von-Laue-Str. 1, 60438 Frankfurt
am Main, Germany 33institutetext: Asia Pacific Center for Theoretical Physics,
Pohang, 37673, Korea 44institutetext: Department of Physics, Pohang University
of Science and Technology, Pohang, 37673, Korea 55institutetext: Frankfurt
Institute for Advanced Studies, Ruth-Moufang-Str. 1,60438 Frankfurt, Germany
66institutetext: School of Mathematics, Trinity College, Dublin 2, Ireland

# Exploring the Phase Diagram of V-QCD with Neutron Star Merger Simulations

Tuna Demircik 11 Christian Ecker 22<EMAIL_ADDRESS>Matti Järvinen
3344 Luciano Rezzolla 225566 Samuel Tootle 22 Konrad Topolski 22

###### Abstract

Determining the phase structure of Quantum Chromodynamics (QCD) and its
Equation of State (EOS) at densities and temperatures realized inside neutron
stars and their mergers is a long-standing open problem. The holographic V-QCD
framework provides a model for the EOS of dense and hot QCD, which describes
the deconfinement phase transition between a dense baryonic and a quark matter
phase. We use this model in fully general relativistic hydrodynamic (GRHD)
simulations to study the formation of quark matter and the emitted
gravitational wave signal of binary systems that are similar to the first ever
observed neutron star merger event GW170817.

APCTP Pre2022 - 026

##  1 Introduction

Multi-messenger observations from binary neutron-star mergers are expected to
provide new insights into properties of dense and hot Quantum Chromodynamics
(QCD) that are currently inaccessible with first-principle methods in QCD and
collider experiments on earth. A key quantity to establish the connection
between the properties of dense QCD matter and neutron stars is the Equation
of State (EOS), which in the zero temperature limit determines the pressure as
a function of the energy density as it is required to solve for the mass-
radius relation of static and cold stars in General Relativity. Robust
predictions for the EOS are only available from chiral effective theory (CET)
calculations Hebeler:2013nza ; Gandolfi2019 ; Keller:2020qhx ;
Drischler:2020yad  at densities below and close to the nuclear saturation
density $n_{s}=0.16\,\rm fm^{-3}$ and from perturbative QCD Fraga2014 ;
Gorda:2021kme ; Gorda:2021znl  at densities that are much larger than those
reached in even the most massive neutron stars. At densities a few times
larger than $n_{s}$, which are realized in the core regions of neutron stars
and their mergers, the available options are either model-building or model
agnostic parameterizations of the EOS. The latter approach, when combined with
theory input from CET and perturbative QCD together with observational data of
neutron star masses and radii NANOGrav:2019jur ; Fonseca:2021wxt ;
Riley:2019yda ; Miller:2019cac ; Nattila:2017wtj ; Miller:2021qha ;
Riley:2021pdl  and bounds on the tidal deformability derived from direct
gravitational wave detections of binary neutron star mergers (Abbott2018a, ),
has lead to a number of important insights on the cold EOS at neutron star
densities Annala2019 ; Annala:2022 ; Altiparmak:2022  and on the properties of
isolated neutron stars (see Ecker:2022 ; Ecker:2022dlg  and references
therein). However, in neutron star mergers, where strong shocks during and
after the collision can locally heat matter to temperatures $T>50$ MeV, the
zero temperature approximation of the EOS becomes insufficient, which leaves
model-building as the only feasible option.

An important open question is whether a phase transition between confined
hadrons and deconfined quark matter can happen during and after the collision
of two neutron stars. It is for example expected that imprints of the
deconfinement phase transition might be visible in the kHz gravitational-wave
signal emitted during a possible post-merger hyper-massive neutron star (HMNS)
stage that encodes information about the hot and dense part of the EOS that is
inaccessible during the inspiral phase Bauswein:2018bma ; Weih:2019xvw .
Addressing the question of the detectability of quark matter in binary
neutron-star mergers requires two important ingredients. The first one is a
model for dense matter that is able to make predictions for the phase
structure of QCD at all temperatures and densities realized in neutron-star
mergers and that can consistently describe the deconfinement phase transition.
In addition, the EOS derived from such a model needs to be consistent with the
aforementioned theoretical constraints obtained from CET and perturbative QCD
in their respective regimes of validity, as well as with available
observational data of neutron stars and their mergers. The second crucial
ingredient is a numeric framework that is able to use such a EOS model as an
input and that can solve complicated GRHD equations of a neutron star merger
as well as extract the gravitational waves produced in such events. In the
following, we shall use EOSs that were constructed with the recently developed
finite-temperature framework based on the holographic V-QCD model
Demircik:2021zll  together with a set of state-of-the-art general relativistic
hydrodynamics codes Most2019b ; Papenfort2021b  to simulate the formation of
quark matter in the HMNS phase of binary neutron star mergers that resemble
the first-ever detected event GW170817.

This proceedings contribution summarises material published in
Demircik:2021zll ; Tootle:2022pvd  and was presented at The XVth Quark
confinement and the Hadron spectrum conference.

##  2 The V-QCD Equation of State Framework

Figure 1: Left: Cold, beta equilibrium slices of the V-QCD EOS. Shown are the
soft (red), intermediate (green) and stiff (blue) versions of the model.
Right: Corresponding mass-radius relations where dashed lines mark unstable
quark matter branches. See Tootle:2022pvd  for details.

In this section we briefly summarise the V-QCD framework for the hybrid EOS
model used in our merger simulations (for details see Jarvinen:2021jbd ;
Demircik:2021zll ). The important novel component of this construction is the
use of a gauge/gravity duality calculation for the dense nucleonic phase and
its transition to deconfined quark matter. The gauge/gravity duality
Maldacena:1997re ; Witten:1998qj  allows one to translate intractable problems
in quantum field theory at strong coupling and finite density to tractable
problems in classical five-dimensional gravity and has been recently applied
by various groups to neutron-star physics Hoyos:2016zke ; Ghoroku:2019trx ;
BitaghsirFadafan:2019ofb ; Ecker:2019xrw ; Demircik:2020jkc ; Kovensky:2020xif
; Pinkanjanarod:2020mgi ; Mamani:2020pks ; Ghoroku:2021fos ; Jarvinen:2021jbd
; Hoyos:2021uff ; Kovensky:2021kzl ; Bartolini:2022rkl .

V-QCD is a holographic bottom-up model with a number of free parameters that
are tuned to reproduce characteristic features of QCD, such as confinement and
chiral symmetry breaking and are further constrained by demanding consistency
with lattice QCD data and astrophysical measurements Jarvinen:2021jbd . The
V-QCD action consists of a gluon part that is given by improved holographic
QCD Gursoy:2007cb ; Gursoy:2007er , a quark (flavour) part implemented with a
tachyonic Dirac-Born-Infeld action Bigazzi:2005md ; Casero:2007ae  and a
component for homogeneous nuclear matter derived in Ishii:2019gta . The so-
called Veneziano limit is assumed in which both the number of colours $N_{c}$
and flavours $N_{f}$ is large but their ratio is $\mathcal{O}(1)$
Jarvinen:2011qe  in order to maintain the backreaction between quarks and
gluons in the holographic large $N_{c}$ limit. This calculation is then
extended to densities around and below $n_{s}$, where the homogeneous
approximation for V-QCD nuclear matter becomes unreliable, with a combination
of the traditional nuclear theory APR Akmal:1998cf  and HS(DD2) Hempel2010 ;
Typel:2009sy  EOSs including their temperature dependence up to $T\approx
160\leavevmode\nobreak\ $MeV, as well as with bosonic contributions given by a
black body photon gas and all mesons of mass $\leq 1$ GeV from the particle
data group listings ParticleDataGroup:2020ssz . The latter contributions are
important to model the QCD critical end point of the deconfinement transition.
The finite temperature extension of the cold V-QCD baryon sector is realised
by matching it with a van der Waals hadron gas model Rischke:1991ke , whose
details can be found in Demircik:2021zll . Finally, the quark fraction and the
phase boundaries of the mixed phase at the first order deconfinement
transition is implemented with a Maxwell construction. In this way the
combined model is able to provide predictions for the finite temperature EOS
on the part of the QCD phase diagram that is relevant to binary neutron star
mergers, core-collapse supernovae and heavy-ion collision experiments.

We employ three variants of the EOS, which were introduced in Demircik:2021zll
in order to represent the parameter dependence of the model. These are called
the soft, intermediate, and stiff variants; and reflect three different
choices of the parameters in the action of V-QCD Jokela:2018ers . In Fig. 1
(left) we show cold beta-equilibrium slices of these three EOSs with
uncertainty bands from nuclear theory Hebeler:2013nza  (blue) and perturbative
QCD Fraga2014  (orange). Red, green and blue lines are the soft, intermediate
and stiff versions of V-QCD where the dotted part of these curves mark the
first order phase transition between the baryonic (solid) the quark matter
(dashed) phase. In addition we show markers for the central densities reached
in various isolated non-rotating stars that we use to initialise the binary
systems in our simulations. Light blue and pink regions mark the residual
freedom of model agnostic parametrizations of the EOS Annala2019 ; Annala:2022
and of the V-QCD model Jokela:2020piw , respectively, after imposing
constraints from neutron star observations. In Fig. 1 (right) we show the
corresponding mass-radius relations of non-rotating stars together with
various error bands of the direct mass measurement (green area) of the heavy
pulsar PSR J0740+6620 NANOGrav:2019jur ; Fonseca:2021wxt  ($M=2.08\pm
0.07M_{\odot}$) and direct radius measurements of PSR J0740+6620
Miller:2021qha ; Riley:2021pdl  and PSR J0030+0451 Riley:2019yda ;
Miller:2019cac  obtained by the NICER experiment (pink ellipses) as well as
from the measurement of the X-ray binary 4U 1702-429 Nattila:2017wtj  (cyan
area; see Jokela:2021vwy  for a more detailed analysis of the impact of the
NICER results). Finally, shown with a black dashed line is a theoretical lower
bound for the radii as computed from considerations on the threshold mass by
Koeppel et al. Koeppel2019 . We note that we have checked that our EOSs
respect the upper bound on the binary tidal deformability
$\tilde{\Lambda}<720$ (low-spin priors) obtained from GW171817 and are well
within the bounds recently proposed on general parametrization of the sound
speed in neutron stars Altiparmak:2022 .

##  3 Numerical Methods

A crucial component in the numerical simulation of binary neutron star mergers
is accurate initial data for the subsequent time evolution. For this we use
the recently developed Frankfurt University/Kadath (FUKA) (Papenfort2021b, )
code suite, based on an extended version of the KADATH spectral solver library
(Grandclement09, ). FUKA is able to generate realistic compact binary systems
even with large mass asymmetry and sizeable individual spins. For the binary
evolution we make use of the Einstein Toolkit EinsteinToolkit_etal:2020_11
infrastructure that includes the fixed-mesh box-in-box refinement framework
Carpet Schnetter-etal-03b , which keeps track of the location of the
individual stars and the merger remnant in order to provide a fine
computational grid where this is needed. For example, in the simulations
discussed below, six refinement levels with finest grid-spacing of $221\,{\rm
m}$ were used. To solve the fully general relativistic hydrodynamic equations
in 3+1 spacetime dimensions, we use a number of code components developed by
the Frankfurt group and that are implemented in the Einstein Toolkit. For the
spacetime evolution we use Antelope (Most2019b, ), which solves a constraint
damping formulation of the Z4 system (Bernuzzi:2009ex, ; Alic:2011a, ).
Furthermore, to evolve the hydrodynamic part we use the Frankfurt/Illinois
(FIL) general-relativistic magnetohydrodynamic code (Most2019b, ), which is
based on the IllinoisGRMHD code (Etienne2015, ) and implements fourth order
conservative finite-differencing methods (DelZanna2007, ), enabling a precise
hydrodynamic evolution even at relatively low resolution.

Figure 2:  From top to bottom rows we show snapshots of the baryon-number
density, temperature, and quark fraction in the orbital plane at three
different times $t-t_{\rm mer}=1.4,\,3.3,\,5.6\,\,{\rm ms}$ (from left to
right) that are representative for hot, warm and cold quark stages of a
$q=0.7$ simulation with soft V-QCD EOS. In addition, we mark contours for the
quark fraction $Y_{\rm quark}=0.01,\,0.1,\,0.5$ by dotted, dashed and solid
lines, respectively. See Tootle:2022pvd  for details.

##  4 Results

We have performed a set of fully general relativistic hydrodynamics
simulations of binary neutron star systems with the measured chirp mass
$\mathcal{M}_{c}=\frac{(M_{1}\,M_{2})^{3/5}}{(M_{1}+M_{2})^{1/5}}=1.186\,M_{\odot}$
of GW170817 and two different ratios of the individual neutron star masses
$q=\frac{M_{2}}{M_{1}}=0.7,1$. We simulate these binaries using the
aforementioned variants (soft, intermediate, stiff) of the V-QCD model. In
Fig. 3 we show snapshots of various quantities in the orbital plane of the
HMNS at three different times ($1.4-5.6$ ms) after the merger. The case
depicted is the result of an unequal mass binary with $q=0.7$ and soft V-QCD
EOS. The top row shows the baryon number density $n_{b}$ in units of $n_{s}$.
Regions where quark matter is formed are enclosed by dotted, dashed and solid
black lines. The initially highly deformed HMNS becomes more circular and
compact as time proceeds until a dense quark matter core forms, which then
ultimately leads to a phase transition triggered collapse of the merger
remnant into a black hole (not shown). The second row shows the corresponding
temperature distribution. Here one of our main findings becomes visible. Right
after the collision (first column), strong shocks locally lead to temperatures
beyond $50$ MeV inside the HMNS. This is where hot quark matter is formed in
the hottest but not in the densest regions of the HMNS. The second row shows
the situation a few milliseconds later, where a small amount of warm quark
matter is formed at intermediate densities and temperature. In the third row
the situation right before the black hole collapse is shown. Here a massive
core of cold pure quark matter is formed in the central region of the HMNS
where the density is high, but the temperature is low. The decrease in
temperature of the quark matter core can be explained by the sizeable latent
heat that is necessary to complete the transition from baryonic to quark
matter. Finally, in the third row we report the fraction of quark matter
$Y_{\rm quark}$ produced in the early hot, intermediate warm and late cold
quark stages of the HMNS.

In Fig. 3 we show for three representative simulations of $q=0.7$ binaries the
waveform and their frequency spectra that are extrapolated to the $40$ Mpc
luminosity distance of the GW170817 event. First, second and third column are
the plus-polarization of the dominant $\ell=m=2$ mode of the gravitational
wave strain, the power spectral density together with the sensitivities of the
advanced LIGO and the Einstein Telescope detectors, and the gravitational wave
spectrogram, respectively, where coloured dashed lines indicate various
characteristic post-merger peak-frequencies. The first line shows the results
for the soft V-QCD EOS where the quark phase has been excluded, the second
line shows the same configuration where the quark phase is included, and the
third line is the corresponding result for the intermediate V-QCD EOS.
Comparing these results shows that the dominant effect of quark matter is the
early phase transition triggered black hole collapse that leads to an early
termination of the wave form (first versus second line). The post-merger
frequencies in all three cases turn out to be almost identical, showing that
the three quark matter stages in our simulations have negligible impact on the
details of the waveform.

Figure 3: Shown above is the $h^{22}_{+}$ of the strain (left), the power-
spectral density of the effective strain (middle), and the spectrogram (right)
of the $h^{22}_{+}$ of the soft model with (top) and without quark matter
(center) as well as the intermediate model (bottom) for $q=0.7$ and
$\mathcal{M}_{\rm chirp}=1.188$. See Tootle:2022pvd  for details.

##  5 Summary

We carried out fully general relativistic hydrodynamic simulations of neutron
star binary systems that are similar to GW170817. These simulations employ
different variants of holographic V-QCD EOSs allowed us to identify three
different post-merger stages in which hot, warm and cold quark matter is
formed. The dominant impact of quark matter on the HMNS evolution is a phase
transition triggered collapse to a black hole that results in an early
truncation of the gravitational wave signal. The impact of quark matter on the
post-merger kHz gravitational wave spectrum turns out to be small. Future work
focuses on the impact of quark matter on the threshold mass of promptly
collapsing heavy binaries and the long-time stability of lighter binary
systems. The EOS files used in the simulations are publicly available on the
CompOSE database Typel:2013rza ; Typel:2022lcx  under these links
DEJ(DD2-VQCD) soft, DEJ(DD2-VQCD) intermediate, DEJ(DD2-VQCD) stiff.

## Acknowlegments

CE acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) through the CRC-TR 211 ’Strong-interaction matter under
extreme conditions’– project number 315477589 – TRR 211. ST, KT, and LR
acknowledge the support by the State of Hesse within the Research Cluster
ELEMENTS (Project ID 500/10.006). TD and MJ have been supported by an
appointment to the JRG Program at the APCTP through the Science and Technology
Promotion Fund and Lottery Fund of the Korean Government. TD and MJ have also
been supported by the Korean Local Governments – Gyeongsangbuk-do Province and
Pohang City – and by the National Research Foundation of Korea (NRF) funded by
the Korean government (MSIT) (grant number 2021R1A2C1010834). T.D. also
acknowledges the support of the Narodowe Centrum Nauki (NCN) Sonata Bis Grant
No. 2019/34/E/ST3/00405. LR acknowledges funding by the ERC Advanced Grant
“JETSET: Launching, propagation and emission of relativistic jets from binary
mergers and across mass scales” (Grant No. 884631). The simulations were
performed on HPE Apollo HAWK at the High Performance Computing Center
Stuttgart (HLRS) under the grant BNSMIC. For our visualisations we made an
extensive use of the Kuibit library kuibit21 .

## References

  * (1) K. Hebeler, J.M. Lattimer, C.J. Pethick, A. Schwenk, Astrophys. J. 773, 11 (2013), 1303.4662
  * (2) S. Gandolfi, J. Lippuner, A.W. Steiner, I. Tews, X. Du, M. Al-Mamun, Journal of Physics G Nuclear Physics 46, 103001 (2019), 1903.06730
  * (3) J. Keller, C. Wellenhofer, K. Hebeler, A. Schwenk, Phys. Rev. C 103, 055806 (2021), 2011.05855
  * (4) C. Drischler, J.A. Melendez, R.J. Furnstahl, D.R. Phillips, Phys. Rev. C 102, 054315 (2020), 2004.07805
  * (5) E.S. Fraga, A. Kurkela, A. Vuorinen, Astrophys. J. Lett. 781, L25 (2014), 1311.5154
  * (6) T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi, A. Vuorinen, Phys. Rev. D 104, 074015 (2021), 2103.07427
  * (7) T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi, A. Vuorinen, Phys. Rev. Lett. 127, 162003 (2021), 2103.05658
  * (8) H.T. Cromartie et al. (NANOGrav), Nature Astron. 4, 72 (2019), 1904.06759
  * (9) E. Fonseca et al., Astrophys. J. Lett. 915, L12 (2021), 2104.00880
  * (10) T.E. Riley et al., Astrophys. J. Lett. 887, L21 (2019), 1912.05702
  * (11) M.C. Miller et al., Astrophys. J. Lett. 887, L24 (2019), 1912.05705
  * (12) J. Nättilä, M.C. Miller, A.W. Steiner, J.J.E. Kajava, V.F. Suleimanov, J. Poutanen, Astron. Astrophys. 608, A31 (2017), 1709.09120
  * (13) M.C. Miller et al., Astrophys. J. Lett. 918, L28 (2021), 2105.06979
  * (14) T.E. Riley et al., Astrophys. J. Lett. 918, L27 (2021), 2105.06980
  * (15) The LIGO Scientific Collaboration, the Virgo Collaboration, B.P. Abbott, R. Abbott, T.D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso et al. (LIGO Scientific Collaboration and Virgo Collaboration), Physical Review X 9, 011001 (2019), 1805.11579
  * (16) E. Annala, T. Gorda, A. Kurkela, J. Nättilä, A. Vuorinen, Nature Physics 16, 907 (2020)
  * (17) E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis, A. Vuorinen, Physical Review X 12, 011058 (2022), 2105.05132
  * (18) S. Altiparmak, C. Ecker, L. Rezzolla, The Astrophysical Journal Letters 939, L34 (2022)
  * (19) C. Ecker, L. Rezzolla, The Astrophysical Journal Letters 939, L35 (2022)
  * (20) C. Ecker, L. Rezzolla (2022), 2209.08101
  * (21) A. Bauswein, N.U.F. Bastian, D.B. Blaschke, K. Chatziioannou, J.A. Clark, T. Fischer, M. Oertel, Phys. Rev. Lett. 122, 061102 (2019), 1809.01116
  * (22) L.R. Weih, M. Hanauske, L. Rezzolla, Phys. Rev. Lett. 124, 171103 (2020), 1912.09340
  * (23) T. Demircik, C. Ecker, M. Järvinen, Phys. Rev. X 12, 041012 (2022), 2112.12157
  * (24) E.R. Most, L.J. Papenfort, L. Rezzolla, Mon. Not. R. Astron. Soc. 490, 3588 (2019), 1907.10328
  * (25) L.J. Papenfort, S.D. Tootle, P. Grandclément, E.R. Most, L. Rezzolla, Phys. Rev. D 104, 024057 (2021)
  * (26) S. Tootle, C. Ecker, K. Topolski, T. Demircik, M. JÃrvinen, L. Rezzolla, SciPost Phys. 13, 109 (2022)
  * (27) M. Järvinen, Eur. Phys. J. C 82, 282 (2022), 2110.08281
  * (28) J.M. Maldacena, Adv. Theor. Math. Phys. 2, 231 (1998), hep-th/9711200
  * (29) E. Witten, Adv. Theor. Math. Phys. 2, 253 (1998), hep-th/9802150
  * (30) C. Hoyos, D. Rodríguez Fernández, N. Jokela, A. Vuorinen, Phys. Rev. Lett. 117, 032501 (2016), 1603.02943
  * (31) K. Ghoroku, K. Kashiwa, Y. Nakano, M. Tachibana, F. Toyoda, Phys. Rev. D 99, 106011 (2019), 1902.01093
  * (32) K. Bitaghsir Fadafan, J. Cruz Rojas, N. Evans, Phys. Rev. D 101, 126005 (2020), 1911.12705
  * (33) C. Ecker, M. Järvinen, G. Nijs, W. van der Schee, Phys. Rev. D 101, 103006 (2020), 1908.03213
  * (34) T. Demircik, C. Ecker, M. Järvinen, Astrophys. J. Lett. 907, L37 (2021), 2009.10731
  * (35) N. Kovensky, A. Schmitt, JHEP 09, 112 (2020), 2006.13739
  * (36) S. Pinkanjanarod, P. Burikham, Eur. Phys. J. C 81, 705 (2021), 2007.10615
  * (37) L.A.H. Mamani, C.V. Flores, V.T. Zanchin, Phys. Rev. D 102, 066006 (2020), 2006.09401
  * (38) K. Ghoroku, K. Kashiwa, Y. Nakano, M. Tachibana, F. Toyoda (2021), 2107.14450
  * (39) C. Hoyos, N. Jokela, A. Vuorinen (2021), 2112.08422
  * (40) N. Kovensky, A. Poole, A. Schmitt, Phys. Rev. D 105, 034022 (2022), 2111.03374
  * (41) L. Bartolini, S.B. Gudnason, J. Leutgeb, A. Rebhan, Phys. Rev. D 105, 126014 (2022), 2202.12845
  * (42) U. Gursoy, E. Kiritsis, JHEP 02, 032 (2008), 0707.1324
  * (43) U. Gursoy, E. Kiritsis, F. Nitti, JHEP 02, 019 (2008), 0707.1349
  * (44) F. Bigazzi, R. Casero, A.L. Cotrone, E. Kiritsis, A. Paredes, JHEP 10, 012 (2005), hep-th/0505140
  * (45) R. Casero, E. Kiritsis, A. Paredes, Nucl. Phys. B 787, 98 (2007), hep-th/0702155
  * (46) T. Ishii, M. Järvinen, G. Nijs, JHEP 07, 003 (2019), 1903.06169
  * (47) M. Jarvinen, E. Kiritsis, JHEP 03, 002 (2012), 1112.1261
  * (48) A. Akmal, V.R. Pandharipande, D.G. Ravenhall, Phys. Rev. C 58, 1804 (1998), nucl-th/9804027
  * (49) M. Hempel, J. Schaffner-Bielich, Nuclear Physics A 837, 210 (2010), 0911.4073
  * (50) S. Typel, G. Ropke, T. Klahn, D. Blaschke, H.H. Wolter, Phys. Rev. C 81, 015803 (2010), 0908.2344
  * (51) P.A. Zyla et al. (Particle Data Group), PTEP 2020, 083C01 (2020)
  * (52) D.H. Rischke, M.I. Gorenstein, H. Stoecker, W. Greiner, Z. Phys. C 51, 485 (1991)
  * (53) N. Jokela, M. Järvinen, J. Remes, JHEP 03, 041 (2019), 1809.07770
  * (54) N. Jokela, M. Järvinen, G. Nijs, J. Remes, Phys. Rev. D 103, 086004 (2021), 2006.01141
  * (55) N. Jokela, M. Järvinen, J. Remes, Phys. Rev. D 105, 086005 (2022), 2111.12101
  * (56) S. Koeppel, L. Bovard, L. Rezzolla, Astrophys. J. Lett. 872, L16 (2019), 1901.09977
  * (57) P. Grandclement, J. Comput. Phys. 229, 3334 (2010), 0909.1228
  * (58) R. Haas, et al., _The Einstein Toolkit_ (2020), to find out more, visit http://einsteintoolkit.org, https://doi.org/10.5281/zenodo.4298887
  * (59) E. Schnetter, S.H. Hawley, I. Hawke, Class. Quantum Grav. 21, 1465 (2004), gr-qc/0310042
  * (60) S. Bernuzzi, D. Hilditch, Phys. Rev. D 81, 084003 (2010), 0912.2920
  * (61) D. Alic, C. Bona-Casas, C. Bona, L. Rezzolla, C. Palenzuela, Phys. Rev. D 85, 064040 (2012), 1106.2254
  * (62) Z.B. Etienne, V. Paschalidis, R. Haas, P. Mösta, S.L. Shapiro, Class. Quantum Grav. 32, 175009 (2015), 1501.07276
  * (63) L. Del Zanna, O. Zanotti, N. Bucciantini, P. Londrillo, Astron. Astrophys. 473, 11 (2007), 0704.3206
  * (64) S. Typel, M. Oertel, T. Klähn, Phys. Part. Nucl. 46, 633 (2015), 1307.5715
  * (65) S. Typel et al. (2022), 2203.03209
  * (66) G. Bozzola, The Journal of Open Source Software 6, 3099 (2021), 2104.06376


# An optimization-based equilibrium measure describes non-equilibrium steady
state dynamics: application to edge of chaos

Junbin Qiu1 Haiping Huang1,2<EMAIL_ADDRESS>1PMI Lab, School of
Physics, Sun Yat-sen University, Guangzhou 510275, People’s Republic of China
2Guangdong Provincial Key Laboratory of Magnetoelectric Physics and Devices,
Sun Yat-sen University, Guangzhou 510275, People’s Republic of China

###### Abstract

Understanding neural dynamics is a central topic in machine learning, non-
linear physics and neuroscience. However, the dynamics is non-linear,
stochastic and particularly non-gradient, i.e., the driving force can not be
written as gradient of a potential. These features make analytic studies very
challenging. The common tool is to use path integral approach or dynamical
mean-field theory, but the drawback is one has to solve the integro-
differential or dynamical mean-field equations, which is computationally
expensive and has no closed form solutions in general. From the aspect of
associated Fokker-Planck equation, the steady state solution is generally
unknown. Here, we treat searching for the steady state as an optimization
problem, and construct an approximate potential closely related to the speed
of the dynamics, and find that searching for the ground state of this
potential is equivalent to running a stochastic gradient dynamics. The
resultant stationary state follows exactly the canonical Boltzmann measure.
Within this framework, the quenched disorder intrinsic in the neural networks
can be averaged out by applying the replica method. Our theory reproduces the
well-known result of edge-of-chaos, and further the order parameters
characterizing the continuous transition are derived, and different scaling
behavior with respect to inverse temperature in both sides of the transition
is also revealed. Our method opens the door to analytically study the steady
state landscape of the deterministic or stochastic high dimensional dynamics.

##  I Introduction

Non-equilibrium dynamics is a central topic of statistical physics [1, 2, 3],
from which general fluctuation-dissipation relations can be established [4, 5,
6], and path integral approach was developed to derive the relationship
between correlation and response function of the dynamics [7, 8, 9]. One
benchmark for testing or deriving theoretical hypotheses is the Langevin
dynamics with conservative (expressed as gradient of a potential) or
nonconservative forces, including Ising spin glass dynamics [10], steady state
thermodynamics [11], and neural dynamics with partially symmetric couplings
[12]. These studies always resort to dynamical mean-field theory (DMFT) [13,
14], where all possible trajectory paths are integrated to give a physical
action, from which correlation and response functions are derived as an
optimization of the action [15]. In a general context, the dynamics would be
very complicated, e.g., non-gradient and highly coupled, which makes the DMFT
very complicated, and even it remains unknown how to capture the steady state
properties of these dynamics.

It is commonly observed that the speed of dynamics determines whether the
system enters a steady state and thus can be identified [16]. In this paper,
we propose a new angle that only the low-speed or zero speed region of the
phase space is considered. In this sense, one can write down an energy cost
function to optimize for satisfying this speed constraint. This optimization
can be formulated under the stochastic gradient dynamics, where a temperature
like parameter is introduced. This temperature plays the same role as in the
canonical Boltzmann measure. In fact, in the zero temperature and stationary
limit, the gradient dynamics are self-consistent with the original general
dynamics. Therefore, if the ground state of the energy cost is concentrated
on, one can capture the steady state properties of the original phase space.
This has an additional benefit that powerful equilibrium concepts such as
Boltzmann measure, phase transitions, and order parameters can be leveraged
[17].

We verify our framework in recurrent neural networks with asymmetric
couplings, where the deterministic dynamics yields order-to-chaos transition
[18, 19]. We are able to handle the quenched disorder due to the coupling
statistics with the help of replica method, a powerful trick in spin glass
theory [20]. Furthermore, two types of order parameters are derived; one
describes the network activity, while the other describes the response
property of the dynamic system. The former order parameter exhibits a
continuous phase transition from null to positive value as the network
parameters are tuned, while the latter one shows a peak at the exact location
of transition, thereby offering rich physical information about the steady
state of the complicated non-gradient dynamics. Different scaling behavior
with increasing inverse temperature was found on both sides of the transition
point. Our framework also opens the door to analytically study the quasi-
potential landscape of a general dynamics, using large-deviation principle and
landscape analysis in equilibrium statistical mechanics [21, 17], which is
left for a future work.

Our approximate potential analysis is not restricted to simple type of
dynamics, and can be generalized to treat competition dynamics of many species
in ecosystems [22, 23, 24], high-dimensional dynamics of random optimization
[25, 26], high-dimensional gradient flow in machine learning [27], and
adaptive, hierarchical and noisy neural dynamics in neuroscience [28, 29, 30,
31, 32], as long as the dynamics speed can be identified and the steady state
exists. The only limitation is that transient behavior before the steady state
is reached could not be captured by our equilibrium approximation, and in that
case, integro-differential equations are required to solve with expensive
computational cost [15]. We next show how physical insights about a
complicated non-gradient dynamics can be derived using our approximate
potential analysis.

##  II Model

We consider a recurrent neural network (RNN) composed of $N$ interacting
neurons, whose architecture is shown in Fig. 1 (a), illustrated as a simple
example of $N=3$. The network activity is described by $N$ real-valued
synaptic currents ($\mathbf{x}$). The connection strength $J_{ij}$ measures
the asymmetric impact from neuron $j$ to neuron $i$. Note that in our model,
two neurons do not have to interact with each other equally, i.e. $J_{ij}\neq
J_{ji}$. In addition, $J_{ii}=0$. All neurons are connected to each other. The
synaptic currents are transferred to firing rates through the current-to-rate
transfer function $\phi(\cdot)$. Each neuron integrates neighboring activity
via the incoming asymmetric synapses, and therefore the dynamics of the
current is expressed as the following first-order non-linear differential
equation for all $i$ [33, 18]

$\derivative{x_{i}}{t}=-x_{i}+\sum_{j=1}^{N}J_{ij}\phi(x_{j}),$ (1)

where we consider unit time scale. Unless otherwise stated, $\phi=\tanh$. In
this case, the transferred neural activity may become negative, while for
biological reality, the non-negative function form should be chosen. But our
following analysis is still applicable to this type of non-negative transfer
function.

In real neural circuits, $J_{ij}$ may not be identically independently
distributed (i.i.d.). We therefore consider the correlated case whose
statistics is specified below

$\displaystyle\left\langle\quantity(J_{ij})^{2}\right\rangle=$
$\displaystyle\frac{g^{2}}{N},$ (2) $\displaystyle\left\langle
J_{ij}J_{ji}\right\rangle=$ $\displaystyle\frac{g^{2}}{N}\gamma,$ (3)

where $g>0$ is called the gain parameter, which controls the variance of
connections, and $\gamma\in[-1,1]$ controls the degree of connection symmetry.
The special cases of $\gamma=1$ and $\gamma=-1$ correspond to symmetric and
anti-symmetric neural networks, respectively, while $\gamma=0$ corresponds to
fully asymmetric networks [18].

The above dynamics model exhibits an order-to-chaos transition [19]. The
typical number of fixed points explodes as the synaptic gain exceeds a
critical value [34], which can be derived from the rightmost boundary of the
asymptotic eigenvalue spectrum of $\mathbf{J}$ [35]. When $g(1+\gamma)<1$, the
network stays in the trivial fixed-point phase, where all trajectories flow to
zero activity stable state. Once $g(1+\gamma)>1$, the network activity enters
a chaotic phase, where an infinitesimal separation of two trajectories would
be eventually amplified, and individual neurons exhibit slow to fast time
scales as the network parameters go deep into the chaotic phase. We illustrate
the phase diagram in Fig. 1 (b) and show the representative dynamic behaviors
of these two phases in Fig. 1 (c).

Intuitively, the phase space is partitioned into different regions of
different speeds ($\|\derivative{\mathbf{x}}{t}\|^{2}$). We are interested in
the zero-speed or low-speed regions [16]. Focusing on the speed landscape, we
design the following approximate potential (quasi-potential) of the dynamics,

$E(\mathbf{x})=\frac{1}{2}\sum_{i}\quantity(-x_{i}+\sum_{j=1}^{N}J_{ij}\phi(x_{j}))^{2}+\eta\norm*{\mathbf{x}}^{2},$
(4)

where the first term is the kinetic energy, while the second term is a
regularization term constraining the $\ell_{2}$ norm of the network activity.
$\eta$ plays the role of tuning the relative strength with respect to the
kinetic energy.

Figure 1: Key idea of our quasi-potential method and phase diagram of the
considered RNN model. (a) A three-neuron recurrent neural network is shown in
the left panel. Neurons are bidirectionally connected with the strength
displayed near the arrows (synaptic directions). Both connections in opposite
directions are correlated once $\gamma\neq 0$. Here, we simulate the RNN
dynamics [Eq. (1)], whose phase space displays two fixed points
($g=1.2,\gamma=0.1$). In the right panel, we show the trajectories flowing to
these fixed points, where the color and size of the dots represent the energy
of the corresponding position in the phase space, and the arrows represent the
direction and the amount of the velocity (the length of the arrow). We also
show an example of the trivial fixed-point for comparison
($g=0.8,\gamma=0.1$), where nearby trajectories are shown in the light grey
color. The cross symbols denote these corresponding fixed points. (b) The
phase diagram as a function of gain parameter $g$ and symmetric degree
$\gamma$ is divided into two regions: null fixed-point and chaotic phases,
separated by an analytic transition line $g(1+\gamma)=1$ [19]. (c) We show the
representative behavior of the dynamics [$\phi(x)$] in the two phases,
corresponding to two colored points in the phase diagram of (b).

According to Eq. (4), we can write down the following stochastic gradient
dynamics (or Langevin dynamics)

$\derivative{\mathbf{x}}{t}=-\nabla_{\mathbf{x}}E(\mathbf{x})+\sqrt{2T}\bm{\epsilon},$
(5)

where $\bm{\epsilon}(t)$ is a time-dependent white noise, whose statistics is
given by $\left\langle\epsilon_{i}(t)\right\rangle=0$,
$\left\langle\epsilon_{i}(t)\epsilon_{j}(t^{\prime})\right\rangle=\delta_{ij}\delta(t-t^{\prime})$.
The temperature $T$ is used to tune the energy level, playing the same role as
in the equilibrium Boltzmann measure. We can write the gradient (force) in Eq.
(5) in the component wise,

$F_{i}\equiv-x_{i}+h_{i}-\phi^{\prime}(x_{i})\sum_{j:j\neq
i}J_{ji}(h_{j}-x_{j}),$ (6)

where $h_{i}\equiv\sum_{j:j\neq i}J_{ij}\phi(x_{j})$. From this expansion, we
can see that only in the stationary limit ($T\to 0$ as well), the above
Langevin dynamics closely resembles the original dynamics (the same magnitude
order of small speeds). In this sense, the energy defined in Eq. (4) is not a
Lyapunov function whose analytic form for the current out-of-equilibrium
dynamics with asymmetric couplings is unknown, although the driving force is
argued to be decomposed into gradient of a potential related to the steady
state distribution and a divergent-free curl flux [36]. The decomposition
takes into account the probability conservation of neural state evolution. In
contrast, we start from an optimization angle that concentrates on the low-
speed region of the phase space, and formulate the steady state behavior as
the canonical ensemble of stationary fixed points.

It is well-known that the steady state of the above Langevin dynamics has an
equilibrium Boltzmann measure [1],

$P(\mathbf{x})=\frac{1}{Z}e^{-\beta E(\mathbf{x})},$ (7)

where the partition function $Z$ depends on particular realizations of
$\mathbf{J}$, and $\beta=1/T$. An illustration of $E(\mathbf{x})$ is given in
Fig. 1 (a). Our main focus is to compute the free energy of this canonical
ensemble, and derive the steady state properties of the dynamics model. The
theoretical predictions can be confirmed by simulating a discretized Langevin
dynamics, assuming the small discretization step as learning rate, akin to a
Monte-Carlo simulation. We remark that this equilibrium analysis reproduces
the continuous nature of the dynamics transition, more insights and
possibility of exploring the speed landscape would be discussed in the
remaining sections.

##  III Replica method

To derive the quenched disorder average of the free energy, we apply the
replica method, which is popular in studying theory of spin glasses [20] and
statistical mechanics of neural networks [17]. We first write down the
following replicated partition function

$Z^{n}=\int
d\mathbf{x}e^{-\beta\pqty{\frac{1}{2}\sum_{ia}\pqty{-x_{i}^{a}+\sum_{j}J_{ij}\phi(x_{j}^{a})}^{2}+\eta\sum_{a}||\mathbf{x}^{a}||^{2}}},$
(8)

where $i,a$ denote the neuron and replica index, respectively. We denote
$d\mathbf{x}$ as a shorthand for $\prod_{i=1}^{N}\prod_{a=1}^{n}dx_{i}^{a}$,
where $x_{i}^{a}\in\mathbb{R}$. Therefore, $-\beta f=\lim_{n\to
0}\frac{\ln\left\langle Z^{n}\right\rangle}{Nn}$, where
$\left\langle\cdot\right\rangle$ denotes the disorder average over
$\mathbf{J}$. The neural dynamics is much faster than the coupling dynamics,
which we assume a stationary limit of fixed $\mathbf{J}$, and thus the
quenched disorder average is required. But if both dynamics evolve in similar
time scales, an enlarged state space may be considered, and we leave this
issue to future works.

In the course of replica analysis (see details in Appendix A), two order
parameters are naturally introduced as follows,

$\displaystyle Q^{ab}=$
$\displaystyle\frac{1}{N}\sum_{i}\phi(x_{i}^{a})\phi(x_{i}^{b}),$ (9)
$\displaystyle R^{ab}=$
$\displaystyle\frac{1}{N}\sum_{i}\hat{x}_{i}^{a}\phi(x_{i}^{b}),$ (10)

where $\hat{x}_{i}^{a}$ acts as a response field, whose name is inspired by
the analysis using DMFT [15] and explained in the next section. As the first
level of approximation, we assume the replica symmetric ansätz to order
parameters $Q^{ab}=q\delta_{ab}+Q(1-\delta_{ab})$ and
$R^{ab}=r\delta_{ab}+R(1-\delta_{ab})$, whose physical meaning is explained in
more details in Appendix A. Under this approximation, the free energy is given
by

$\displaystyle-\beta f=$
$\displaystyle\frac{1}{2}Q\hat{Q}-q\hat{q}+R\hat{R}-r\hat{r}-\ln\sigma+\frac{1}{2}\beta
g^{2}\gamma\quantity(r^{2}-R^{2})$ (11) $\displaystyle+\int\pqty{DuDv}\ln I,$

where $\sigma\equiv\sqrt{1+g^{2}\beta(q-Q)}$, and $\hat{Q}$, $\hat{q}$,
$\hat{R}$ and $\hat{r}$ are conjugated order parameters. The integral
$I\equiv\int dxe^{\mathcal{H}(x)}$, where

$\displaystyle\mathcal{H}(x)\equiv$ $\displaystyle-\beta\eta
x^{2}+\frac{1}{2}\quantity(2\hat{q}-\hat{Q})\phi^{2}(x)$ (12)
$\displaystyle+\quantity(\sqrt{\frac{g^{2}\beta
Q\hat{Q}-\hat{R}^{2}}{g^{2}\beta Q}}u+\frac{\hat{R}}{g\sqrt{\beta
Q}}v)\phi(x)$ $\displaystyle-\frac{1}{2\sigma^{2}}\quantity(g\sqrt{\beta
Q}v+(\hat{r}-\hat{R})\phi(x)-\sqrt{\beta}x)^{2}.$

In the large $N$ limit, we can use the saddle-point approximation (also called
Laplace method) to derive the closed equations the order parameters must obey.
Therefore, by setting the derivatives of the replica free energy with respect
to the order parameters zero, we obtain the following compact self-consistent
equations (see Appendix B for details)

$\displaystyle q=$
$\displaystyle\quantity[\left\langle\phi^{2}\right\rangle],$ (13a)
$\displaystyle Q=$
$\displaystyle\quantity[\left\langle\phi\right\rangle^{2}],$ (13b)
$\displaystyle\hat{q}=$
$\displaystyle-\frac{gk}{2}+\frac{g^{2}k^{2}Q}{2}+\frac{gk}{2\sigma^{2}}\quantity(\hat{r}-\hat{R})f(1,1,-2)\quantity[\left\langle\phi^{2}\right\rangle]$
(13c) $\displaystyle+$
$\displaystyle\frac{gk}{\sigma^{2}}\quantity(\hat{r}-\hat{R})f(0,-1,1)\quantity[\left\langle\phi\right\rangle^{2}]+\frac{k^{2}}{2}\quantity(1-2gkQ)\quantity[\left\langle
x^{2}\right\rangle]$ (13d) $\displaystyle+$
$\displaystyle\frac{gk\sqrt{\beta}}{\sigma^{2}}f(0,1,-2)\quantity[\left\langle
x\right\rangle\left\langle\phi\right\rangle]+\frac{gk\sqrt{\beta}}{\sigma^{2}}f(-1,0,2)\quantity[\left\langle
x\phi\right\rangle]$ (13e) $\displaystyle+$ $\displaystyle
gk^{3}Q\quantity[\left\langle x\right\rangle^{2}],$ (13f)
$\displaystyle\hat{Q}=$ $\displaystyle
g^{2}k^{2}Q+\frac{2gk}{\sigma^{2}}\quantity(\hat{r}-\hat{R})f(0,1,-1)\quantity[\left\langle\phi^{2}\right\rangle]$
(13g) $\displaystyle+$
$\displaystyle\frac{gk}{\sigma^{2}}\quantity(\hat{r}-\hat{R})f(1,-3,2)\quantity[\left\langle\phi\right\rangle^{2}]-2gk^{3}Q\quantity[\left\langle
x^{2}\right\rangle]$ (13h) $\displaystyle-$
$\displaystyle\frac{2gk\sqrt{\beta}}{\sigma^{2}}f(1,-2,2)\quantity[\left\langle
x\right\rangle\left\langle\phi\right\rangle]+\frac{2gk\sqrt{\beta}}{\sigma^{2}}f(0,-1,2)\quantity[\left\langle
x\phi\right\rangle]$ (13i) $\displaystyle+$ $\displaystyle
k^{2}\quantity(1+2gkQ)\quantity[\left\langle x\right\rangle^{2}],$ (13j)
$\displaystyle r=$
$\displaystyle-\frac{1}{\sigma^{2}}f(1,0,-1)\quantity[\left\langle\phi^{2}\right\rangle]+\frac{1}{\sigma^{2}}f(0,1,-1)\quantity[\left\langle\phi\right\rangle^{2}]$
(13k) $\displaystyle+$
$\displaystyle\frac{\sqrt{\beta}}{\sigma^{2}}\quantity(1-gkQ)\quantity[\left\langle\phi
x\right\rangle]+\frac{\sqrt{\beta}}{\sigma^{2}}gkQ\quantity[\left\langle\phi\right\rangle\left\langle
x\right\rangle],$ (13l) $\displaystyle R=$
$\displaystyle-\frac{1}{\sigma^{2}}f(0,1,-1)\quantity[\left\langle\phi^{2}\right\rangle]-\frac{1}{\sigma^{2}}f(1,-2,1)\quantity[\left\langle\phi\right\rangle^{2}]$
(13m) $\displaystyle-$
$\displaystyle\frac{\sqrt{\beta}}{\sigma^{2}}gkQ\quantity[\left\langle
x\phi\right\rangle]+\frac{\sqrt{\beta}}{\sigma^{2}}\quantity(1+gkQ)\quantity[\left\langle
x\right\rangle\left\langle\phi\right\rangle],$ (13n) $\displaystyle\hat{r}=$
$\displaystyle\beta g^{2}\gamma r,$ (13o) $\displaystyle\hat{R}=$
$\displaystyle\beta g^{2}\gamma R,$ (13p)

where $k\equiv\frac{g\beta}{\sigma^{2}}$ and $f(a,b,c)\equiv
a\hat{r}+b\hat{R}+cgkQ\quantity(\hat{r}-\hat{R})$. In addition, two kinds of
average are specified. One is the disorder average $\quantity[\cdot]\equiv\int
DuDv\cdot$, and the other is the thermal average
$\left\langle\cdot\right\rangle\equiv\frac{\int
dxe^{\mathcal{H}(x)}\cdot}{\int dxe^{\mathcal{H}(x)}}$. $\mathcal{H}(x)$ plays
the role of effective Hamiltonian of our model. When $\gamma=0$ (i.i.d.
scenario), $\hat{r},\hat{R}$ and $f(a,b,c)$ vanish, the above saddle point
equation can be greatly simplified (see Appendix C).

##  IV Results and discussion

###  IV.1 Equilibrium replica analysis reveals the continuous nature of the
chaos transition

The zero temperature limit selects the ground state of the Boltzmann measure
[Eq.(7)], which captures the statistics of the steady state solutions of the
original non-gradient dynamics [Eq. (1)]. An explicit form of Lyapunov
function that is monotonically is decreasing over time can only be identified
provided that the neural coupling is symmetric, and only in this special case,
the Lyapunov function reduces to an Ising-model like Hamiltonian, as
encountered in standard Hopfield networks [37]. We thus remark that the quasi-
potential defined in Eq. (4) is not the physical potential of the original
dynamics, yet capturing only the steady state properties of the dynamics,
e.g., the low- (or zero-)speed landscape. In this sense, not all
characteristics of a non-gradient dynamics are included in our approximate
potential description, e.g., the transient behavior or the temporal evolution
of autocorrelation even in the time translation invariant regime is not
captured. However, the statistics of fixed-point solutions can be well
captured by the quasi-potential approach.

Figure 2: The average energy (a) and $\ell_{2}$ norm (b) of network activity
versus the inverse temperature $\beta$ in the i.i.d. scenario $\gamma=0$ with
$\eta=0$. The orange and the blue curves represent the setting $g=0.8$ and
$g=1.2$, respectively. The red dashed line denotes the ground state energy.
Two lines of $g=0.8$ and $g=1.2$ are nearly indistinguishable.

The average energy is an indicator of whether the RNN dynamics is trapped in
fixed points. In other words, steady states are achievable provided that the
average energy can be optimized to the ground state zero value. Without loss
of generality, we consider the independent setting $\gamma=0$ and set
$\eta=0$. As shown in Fig. 2, with increasing value of $\beta$, the energy
decreases to $0$, implying that the ground state can be achieved in a very
large $\beta$. When $\beta=10000$, the energy becomes less than $10^{-4}$,
suggesting that it is reasonable to set $\beta=10000$ for the zero temperature
limit considered in our numerical solutions of the SDEs and experiments
verifying our theoretical results. We also show the $\ell_{2}$ norm of the
network activity in Fig. 2. The norm displays different magnitudes before and
after the onset of chaos.

With a finite temperature, the low-speed solutions of the original dynamics
could also be studied. For the remaining analysis, we focus on a very large
$\beta$ (e.g., $\beta=10^{4}$) and set $\eta=0$ unless otherwise stated.
Varying the control parameters $g$ and $\gamma$, we plot the profiles of all
order parameters ($q,Q,r,R$) in Fig. 3.

Our replica analysis reveals exactly the order-to-chaos transition [19]. When
$g<\frac{1}{1+\gamma}$, the zero fixed-point uniquely dominates the phase
space, resulting in vanishing $q$ and $Q$, as expected from the physical
definition of the order parameters. Once $g>\frac{1}{1+\gamma}$, $q$ and $Q$
increases rapidly but continually from zero, indicating the emergence of self-
sustained collective activity. This coincides with the known picture of
exponentially growing topological complexity [34] and the positive maximal
Lyapunov exponent [18, 30, 17]. It would be interesting to count the
stationary points of our potential, rather than using the hard-to-compute Kac-
Rice formula [34], which we left for a future thorough work. Our method thus
provides a relatively simple way to access the dynamics landscape in the
steady state limit. We remark that the specific form of steady state
distribution of the network activity for the original dynamics is unknown.

Let us now look at the associated order parameters involving the response
field $\hat{x}$, which emerges due to the linearization of the quadratic terms
in the Hamiltonian (see details in Appendix A). This response field is similar
to that used to enforce the dynamics equation in a field-theoretic formula of
the stochastic dynamics [15]. But in our current setting, the response field
is time-independent, as only steady state properties are focused on in this
work. Therefore, we argue that $R$ and $r$ bear the same physical meaning with
the response function (equal time) in the DMFT language [15]. In fact, $r$ and
$R$ display a peak at the transition point. This peak is an indicator of
continuous dynamics transition, and may be related to the divergence of the
relaxation time-scale of the auto-correlation function [18, 19]. In other
words, the dynamics at the edge of chaos is more responsive than at other
locations in the phase diagram. The edge of chaos has thus great computational
benefits [38, 39, 17]. Surprisingly, this is also consistent with the
experimental data of cortical electrodynamics [40], which claimed the
information richness at the edge of chaos, and the conscious brain states are
observed at the edge of chaos. We conclude that this peak of the response
function is a meaningful byproduct of our theory, and moreover an indicator of
the edge of chaos.

###  IV.2 Different scaling behavior of order parameters around the transition

In this section, we investigate how the order parameters behave in the zero
temperature limit in a special case of $\gamma=0$, where the transition point
locates at $g=1$ [18]. The conclusion holds for other values of $\gamma$. As
shown in Fig. 4, before the onset of chaos ($g=0.8$), $q$ and $Q$ are
inversely proportional to $\beta$, as $\beta$ increases. In addition,
$\hat{q}$ and $\hat{Q}$ are proportional to $\beta$ (either decreasing or
increasing). In contrast, after the chaos transition ($g=1.2$), the
large-$\beta$ asymptotic scaling behavior changes to a different form
deviating strongly from that in the case of trivial (ordered) phase. Note that
the average energy at $\beta=10^{4}$ already converges to zero (the ground
state value) (see Fig. 2).

Figure 3: Order parameter profile as a function of gain parameters $g$ and
symmetric degree $\gamma$. $g$ is varied from $0.4$ to $1.2$, and $\gamma$ is
varied from $-0.2$ to $1.0$ with an interval of $0.1$. The colors of curves
correspond to the different $\gamma$ values. The grey vertical plane specifies
the location of the critical line $g\quantity(\gamma+1)=1$ [19]. Figure 4:
Order parameters as a function of increasing $\beta$. $q$, $Q$, $\hat{q}$ and
$\hat{Q}$ are plotted in (a), (b), (c) and (d), respectively. We consider
$g=0.8$ (blue dots) in the main figure and $g=1.2$ (orange dots) in the inset,
and tune $\beta$ from $100$ to $10000$. The experimental SDE solutions are
fitted in an asymptotic form shown in the plots.

##  V Concluding remarks

In this paper, we provide an relatively simple way to study non-equilibrium
dynamics, by proposing a quasi-potential. Under this potential, we derive a
stochastic dynamics to implement the optimization of searching for the low-
speed region of the phase space. In the zero temperature limit, only the
steady states are considered, and thus our canonical ensemble theory yields a
good description of the original non-potential dynamics in the long time
limit. More precisely, the low (zero)-speed landscape of the phase space can
be analytically studied, which is relevant to our mechanistic understanding of
recurrent neural dynamics, which can be adaptive, hierarchical and noisy. This
equilibrium transformation allows for using the Boltzmann measure to describe
the low-speed landscape, as could be derived from the Fokker-Planck equation
of the corresponding stochastic dynamics of the optimization [1]. Because of
the quenched disorder present in our model, the replica method can be readily
adapted as well.

The new finding is that two physically relevant order parameters are derived
from our canonical ensemble calculation. One is activity related, which
exactly reproduces the order-to-chaos transition phenomenon, previously
derived by using more complicated DMFT and random matrix theory [18, 35, 19,
15]. The other is response related, which shows a peak at the very transition
point, demonstrating the computation benefits of the edge of chaos. The
continuous nature of this chaos transition is confirmed by our theory as well.
This transition is related to supersymmetry breaking in a field-theoretic
representation of the dynamics in an action [41]. Remarkably, we observe that
as the temperature goes to zero, in both sides of the transition point, the
order parameters show different scaling behavior with increasing inverse
temperature. This is an evidence of dynamical complexity in RNNs.

Our work reveals a quasi-potential landscape behind non-gradient Langevin
dynamics (or deterministic dynamics). This could inspire several important
future directions. First, a complete picture of this landscape can be studied,
as powerful equilibrium approaches can be leveraged, such as large deviation
analysis [21], and landscape geometry analysis [42]. Second, one can ask what
is the quasi-potential of a coupled-dynamics system where the dynamics of
constituent elements can be divided into several levels of different time
scales. An enlarged state space is required to introduce. Lastly, one could
modify the quasi-potential dynamics to approximate the original non-gradient
dynamics with increasing precision.

To sum up, our approach does not rely on the special form of the non-linear
dynamics, as long as the speed of the dynamics can be identified, which
already includes a broad range of neural dynamics models [28, 17]. Even if the
dynamics contains slowly evolving couplings, we can also treat this as a
quenched disorder that can be averaged resulting in a self-averaged free
energy, from which transitions can be detected and analyzed. As some types of
the dynamics can not be accessed by DMFT, or the resultant integro-
differential equations for correlation and response functions are hard to
solve, our quasi-potential method is a promising starting point to understand
the steady state landscape of general non-equilibrium dynamics [43], which
includes competition dynamics of many species in ecosystems [22, 23, 24],
high-dimensional dynamics of random optimization [25, 26], high-dimensional
gradient flow in machine learning [27], and adaptive, hierarchical and noisy
neural dynamics in neuroscience [28, 29, 30, 31, 32].

##  Appendix A Details of replica calculation

Starting from Eq. (8), we perform the Hubbard–Stratonovich (HS) transformation
$e^{-\frac{1}{2}ab^{2}}=\int D\hat{x}e^{\mathrm{i}\sqrt{a}b\hat{x}}$ for
$a>0$, to linearize the quadratic terms and have

$\displaystyle Z^{n}=$ $\displaystyle\int
d\mathbf{x}D\mathbf{\hat{x}}\exp\bqty{\mathrm{i}\sqrt{\beta}\sum_{ia}\hat{x}_{i}^{a}\pqty{-x_{i}^{a}+\sum_{j}J_{ij}\phi(x_{j}^{a})}-\beta\eta\sum_{a}\norm*{\mathbf{x}^{a}}^{2}},$
(14)

where we define $D\hat{x}\equiv
e^{-\frac{1}{2}\hat{x}^{2}}d\hat{x}/\sqrt{2\pi}$ as the Gaussian measure, and
we call $\hat{x}$ as a response field [15]. We define $i,j$ as the neuron
index and $a,b$ as the replica index. The neuron index runs from $1$ to $N$,
while the replica index runs from $1$ to $n$. In addition,
$d\mathbf{x}\equiv\prod_{i,a}dx_{i}^{a}$ and
$D\mathbf{\hat{x}}\equiv\prod_{i,a}D\hat{x}_{i}^{a}$.

Before implementation of the quenched disorder average over the coupling
statistics, it is convenient to separate the connection matrix into a
symmetric part $\mathbf{J}^{S}$ and an anti-symmetric part $\mathbf{J}^{A}$
[12],

$\mathbf{J}=\mathbf{J}^{S}+\rho\mathbf{J}^{A},$ (15)

where by definition $J^{S}_{ij}=J^{S}_{ji}$ and $J^{A}_{ij}=-J^{A}_{ji}$, and
their variances are defined as follows,

$\left\langle\quantity(J^{S}_{ij})^{2}\right\rangle=\left\langle\quantity(J^{A}_{ij})^{2}\right\rangle=\frac{g^{2}}{N}\frac{1}{1+\rho^{2}}.$
(16)

Both parts are uncorrelated. It then follows that

$\left\langle
J_{ij}J_{ji}\right\rangle=\frac{g^{2}}{N}\frac{1-\rho^{2}}{1+\rho^{2}}=\frac{g^{2}}{N}\gamma,$
(17)

where the relationship between $\rho$ and $\gamma$ is given below,

$\rho=\sqrt{\frac{1-\gamma}{1+\gamma}}.$ (18)

Defining $\left\langle\cdot\right\rangle$ as the quenched disorder average, we
have

$\left\langle Z^{n}\right\rangle=\int
d\mathbf{x}D\mathbf{\hat{x}}\exp\bqty{-\mathrm{i}\sqrt{\beta}\sum_{ia}x_{i}^{a}\hat{x}_{i}^{a}-\beta\eta\sum_{a}\norm*{\mathbf{x}^{a}}^{2}}\left\langle\exp\bqty{\mathrm{i}\sqrt{\beta}\sum_{ia}\hat{x}_{i}^{a}\sum_{j}J_{ij}\phi(x_{j}^{a})}\right\rangle.$
(19)

The quenched disorder average is first completed as follows,

$\displaystyle\left\langle\exp\quantity[\mathrm{i}\sqrt{\beta}\sum_{ij}J_{ij}\sum_{a}\hat{x}_{i}^{a}\phi(x_{j}^{a})]\right\rangle$
(20) $\displaystyle=$
$\displaystyle\left\langle\exp\quantity[\mathrm{i}\sqrt{\beta}\sum_{ij}\quantity(J^{s}_{ij}+\rho
J^{a}_{ij})\sum_{a}\hat{x}_{i}^{a}\phi(x_{j}^{a})]\right\rangle$
$\displaystyle=$
$\displaystyle\exp\quantity[-\frac{1}{2}\beta\frac{g^{2}}{N}\frac{1}{1+\rho^{2}}\sum_{i<j}\sum_{ab}\quantity(\hat{x}_{i}^{a}\hat{x}_{i}^{b}\phi(x_{j}^{a})\phi(x_{j}^{b})+\hat{x}_{i}^{a}\hat{x}_{j}^{b}\phi(x_{j}^{a})\phi(x_{i}^{b})+\hat{x}_{j}^{a}\hat{x}_{i}^{b}\phi(x_{i}^{a})\phi(x_{j}^{b})+\hat{x}_{j}^{a}\hat{x}_{j}^{b}\phi(x_{i}^{a})\phi(x_{i}^{b}))]$
$\displaystyle\times\exp\quantity[-\frac{1}{2}\beta\frac{g^{2}}{N}\frac{\rho^{2}}{1+\rho^{2}}\sum_{i<j}\sum_{ab}\quantity(\hat{x}_{i}^{a}\hat{x}_{i}^{b}\phi(x_{j}^{a})\phi(x_{j}^{b})-\hat{x}_{i}^{a}\hat{x}_{j}^{b}\phi(x_{j}^{a})\phi(x_{i}^{b})-\hat{x}_{j}^{a}\hat{x}_{i}^{b}\phi(x_{i}^{a})\phi(x_{j}^{b})+\hat{x}_{j}^{a}\hat{x}_{j}^{b}\phi(x_{i}^{a})\phi(x_{i}^{b}))]$
$\displaystyle=$
$\displaystyle\exp\quantity[-\frac{1}{2}\beta\frac{g^{2}}{N}\sum_{i\not=j}\sum_{ab}\quantity(\frac{1}{1+\rho^{2}}\quantity(\hat{x}_{i}^{a}\hat{x}_{i}^{b}\phi(x_{j}^{a})\phi(x_{j}^{b})+\hat{x}_{i}^{a}\hat{x}_{j}^{b}\phi(x_{j}^{a})\phi(x_{i}^{b}))+\frac{\rho^{2}}{1+\rho^{2}}\quantity(\hat{x}_{i}^{a}\hat{x}_{i}^{b}\phi(x_{j}^{a})\phi(x_{j}^{b})-\hat{x}_{i}^{a}\hat{x}_{j}^{b}\phi(x_{j}^{a})\phi(x_{i}^{b})))]$
$\displaystyle=$
$\displaystyle\exp\quantity[-\frac{1}{2}\beta\frac{g^{2}}{N}\sum_{ij}\sum_{ab}\quantity(\hat{x}_{i}^{a}\hat{x}_{i}^{b}\phi(x_{j}^{a})\phi(x_{j}^{b})+\gamma\hat{x}_{i}^{a}\hat{x}_{j}^{b}\phi(x_{j}^{a})\phi(x_{i}^{b}))]$
$\displaystyle=$ $\displaystyle\exp\quantity[-\frac{1}{2}\beta
g^{2}\sum_{ab}Q^{ab}\sum_{i}\hat{x}_{i}^{a}\hat{x}_{i}^{b}-\frac{1}{2}\beta
g^{2}N\gamma\sum_{ab}R^{ab}R^{ba}],$

where we are left with only $\gamma$ (whose value can now takes from $-1$ to
$1$, although $\gamma=-1$ is ill-defined for $\rho$), and we treat in our
model the diagonal elements of the connection matrix negligible compared to
the contribution from off-diagonal elements of the connection matrix, because
of large $N$ limit. To arrive at the last step, we have to introduce the
following order parameters

$\displaystyle Q^{ab}=$
$\displaystyle\frac{1}{N}\sum_{i}\phi(x_{i}^{a})\phi(x_{i}^{b}),$ (21)
$\displaystyle R^{ab}=$
$\displaystyle\frac{1}{N}\sum_{i}\hat{x}_{i}^{a}\phi(x_{i}^{b}).$ (22)

Hence, we need to insert the following identity in Eq. (19),

$\displaystyle 1=$ $\displaystyle\prod_{a\leq b}\int
dQ^{ab}\delta\quantity(Q^{ab}-\frac{1}{N}\sum_{i}\phi(x_{i}^{a})\phi(x_{i}^{b}))\prod_{ab}\int
dR^{ab}\delta\quantity(R^{ab}-\frac{1}{N}\sum_{i}\hat{x}_{i}^{a}\phi(x_{i}^{b}))$
(23) $\displaystyle=$
$\displaystyle\int\frac{d\mathbf{Q}d\mathbf{\hat{Q}}d\mathbf{R}d\mathbf{\hat{R}}}{(2\pi)^{D}}\exp\quantity[-\mathrm{i}\sum_{a\leq
b}Q^{ab}\hat{Q}^{ab}-\mathrm{i}\sum_{ab}R^{ab}\hat{R}^{ab}+\mathrm{i}\frac{1}{N}\sum_{i}\sum_{a\leq
b}\hat{Q}^{ab}\phi(x_{i}^{a})\phi(x_{i}^{b})+\mathrm{i}\frac{1}{N}\sum_{ab}\hat{R}^{ab}\sum_{i}\hat{x}_{i}^{a}\phi(x_{i}^{b})]$
$\displaystyle=$
$\displaystyle\int\frac{d\mathbf{Q}d\mathbf{\hat{Q}}d\mathbf{R}d\mathbf{\hat{R}}}{[2\pi(\mathrm{i}/N)^{2}]^{D}}\exp\quantity[-N\sum_{a\leq
b}Q^{ab}\hat{Q}^{ab}-N\sum_{ab}R^{ab}\hat{R}^{ab}+\sum_{i}\sum_{a\leq
b}\hat{Q}^{ab}\phi(x_{i}^{a})\phi(x_{i}^{b})+\mathrm{i}\sum_{ab}\hat{R}^{ab}\sum_{i}\hat{x}_{i}^{a}\phi(x_{i}^{b})],$

where $D=n^{2}+n(n+1)/2$, and we apply the Fourier transformation of Dirac
delta function by introducing the conjugated order parameters $\hat{Q}^{ab}$,
$\hat{R}^{ab}$, and rescale the order parameters as
$\hat{Q}^{ab}\to-\mathrm{i}N\hat{Q}^{ab}$, $R^{ab}\to-\mathrm{i}R^{ab}$ and
$\hat{R}^{ab}\to N\hat{R}^{ab}$ in the last step.
$d\mathbf{Q}d\mathbf{\hat{Q}}d\mathbf{R}d\mathbf{\hat{R}}$ is a shorthand for
$\prod_{a\leq b}dQ^{ab}d\hat{Q}_{ab}\prod_{ab}dR^{ab}d\hat{R}^{ab}$.

As a result, Eq. (19) becomes

$\displaystyle\left\langle Z^{n}\right\rangle\propto$ $\displaystyle\int
d\mathbf{x}D\mathbf{\hat{x}}d\mathbf{Q}d\mathbf{\hat{Q}}d\mathbf{R}d\mathbf{\hat{R}}\exp\left[-\mathrm{i}\sqrt{\beta}\sum_{ia}x_{i}^{a}\hat{x}_{i}^{a}-\beta\eta\sum_{ia}\pqty{x_{i}^{a}}^{2}-\frac{1}{2}g^{2}\beta\sum_{ab}Q^{ab}\pqty{\sum_{i}\hat{x}_{i}^{a}\hat{x}_{i}^{b}}+\frac{1}{2}\beta
g^{2}N\gamma\sum_{ab}R^{ab}R^{ba}\right.$ (24)
$\displaystyle\left.-N\sum_{a\leq
b}Q^{ab}\hat{Q}^{ab}-N\sum_{ab}R^{ab}\hat{R}^{ab}+\sum_{i}\sum_{a\leq
b}\hat{Q}^{ab}\phi(x_{i}^{a})\phi(x_{i}^{b})+\mathrm{i}\sum_{ab}\hat{R}^{ab}\sum_{i}\hat{x}_{i}^{a}\phi(x_{i}^{b})\right],$

where we have neglected the irrelevant constant due to the variable change in
the integral, which does not affect the later saddle point approximation in
the large $N$ limit.

Equation (24) can be rewritten into a compact form,

$\left\langle Z^{n}\right\rangle\propto\int
d\mathbf{Q}d\mathbf{\hat{Q}}d\mathbf{R}d\mathbf{\hat{R}}\exp\bqty{N\quantity(-\sum_{a\leq
b}Q^{ab}\hat{Q}^{ab}-\sum_{ab}R^{ab}\hat{R}^{ab}+\mathcal{G})},$ (25)

where the model-dependent action $\mathcal{G}$ reads,

$\displaystyle\mathcal{G}=$ $\displaystyle\ln\int
d\mathbf{x}D\mathbf{\hat{x}}\exp\left[-\mathrm{i}\sqrt{\beta}\sum_{a}x^{a}\hat{x}^{a}-\beta\eta\sum_{a}\pqty{x^{a}}^{2}-\frac{1}{2}g^{2}\beta\sum_{ab}Q^{ab}\hat{x}^{a}\hat{x}^{b}+\sum_{a\leq
b}\hat{Q}^{ab}\phi(x^{a})\phi(x^{b})\right.$ (26)
$\displaystyle\left.+\frac{1}{2}\beta
g^{2}\gamma\sum_{ab}R^{ab}R^{ba}+\mathrm{i}\sum_{ab}\hat{R}^{ab}\hat{x}^{a}\phi(x^{b})\right],$

where $d\mathbf{x}$ and $D\mathbf{\hat{x}}$ have been reduced to
$\prod_{a}dx^{a}$ and $\prod_{a}D\hat{x}^{a}$, respectively.

To proceed, we have to assume the replica symmetric (RS) ansätz as the first
level of approximation, i.e., the replica overlap matrix are permutation
invariant, as is intuitively expected from the replicated partition function.
This first level of approximation can be refined by going to multiple steps of
replica symmetry breaking [20], provided that the theoretical prediction is in
significant disagreement with numerical simulations. The RS ansätz reads
$Q^{ab}=q\delta_{ab}+Q(1-\delta_{ab})$,
$\hat{Q}^{ab}=\hat{q}\delta_{ab}+\hat{Q}(1-\delta_{ab})$,
$R^{ab}=r\delta_{ab}+R(1-\delta_{ab})$ and
$\hat{R}^{ab}=\hat{r}\delta_{ab}+\hat{R}(1-\delta_{ab})$. We have then the
following result of the averaged replicated partition function,

$\left\langle Z^{n}\right\rangle\propto\int
dQd\hat{Q}dqd\hat{q}dRd\hat{R}drd\hat{r}\exp\bqty{N\pqty{-\frac{n(n-1)}{2}Q\hat{Q}-nq\hat{q}-n(n-1)R\hat{R}-nr\hat{r}+\mathcal{G}}},$
(27)

where

$\displaystyle\mathcal{G}=$ $\displaystyle\ln\int
d\mathbf{x}D\mathbf{\hat{x}}\exp\left[-\mathrm{i}\sqrt{\beta}\sum_{a}x^{a}\hat{x}^{a}-\beta\eta\sum_{a}(x^{a})^{2}+\frac{1}{2}\beta
g^{2}\gamma\quantity(n(n-1)R^{2}+nr^{2})-\frac{1}{2}g^{2}\beta\pqty{Q\pqty{\sum_{a}\hat{x}^{a}}^{2}+\pqty{q-Q}\sum_{a}\pqty{\hat{x}^{a}}^{2}}\right.$
(28)
$\displaystyle\left.+\frac{1}{2}\pqty{\hat{Q}\pqty{\sum_{a}\phi(x^{a})}^{2}+\pqty{2\hat{q}-\hat{Q}}\sum_{a}\phi^{2}(x^{a})}+\mathrm{i}\quantity(\hat{R}\sum_{ab}\hat{x}^{a}\phi(x^{b})+\quantity(\hat{r}-\hat{R})\sum_{a}\hat{x}^{a}\phi(x^{a}))\right].$

To linearize the quadratic terms in Eq. (28), we use the HS transformation
once again as follows,

$e^{\frac{1}{2}\hat{Q}\pqty{\sum_{a}\phi(x^{a})}^{2}-\frac{1}{2}g^{2}\beta
Q\pqty{\sum_{a}\hat{x}^{a}}^{2}+\mathrm{i}\hat{R}\sum_{ab}\hat{x}^{a}\phi(x^{b})}=\int
DuDve^{\kappa_{1}u\sum_{a}\phi(x^{a})+v\quantity(\kappa_{2}\sum_{a}\phi(x^{a})+\mathrm{i}\kappa_{3}\sum_{a}\hat{x}^{a})},$
(29)

where the auxiliary coefficients $\kappa_{1}=\sqrt{\frac{g^{2}\beta
Q\hat{Q}-\hat{R}^{2}}{g^{2}\beta Q}}$, $\kappa_{2}=\frac{\hat{R}}{g\sqrt{\beta
Q}}$ and $\kappa_{3}=g\sqrt{\beta Q}$.

Therefore, we have

$\mathcal{G}=\ln\int DuDv\exp\left[\frac{1}{2}\beta
g^{2}\gamma\quantity(n(n-1)R^{2}+nr^{2})\right]I^{n},$ (30)

where

$\displaystyle I$ $\displaystyle\equiv\int dxD\hat{x}\exp\left[-\beta\eta
x^{2}-\mathrm{i}\sqrt{\beta}x\hat{x}-\frac{1}{2}g^{2}\beta\pqty{q-Q}\hat{x}^{2}+\frac{1}{2}\pqty{2\hat{q}-\hat{Q}}\phi^{2}(x)+\mathrm{i}\quantity(\hat{r}-\hat{R})\hat{x}\phi(x)\right.$
(31)
$\displaystyle\left.+\mathrm{i}\kappa_{3}v\hat{x}+\quantity(\kappa_{1}u+\kappa_{2}v)\phi(x)\right].$

We next complete the integral over $D\hat{x}$ as

$\int
D\hat{x}\exp\bqty{-\frac{1}{2}g^{2}\beta\pqty{q-Q}\hat{x}^{2}+i\quantity(\kappa_{3}v+\quantity(\hat{r}-\hat{R})-\sqrt{\beta}x)\hat{x}}=\frac{1}{\sigma}\exp\bqty{-\frac{1}{2}\frac{1}{\sigma^{2}}\quantity(\kappa_{3}v+(\hat{r}-\hat{R})\phi(x)-\sqrt{\beta}x)^{2}},$
(32)

where $\sigma\equiv\sqrt{1+g^{2}\beta(q-Q)}$. Thus we obtain,

$\mathcal{G}=-n\ln\sigma+\frac{1}{2}\beta
g^{2}\gamma\quantity(n(n-1)R^{2}+nr^{2})+\ln\int DuDvI^{n},$ (33)

where $I\equiv\int dxe^{\mathcal{H}(x)}$, and

$\mathcal{H}(x)\equiv-\beta\eta
x^{2}+\frac{1}{2}\pqty{2\hat{q}-\hat{Q}}\phi^{2}(x)+\quantity(\kappa_{1}u+\kappa_{2}v)\phi(x)-\frac{1}{2}\frac{1}{\sigma^{2}}\quantity(\kappa_{3}v+(\hat{r}-\hat{R})\phi(x)-\sqrt{\beta}x)^{2},$
(34)

where $\mathcal{H}(x)$ could be considered as the single-variable effective
Hamiltonian of our model, which emerges after the quenched disorder is
averaged out.

Finally, we perform the saddle point approximation because of large $N$, such
that $\left\langle Z^{n}\right\rangle\propto\int
d\mathcal{O}d\hat{\mathcal{O}}e^{N\Phi(\mathcal{O},\hat{\mathcal{O}})}\approx
e^{N\Phi(\mathcal{O}^{*},\hat{\mathcal{O}}^{*})}$, where
$\mathcal{O}=\\{Q,q,R,r\\}$ and the superscript $*$ refers to the saddle point
value. The free energy can be obtained from the replica trick $-\beta
f=\lim_{n\to 0}\frac{\ln\left\langle Z^{n}\right\rangle}{nN}$, and therefore,

$-\beta
f=\frac{1}{2}Q\hat{Q}-q\hat{q}+R\hat{R}-r\hat{r}-\ln\sigma+\frac{1}{2}\beta
g^{2}\gamma\quantity(r^{2}-R^{2})+\int DuDv\ln I.$ (35)

##  Appendix B Saddle Point Equations

In this section, we omit the superscript $*$ on the order parameters, and
derive the self-consistent equations these order parameters must obey, which
is the so-called saddle point equations (SDEs). The principle is to optimize
the action with respect to the order parameters and their conjugated
counterparts. The result is given bellow,

$\displaystyle q=$
$\displaystyle\quantity[\left\langle\phi^{2}\right\rangle],$ (36a)
$\displaystyle Q=$
$\displaystyle\quantity[\left\langle\phi\right\rangle^{2}],$ (36b)
$\displaystyle r=$
$\displaystyle-\frac{1}{\sigma^{2}}f(1,0,-1)\quantity[\left\langle\phi^{2}\right\rangle]+\frac{1}{\sigma^{2}}f(0,1,-1)\quantity[\left\langle\phi\right\rangle^{2}]+\frac{\sqrt{\beta}}{\sigma^{2}}\quantity(1-gkQ)\quantity[\left\langle
x\phi\right\rangle]+\frac{\sqrt{\beta}}{\sigma^{2}}gkQ\quantity[\left\langle
x\right\rangle\left\langle\phi\right\rangle],$ (36c) $\displaystyle R=$
$\displaystyle-\frac{1}{\sigma^{2}}f(0,1,-1)\quantity[\left\langle\phi^{2}\right\rangle]-\frac{1}{\sigma^{2}}f(1,-2,1)\quantity[\left\langle\phi\right\rangle^{2}]-\frac{\sqrt{\beta}}{\sigma^{2}}gkQ\quantity[\left\langle
x\phi\right\rangle]+\frac{\sqrt{\beta}}{\sigma^{2}}\quantity(1+gkQ)\quantity[\left\langle
x\right\rangle\left\langle\phi\right\rangle],$ (36d) $\displaystyle\hat{q}=$
$\displaystyle-\frac{gk}{2}+\frac{g^{2}k^{2}Q}{2}+\frac{gk}{2\sigma^{2}}\quantity(\hat{r}-\hat{R})f(1,1,-2)\quantity[\left\langle\phi^{2}\right\rangle]+\frac{gk}{\sigma^{2}}\quantity(\hat{r}-\hat{R})f(0,-1,1)\quantity[\left\langle\phi\right\rangle^{2}]+\frac{k^{2}}{2}\quantity(1-2gkQ)\quantity[\left\langle
x^{2}\right\rangle]$ (36e) $\displaystyle+$
$\displaystyle\frac{gk\sqrt{\beta}}{\sigma^{2}}f(0,1,-2)\quantity[\left\langle
x\right\rangle\left\langle\phi\right\rangle]+\frac{gk\sqrt{\beta}}{\sigma^{2}}f(-1,0,2)\quantity[\left\langle
x\phi\right\rangle]+gk^{3}Q\quantity[\left\langle x\right\rangle^{2}],$ (36f)
$\displaystyle\hat{Q}=$ $\displaystyle
g^{2}k^{2}Q+\frac{2gk}{\sigma^{2}}\quantity(\hat{r}-\hat{R})f(0,1,-1)\quantity[\left\langle\phi^{2}\right\rangle]+\frac{gk}{\sigma^{2}}\quantity(\hat{r}-\hat{R})f(1,-3,2)\quantity[\left\langle\phi\right\rangle^{2}]-2gk^{3}Q\quantity[\left\langle
x^{2}\right\rangle]$ (36g) $\displaystyle-$
$\displaystyle\frac{2gk\sqrt{\beta}}{\sigma^{2}}f(1,-2,2)\quantity[\left\langle
x\right\rangle\left\langle\phi\right\rangle]+\frac{2gk\sqrt{\beta}}{\sigma^{2}}f(0,-1,2)\quantity[\left\langle
x\phi\right\rangle]+k^{2}\quantity(1+2gkQ)\quantity[\left\langle
x\right\rangle^{2}],$ (36h) $\displaystyle\hat{r}=$ $\displaystyle\beta
g^{2}\gamma r,$ (36i) $\displaystyle\hat{R}=$ $\displaystyle\beta g^{2}\gamma
R,$ (36j)

where $k\equiv\frac{g\beta}{\sigma^{2}}$ and $f(a,b,c)\equiv
a\hat{r}+b\hat{R}+cgkQ\quantity(\hat{r}-\hat{R})$, and we define disorder and
thermal averages, respectively as follows,

$\quantity[\bullet]\equiv\int DuDv\bullet,$ (37)

and

$\left\langle\bullet\right\rangle\equiv\frac{\int
dxe^{\mathcal{H}(x)}\bullet}{\int dxe^{\mathcal{H}(x)}}.$ (38)

Note that we use the integral identity $\int Dxxf(x)=\int Dxf^{\prime}(x)$ to
write the above SDEs in the compact form.

We are also interested in the average energy, whose ground state value is
related to the speed level of the original dynamics. The thermodynamic
relation $\left\langle E\right\rangle=\frac{\partial[\beta f]}{\partial\beta}$
is used to derive the following formula,

$\displaystyle U=$
$\displaystyle\frac{g^{2}}{2\sigma^{2}}\quantity(q-gkQ(q-Q))-\frac{1}{2}g^{2}(r^{2}-R^{2})\gamma+\frac{1}{2\sigma^{2}}\quantity(1+2\eta\sigma^{2}-gk(q-Q)-\frac{2gkQ}{\sigma^{2}})\quantity[\left\langle
x^{2}\right\rangle]+\frac{gkQ}{\sigma^{4}}\quantity[\left\langle
x\right\rangle^{2}]$ (39)
$\displaystyle-\frac{1}{2\sqrt{\beta}\sigma^{2}}\quantity(f(1,-1,0)+\frac{1}{\sigma^{2}}f(0,1,-3)+gk(q-Q)f(-2,1,1))\quantity[\left\langle
x\phi\right\rangle]$
$\displaystyle-\frac{1}{2\sqrt{\beta}\sigma^{2}}\quantity(\frac{1}{\sigma^{2}}f(0,-1,3)+gk(q-Q)f(0,1,-1))\quantity[\left\langle
x\right\rangle\left\langle\phi\right\rangle]$
$\displaystyle-\frac{1}{2g^{2}Q\beta^{2}}\quantity(\hat{R}^{2}+g^{2}k^{2}Q(q-Q)\quantity(\hat{r}-\hat{R})^{2}+f(0,-1,1)\quantity(f(0,1,1)-2g^{2}k^{2}Q(q-Q)(\hat{r}-\hat{R})))\quantity[\left\langle\phi^{2}\right\rangle]$
$\displaystyle+\frac{1}{2g^{2}Q\beta^{2}}\quantity(\hat{R}^{2}+f(0,-1,1)\quantity(f(0,1,1)-2g^{2}k^{2}Q(q-Q)(\hat{r}-\hat{R})))\quantity[\left\langle\phi\right\rangle^{2}].$

To characterize the average activity level in the neural population, we also
derive the $\ell_{2}$ norm of the activity as follows

$\frac{1}{N}\left\langle\|\mathbf{x}\|^{2}\right\rangle=-\frac{1}{\beta}\partialderivative{\quantity(-\beta
f)}{\eta}=\quantity[\left\langle x^{2}\right\rangle].$ (40)

##  Appendix C Independent case $\gamma=0$

In the case $\gamma=0$, all the connection weights are independent and
uncorrelated. Hence, the expression of the SDEs gets greatly simplified as

$\displaystyle q=$
$\displaystyle\quantity[\left\langle\phi^{2}\right\rangle],$ (41)
$\displaystyle Q=$
$\displaystyle\quantity[\left\langle\phi\right\rangle^{2}],$ (42)
$\displaystyle r=$
$\displaystyle\frac{\sqrt{\beta}}{\sigma^{2}}\quantity(1-gkQ)\quantity[\left\langle
x\phi\right\rangle]+\frac{\sqrt{\beta}}{\sigma^{2}}gkQ\quantity[\left\langle
x\right\rangle\left\langle\phi\right\rangle],$ (43) $\displaystyle R=$
$\displaystyle-\frac{\sqrt{\beta}}{\sigma^{2}}gkQ\quantity[\left\langle
x\phi\right\rangle]+\frac{\sqrt{\beta}}{\sigma^{2}}\quantity(1+gkQ)\quantity[\left\langle
x\right\rangle\left\langle\phi\right\rangle],$ (44) $\displaystyle\hat{q}=$
$\displaystyle-\frac{gk}{2}+\frac{g^{2}k^{2}Q}{2}+\frac{k^{2}}{2}\quantity(1-2gkQ)\quantity[\left\langle
x^{2}\right\rangle]+gk^{3}Q\quantity[\left\langle x\right\rangle^{2}],$ (45)
$\displaystyle\hat{Q}=$ $\displaystyle
g^{2}k^{2}Q-2gk^{3}Q\quantity[\left\langle
x^{2}\right\rangle]+k^{2}\quantity(1+2gkQ)\quantity[\left\langle
x\right\rangle^{2}],$ (46) $\displaystyle\hat{r}=$ $\displaystyle 0,$ (47)
$\displaystyle\hat{R}=$ $\displaystyle 0.$ (48)

Because $\hat{r}=\hat{R}=0$, we have $f(a,b,c)=0$. The free energy is thus
reduced to

$-\beta
f=\frac{1}{2}Q\hat{Q}-q\hat{q}-\ln\sigma-\frac{g^{2}\beta}{2\sigma^{2}}Q+\int
DuDv\ln I,$ (49)

where

$I\equiv\int
dx\exp\bqty{-\frac{1}{2}\beta\pqty{\frac{1}{\sigma^{2}}+2\eta}x^{2}+\frac{1}{2}\pqty{2\hat{q}-\hat{Q}}\phi^{2}(x)+\sqrt{\hat{Q}}u\phi(x)+\frac{g\beta}{\sigma^{2}}\sqrt{Q}vx}.$
(50)

This free energy does not depend on the order parameters $r$ and $R$ anymore.
The average energy also yields the following simple form

$U=\frac{g^{2}}{2\sigma^{2}}\quantity(q-gkQ(q-Q))+\frac{1}{2\sigma^{2}}\quantity(1+2\eta\sigma^{2}-gk(q-Q)-\frac{2gkQ}{\sigma^{2}})\quantity[\left\langle
x^{2}\right\rangle]+\frac{gkQ}{\sigma^{4}}\quantity[\left\langle
x\right\rangle^{2}].$ (51)

##  Appendix D Additional experiment results

To illustrate the phase transition in more details, we choose some curves
($\gamma=0.2$, $0.5$, $0.8$) and plot them in the 2-dimension plane (Fig. 5).
$q$ and $Q$ are close to zero before the critical point, while growing rapidly
after the critical point. $r$ and $R$ achieve their maximum values at the
critical point. This is also supported by the iteration dynamics of the SDEs
(Fig. 7).

Next we explore the effect of $\beta$ in optimizing the energy and the order
parameters. When $g$ is close to the critical point, we show how $q$ varies
with increasing $\beta$ (see Fig. 6 for $g=0.98$ and $g=1.0$). We observe that
$q$ converges more rapidly towards zero when $g$ is below the critical value.
When $g$ is approaching the critical value from below, the decaying towards
zero becomes much slower, signaling a continuous phase transition.

Figure 5: Phase transitions as $g$ increases. We plot three different curves
($\gamma=0.2$, $0.5$, $0.8$), and different colors represent different
$\gamma$ values. The dashed line represents the DMFT result of critical points
in Ref. [19]. Figure 6: The relationship between $q$ and $\beta$ when $g=0.98$
and $g=1.0$ ($\gamma=0$). The blue dots denote $g=0.98$ case, and the orange
dots denote $g=1.0$. Both axes are in the log scale. Figure 7: The iterative
dynamics of $q$ and $Q$ for different values of $g$. All initial values are
set to $10^{-5}$ for $q$ and $10^{-6}$ for $Q$. The data points are obtained
by using the Algo. 1.

##  Appendix E Numerical details of solving SDEs

The most expensive cost of solving SDEs comes from the double average terms,
such as $\quantity[\left\langle x\right\rangle^{2}]$ and
$\quantity[\left\langle\phi\right\rangle^{2}]$, where the outer integral is
related to $DuDv$ and the inner integral is related $dx$. These integrals are
calculated by Monte Carlo (MC) method with important sampling technique [17].
Therefore, we need to generate $100k\times 6k$ Gaussian samples, where $100k$
MC samples of $u$ and $v$ is used for the outer integral and $6k$ MC samples
of $x$ for each pair of $u$ and $v$, is used for the inner integral. Note that
the inner integral in Eq. (38) does not contain Gaussian measure at the first
sight. However, we can take a reorganization of terms in the effective
Hamiltonian. We have then the following equivalent transformation,

$\left\langle\bullet\right\rangle\equiv\frac{\int
Dxe^{\tilde{\mathcal{H}}(x)}\bullet}{\int Dxe^{\tilde{\mathcal{H}}(x)}},$ (52)

where the Gaussian measure is introduced as $Dx\equiv
e^{-\frac{1}{2}\beta\quantity(\frac{1}{\sigma^{2}}+2\eta)x^{2}}dx/\sqrt{2\pi}$
and

$\tilde{\mathcal{H}}(x)\equiv\frac{1}{2}\pqty{2\hat{q}-\hat{Q}}\phi^{2}(x)+\quantity(\kappa_{1}u+\kappa_{2}v)\phi(x)-\frac{1}{2}\frac{1}{\sigma^{2}}\quantity(\kappa_{3}v+(\hat{r}-\hat{R})\phi(x))^{2}+\frac{\sqrt{\beta}}{\sigma^{2}}\quantity(\kappa_{3}v+(\hat{r}-\hat{R})\phi(x))x.$
(53)

After initialization of the order parameters, we complete the double average
using the above MC sampling. Then we update the order parameters as follows,

$\mathcal{O}_{t+1}=\alpha\mathcal{O}_{t}+(1-\alpha)f(\mathcal{O}_{t}),$ (54)

where $\alpha$ is the damping parameter (a value close to one) to speed up the
convergence. $f(\mathcal{O})$ represents the right-hand-side of Eqs. (36). We
stop the iteration when
$\absolutevalue{\mathcal{O}_{t+1}-\mathcal{O}_{t}}<10^{-3}$ for all order
parameters. We give a pseudocode in Algorithm. 1 for solving the SDEs.

Algorithm 1 Solver of SDEs

1:The initial values of order parameters $\mathcal{O}_{0}$,
$\mathcal{O}\in\\{q,Q,r,R\\}$

2:The convergent values of order parameters $\mathcal{O}^{*}$

3:Generate $100k\times 6k$ Gaussian Monte Carlo samples

4:repeat

5: Calculate all the double averages by the MC method in Eqs. (36)

6: Let
$\mathcal{O}_{t+1}\leftarrow\alpha\mathcal{O}_{t}+(1-\alpha)f(\mathcal{O}_{t})$,
where $f(\mathcal{O})$ are the functions specified by the SDEs [see Eqs. (36)]

7:until $\absolutevalue{\mathcal{O}_{t+1}-\mathcal{O}_{t}}<10^{-3}$

8:Let $\mathcal{O}^{*}\leftarrow\mathcal{O}_{t}$

###### Acknowledgements.

We would like to thank Zijian Jiang, Yang Zhao and Wenkang Du for helpful
discussions at earlier stages of this project. This research was supported by
the National Natural Science Foundation of China for Grant Number 12122515
(H.H.), and Guangdong Provincial Key Laboratory of Magnetoelectric Physics and
Devices (No. 2022B1212010008), and Guangdong Basic and Applied Basic Research
Foundation (Grant No. 2023B1515040023).

## References

  * [1] Hannes Risken.  The Fokker-Planck Equation: Methods of Solution and Applications.  Springer-Verlag Berlin, Berlin, 1996.
  * [2] Steven H Strogatz.  Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering.  CRC Press, 2015.
  * [3] Mikhail I. Rabinovich, Pablo Varona, Allen I. Selverston, and Henry D. I. Abarbanel.  Dynamical principles in neuroscience.  Rev. Mod. Phys., 78:1213–1265, 2006.
  * [4] Jorge Kurchan.  Fluctuation theorem for stochastic dynamics.  Journal of Physics A: Mathematical and General, 31(16):3719, 1998\.
  * [5] Takahiro Ohkuma and Takao Ohta.  Fluctuation theorems for non-linear generalized langevin systems.  Journal of Statistical Mechanics: Theory and Experiment, 2007(10):P10010, 2007.
  * [6] Udo Seifert.  Stochastic thermodynamics, fluctuation theorems and molecular machines.  Reports on Progress in Physics, 75(12):126001, 2012.
  * [7] P. C. Martin, E. D. Siggia, and H. A. Rose.  Statistical dynamics of classical systems.  Phys. Rev. A, 8:423–437, 1973.
  * [8] Hans-Karl Janssen.  On a lagrangean for classical field dynamics and renormalization group calculations of dynamical critical properties.  Zeitschrift für Physik B Condensed Matter, 23(4):377–380, 1976\.
  * [9] C. De Dominicis.  Dynamics as a substitute for replicas in systems with quenched random impurities.  Phys. Rev. B, 18:4913–4919, 1978.
  * [10] Hans-Jürgen Sommers.  Path-integral approach to ising spin-glass dynamics.  Phys. Rev. Lett., 58:1268–1271, 1987.
  * [11] Takahiro Hatano and Shin-ichi Sasa.  Steady-state thermodynamics of langevin systems.  Phys. Rev. Lett., 86:3463–3466, 2001.
  * [12] A. Crisanti and H. Sompolinsky.  Dynamics of spin systems with randomly asymmetric bonds: Langevin dynamics and a spherical model.  Phys. Rev. A, 36:4922–4939, 1987.
  * [13] Carson C. Chow and Michael A. Buice.  Path integral methods for stochastic differential equations.  The Journal of Mathematical Neuroscience (JMN), 5(1):8, 2015.
  * [14] John A Hertz, Yasser Roudi, and Peter Sollich.  Path integral methods for the dynamics of stochastic and disordered systems.  Journal of Physics A: Mathematical and Theoretical, 50(3):033001, 2017.
  * [15] Wenxuan Zou and Haiping Huang.  Introduction to dynamical mean-field theory of randomly connected neural networks with bidirectionally correlated couplings.  arXiv:2305.08459, 2023.
  * [16] David Sussillo and Omri Barak.  Opening the black box: Low-dimensional dynamics in high-dimensional recurrent neural networks.  Neural Computation, 25:626–649, 2013.
  * [17] Haiping Huang.  Statistical Mechanics of Neural Networks.  Springer, Singapore, 2022.
  * [18] H. Sompolinsky, A. Crisanti, and H. J. Sommers.  Chaos in random neural networks.  Phys. Rev. Lett., 61:259–262, 1988.
  * [19] Daniel Marti, Nicolas Brunel, and Srdjan Ostojic.  Correlations between synapses in pairs of neurons slow down dynamics in randomly connected neural networks.  Phys. Rev. E, 97:062314, 2018.
  * [20] M. Mézard, G. Parisi, and M. A. Virasoro.  Spin Glass Theory and Beyond.  World Scientific, Singapore, 1987.
  * [21] Hugo Touchette.  The large deviation approach to statistical mechanics.  Physics Reports, 478(1-3):1–69, 2009.
  * [22] Goel NS, Maitra SC, and Montroll EW.  On the volterra and other nonlinear models of interacting populations.  Rev. Mod. Phys., 43:231–276, 1971.
  * [23] Robert May.  Will a large complex system be stable?  Nature, 238(5364):413–414, 1972.
  * [24] Yan V. Fyodorov and Boris A. Khoruzhenko.  Nonlinear analogue of the may-wigner instability transition.  Proceedings of the National Academy of Sciences, 113:6827–6832, 2016.
  * [25] Bertrand Lacroix-A-Chez-Toine and Yan V Fyodorov.  Counting equilibria in a random non-gradient dynamics with heterogeneous relaxation rates.  Journal of Physics A: Mathematical and Theoretical, 55(14):144001, 2022.
  * [26] Pierpaolo Vivo.  Random linear systems with quadratic constraints: from random matrix theory to replicas and back.  arXiv:2401.03209, 2024.
  * [27] Max Welling and Yee Whye Teh.  Bayesian learning via stochastic gradient langevin dynamics.  In International Conference on Machine Learning, 2011.
  * [28] Wulfram Gerstner, Werner M. Kistler, Richard Naud, and Liam Paninski.  Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition.  Cambridge University Press, United Kingdom, 2014.
  * [29] Jannis Schuecker, Sven Goedeke, and Moritz Helias.  Optimal sequence memory in driven random networks.  Phys. Rev. X, 8:041029, Nov 2018.
  * [30] M. Helias and D. Dahmen.  Statistical field theory for neural networks.  Springer, Berlin, 2020.
  * [31] David G. Clark, L. F. Abbott, and Ashok Litwin-Kumar.  Dimension of activity in random neural networks.  Phys. Rev. Lett., 131:118401, 2023.
  * [32] Zijian Jiang, Ziming Chen, Tianqi Hou, and Haiping Huang.  Spectrum of non-hermitian deep-hebbian neural networks.  Phys. Rev. Res., 5:013090, 2023.
  * [33] Shun-Ichi Amari.  Characteristics of random nets of analog neuron-like elements.  IEEE Transactions on Systems, Man, and Cybernetics, SMC-2(5):643–657, 1972.
  * [34] Gilles Wainrib and Jonathan Touboul.  Topological and dynamical complexity of random neural networks.  Phys. Rev. Lett., 110:118101, 2013.
  * [35] H. J. Sommers, A. Crisanti, H. Sompolinsky, and Y. Stein.  Spectrum of large random asymmetric matrices.  Phys. Rev. Lett., 60:1895–1898, 1988.
  * [36] Han Yan, Lei Zhao, Liang Hu, Xidi Wang, Erkang Wang, and Jin Wang.  Nonequilibrium landscape theory of neural networks.  Proceedings of the National Academy of Sciences of the United States of America, 110(45):E4185–E4194, 2013.
  * [37] ACC Coolen.  Statistical mechanics of recurrent neural networks i—statics.  In Handbook of biological physics, volume 4, pages 553–618. Elsevier, 2001.
  * [38] Chris G. Langton.  Computation at the edge of chaos: Phase transitions and emergent computation.  Physica D: Nonlinear Phenomena, 42(1):12–37, 1990.
  * [39] Nils Bertschinger and Thomas Natschläger.  Real-time computation at the edge of chaos in recurrent neural networks.  Neural computation, 16(7):1413–1436, 2004.
  * [40] Daniel Toker, Ioannis Pappas, Janna D Lendner, Joel Frohlich, Diego M Mateos, Suresh Muthukumaraswamy, Robin Carhart-Harris, Michelle Paff, Paul M Vespa, Martin M Monti, et al.  Consciousness is supported by near-critical slow cortical electrodynamics.  Proceedings of the National Academy of Sciences, 119(7):e2024455119, 2022.
  * [41] G. Parisi and N. Sourlas.  Supersymmetric field theories and stochastic differential equations.  Nuclear Physics B, 206(2):321–332, 1982.
  * [42] Haiping Huang and Yoshiyuki Kabashima.  Origin of the computational hardness for learning with binary synapses.  Phys. Rev. E, 90:052813, 2014.
  * [43] Haiping Huang.  Eight challenges in developing theory of intelligence.  arXiv:2306.11232, 2023.